mesh_io_msh.cc
No OneTemporary
Actions

Subscribers

None

File Metadata

Created: Sun, May 26, 23:06

mesh_io_msh.cc
View Options

	/**
	* @file mesh_io_msh.cc
	*
	* @author Dana Christen <dana.christen@gmail.com>
	* @author Mauro Corrado <mauro.corrado@epfl.ch>
	* @author David Simon Kammer <david.kammer@epfl.ch>
	* @author Nicolas Richart <nicolas.richart@epfl.ch>
	*
	* @date creation: Fri Jun 18 2010
	* @date last modification: Tue Feb 20 2018
	*
	* @brief Read/Write for MSH files generated by gmsh
	*
	* @section LICENSE
	*
	* Copyright (©) 2010-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
	* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
	*
	* Akantu is free software: you can redistribute it and/or modify it under the
	* terms of the GNU Lesser General Public License as published by the Free
	* Software Foundation, either version 3 of the License, or (at your option) any
	* later version.
	*
	* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
	* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
	* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
	* details.
	*
	* You should have received a copy of the GNU Lesser General Public License
	* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
	*
	*/

	/* -----------------------------------------------------------------------------
	Version (Legacy) 1.0

	$NOD
	number-of-nodes
	node-number x-coord y-coord z-coord
	...
	$ENDNOD
	$ELM
	number-of-elements
	elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list
	...
	$ENDELM
	-----------------------------------------------------------------------------
	Version 2.1

	$MeshFormat
	version-number file-type data-size
	$EndMeshFormat
	$Nodes
	number-of-nodes
	node-number x-coord y-coord z-coord
	...
	$EndNodes
	$Elements
	number-of-elements
	elm-number elm-type number-of-tags < tag > ... node-number-list
	...
	$EndElements
	$PhysicalNames
	number-of-names
	physical-dimension physical-number "physical-name"
	...
	$EndPhysicalNames
	$NodeData
	number-of-string-tags
	< "string-tag" >
	...
	number-of-real-tags
	< real-tag >
	...
	number-of-integer-tags
	< integer-tag >
	...
	node-number value ...
	...
	$EndNodeData
	$ElementData
	number-of-string-tags
	< "string-tag" >
	...
	number-of-real-tags
	< real-tag >
	...
	number-of-integer-tags
	< integer-tag >
	...
	elm-number value ...
	...
	$EndElementData
	$ElementNodeData
	number-of-string-tags
	< "string-tag" >
	...
	number-of-real-tags
	< real-tag >
	...
	number-of-integer-tags
	< integer-tag >
	...
	elm-number number-of-nodes-per-element value ...
	...
	$ElementEndNodeData

	-----------------------------------------------------------------------------
	elem-type

	1: 2-node line.
	2: 3-node triangle.
	3: 4-node quadrangle.
	4: 4-node tetrahedron.
	5: 8-node hexahedron.
	6: 6-node prism.
	7: 5-node pyramid.
	8: 3-node second order line
	9: 6-node second order triangle
	10: 9-node second order quadrangle
	11: 10-node second order tetrahedron
	12: 27-node second order hexahedron
	13: 18-node second order prism
	14: 14-node second order pyramid
	15: 1-node point.
	16: 8-node second order quadrangle
	17: 20-node second order hexahedron
	18: 15-node second order prism
	19: 13-node second order pyramid
	20: 9-node third order incomplete triangle
	21: 10-node third order triangle
	22: 12-node fourth order incomplete triangle
	23: 15-node fourth order triangle
	24: 15-node fifth order incomplete triangle
	25: 21-node fifth order complete triangle
	26: 4-node third order edge
	27: 5-node fourth order edge
	28: 6-node fifth order edge
	29: 20-node third order tetrahedron
	30: 35-node fourth order tetrahedron
	31: 56-node fifth order tetrahedron
	-------------------------------------------------------------------------- */

	/* -------------------------------------------------------------------------- */
	#include "element_group.hh"
	#include "mesh_io.hh"
	#include "mesh_utils.hh"
	#include "node_group.hh"
	/* -------------------------------------------------------------------------- */
	#include <fstream>
	/* -------------------------------------------------------------------------- */

	namespace akantu {

	/* -------------------------------------------------------------------------- */
	/* Methods Implentations */
	/* -------------------------------------------------------------------------- */
	MeshIOMSH::MeshIOMSH() {
	canReadSurface = true;
	canReadExtendedData = true;

	_msh_nodes_per_elem[_msh_not_defined] = 0;
	_msh_nodes_per_elem[_msh_segment_2] = 2;
	_msh_nodes_per_elem[_msh_triangle_3] = 3;
	_msh_nodes_per_elem[_msh_quadrangle_4] = 4;
	_msh_nodes_per_elem[_msh_tetrahedron_4] = 4;
	_msh_nodes_per_elem[_msh_hexahedron_8] = 8;
	_msh_nodes_per_elem[_msh_prism_1] = 6;
	_msh_nodes_per_elem[_msh_pyramid_1] = 1;
	_msh_nodes_per_elem[_msh_segment_3] = 3;
	_msh_nodes_per_elem[_msh_triangle_6] = 6;
	_msh_nodes_per_elem[_msh_quadrangle_9] = 9;
	_msh_nodes_per_elem[_msh_tetrahedron_10] = 10;
	_msh_nodes_per_elem[_msh_hexahedron_27] = 27;
	_msh_nodes_per_elem[_msh_hexahedron_20] = 20;
	_msh_nodes_per_elem[_msh_prism_18] = 18;
	_msh_nodes_per_elem[_msh_prism_15] = 15;
	_msh_nodes_per_elem[_msh_pyramid_14] = 14;
	_msh_nodes_per_elem[_msh_point] = 1;
	_msh_nodes_per_elem[_msh_quadrangle_8] = 8;

	_msh_to_akantu_element_types[_msh_not_defined] = _not_defined;
	_msh_to_akantu_element_types[_msh_segment_2] = _segment_2;
	_msh_to_akantu_element_types[_msh_triangle_3] = _triangle_3;
	_msh_to_akantu_element_types[_msh_quadrangle_4] = _quadrangle_4;
	_msh_to_akantu_element_types[_msh_tetrahedron_4] = _tetrahedron_4;
	_msh_to_akantu_element_types[_msh_hexahedron_8] = _hexahedron_8;
	_msh_to_akantu_element_types[_msh_prism_1] = _pentahedron_6;
	_msh_to_akantu_element_types[_msh_pyramid_1] = _not_defined;
	_msh_to_akantu_element_types[_msh_segment_3] = _segment_3;
	_msh_to_akantu_element_types[_msh_triangle_6] = _triangle_6;
	_msh_to_akantu_element_types[_msh_quadrangle_9] = _not_defined;
	_msh_to_akantu_element_types[_msh_tetrahedron_10] = _tetrahedron_10;
	_msh_to_akantu_element_types[_msh_hexahedron_27] = _not_defined;
	_msh_to_akantu_element_types[_msh_hexahedron_20] = _hexahedron_20;
	_msh_to_akantu_element_types[_msh_prism_18] = _not_defined;
	_msh_to_akantu_element_types[_msh_prism_15] = _pentahedron_15;
	_msh_to_akantu_element_types[_msh_pyramid_14] = _not_defined;
	_msh_to_akantu_element_types[_msh_point] = _point_1;
	_msh_to_akantu_element_types[_msh_quadrangle_8] = _quadrangle_8;

	_akantu_to_msh_element_types[_not_defined] = _msh_not_defined;
	_akantu_to_msh_element_types[_segment_2] = _msh_segment_2;
	_akantu_to_msh_element_types[_segment_3] = _msh_segment_3;
	_akantu_to_msh_element_types[_triangle_3] = _msh_triangle_3;
	_akantu_to_msh_element_types[_triangle_6] = _msh_triangle_6;
	_akantu_to_msh_element_types[_tetrahedron_4] = _msh_tetrahedron_4;
	_akantu_to_msh_element_types[_tetrahedron_10] = _msh_tetrahedron_10;
	_akantu_to_msh_element_types[_quadrangle_4] = _msh_quadrangle_4;
	_akantu_to_msh_element_types[_quadrangle_8] = _msh_quadrangle_8;
	_akantu_to_msh_element_types[_hexahedron_8] = _msh_hexahedron_8;
	_akantu_to_msh_element_types[_hexahedron_20] = _msh_hexahedron_20;
	_akantu_to_msh_element_types[_pentahedron_6] = _msh_prism_1;
	_akantu_to_msh_element_types[_pentahedron_15] = _msh_prism_15;
	_akantu_to_msh_element_types[_point_1] = _msh_point;

	#if defined(AKANTU_STRUCTURAL_MECHANICS)
	_akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2;
	_akantu_to_msh_element_types[_bernoulli_beam_3] = _msh_segment_2;
	_akantu_to_msh_element_types[_discrete_kirchhoff_triangle_18] =
	_msh_triangle_3;
	#endif

	std::map<ElementType, MSHElementType>::iterator it;
	for (it = _akantu_to_msh_element_types.begin();
	it != _akantu_to_msh_element_types.end(); ++it) {
	UInt nb_nodes = _msh_nodes_per_elem[it->second];

	std::vector<UInt> tmp(nb_nodes);
	for (UInt i = 0; i < nb_nodes; ++i) {
	tmp[i] = i;
	}

	switch (it->first) {
	case _tetrahedron_10:
	tmp[8] = 9;
	tmp[9] = 8;
	break;
	case _pentahedron_6:
	tmp[0] = 2;
	tmp[1] = 0;
	tmp[2] = 1;
	tmp[3] = 5;
	tmp[4] = 3;
	tmp[5] = 4;
	break;
	case _pentahedron_15:
	tmp[0] = 2;
	tmp[1] = 0;
	tmp[2] = 1;
	tmp[3] = 5;
	tmp[4] = 3;
	tmp[5] = 4;
	tmp[6] = 8;
	tmp[8] = 11;
	tmp[9] = 6;
	tmp[10] = 9;
	tmp[11] = 10;
	tmp[12] = 14;
	tmp[14] = 12;
	break;
	case _hexahedron_20:
	tmp[9] = 11;
	tmp[10] = 12;
	tmp[11] = 9;
	tmp[12] = 13;
	tmp[13] = 10;
	tmp[17] = 19;
	tmp[18] = 17;
	tmp[19] = 18;
	break;
	default:
	// nothing to change
	break;
	}
	_read_order[it->first] = tmp;
	}
	}

	/* -------------------------------------------------------------------------- */
	MeshIOMSH::~MeshIOMSH() = default;

	/* -------------------------------------------------------------------------- */
	namespace {
	struct File {
	std::string filename;
	std::ifstream infile;
	std::string line;
	size_t current_line{0};

	size_t first_node_number{std::numeric_limits<UInt>::max()},
	last_node_number{0};
	bool has_physical_names{false};

	std::unordered_map<size_t, size_t> node_tags;
	std::unordered_map<size_t, Element> element_tags;
	double version;

	Mesh & mesh;
	MeshAccessor mesh_accessor;

	std::multimap<int, int> entity_tag_to_physical_tags;

	File(const std::string & filename, Mesh & mesh)
	: filename(filename), mesh(mesh), mesh_accessor(mesh) {
	infile.open(filename.c_str());
	if (not infile.good()) {
	AKANTU_EXCEPTION("Cannot open file " << filename);
	}
	}

	~File() { infile.close(); }

	auto good() { return infile.good(); }

	std::stringstream get_line() {
	std::string tmp_str;
	if (infile.eof()) {
	AKANTU_EXCEPTION("Reached the end of the file " << filename);
	}
	std::getline(infile, tmp_str);
	line = trim(tmp_str);
	++current_line;

	return std::stringstream(line);
	}

	template <typename... Ts> void read_line(Ts &&... ts) {
	auto && sstr = get_line();
	(void)std::initializer_list<int>{
	(sstr >> std::forward<decltype(ts)>(ts), 0)...};
	}
	};
	} // namespace

	/* -------------------------------------------------------------------------- */
	template <typename File, typename Readers>
	void MeshIOMSH::populateReaders2(File & file, Readers & readers) {
	readers["$NOD"] = readers["$Nodes"] = [&](const std::string &) {
	UInt nb_nodes;
	file.read_line(nb_nodes);

	Array<Real> & nodes = file.mesh_accessor.getNodes();
	nodes.resize(nb_nodes);
	file.mesh_accessor.setNbGlobalNodes(nb_nodes);

	size_t index;
	double coord[3];
	UInt spatial_dimension = nodes.getNbComponent();
	/// for each node, read the coordinates
	for (UInt i = 0; i < nb_nodes; ++i) {
	UInt offset = i * spatial_dimension;

	file.read_line(index, coord[0], coord[1], coord[2]);

	file.first_node_number = std::min(file.first_node_number, index);
	file.last_node_number = std::max(file.last_node_number, index);

	/// read the coordinates
	for (UInt j = 0; j < spatial_dimension; ++j)
	nodes.storage()[offset + j] = coord[j];
	}
	};

	readers["$ELM"] = readers["$Elements"] = [&](const std::string &) {
	UInt nb_elements;
	file.read_line(nb_elements);

	Int index;
	UInt msh_type;
	ElementType akantu_type;

	for (UInt i = 0; i < nb_elements; ++i) {
	auto && sstr_elem = file.get_line();

	sstr_elem >> index;
	sstr_elem >> msh_type;

	/// get the connectivity vector depending on the element type
	akantu_type =
	this->_msh_to_akantu_element_types[MSHElementType(msh_type)];

	if (akantu_type == _not_defined) {
	AKANTU_DEBUG_WARNING("Unsuported element kind "
	<< msh_type << " at line " << file.current_line);
	continue;
	}

	auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type);
	auto node_per_element = connectivity.getNbComponent();
	auto & read_order = this->_read_order[akantu_type];

	/// read tags informations
	if (file.version == 2) {
	UInt nb_tags;
	sstr_elem >> nb_tags;
	for (UInt j = 0; j < nb_tags; ++j) {
	Int tag;
	sstr_elem >> tag;

	auto & data = file.mesh_accessor.template getData<UInt>(
	"tag_" + std::to_string(j), akantu_type);
	data.push_back(tag);
	}
	} else if (file.version == 1) {
	Int tag0, tag1, nb_nodes; // reg-phys, reg-elem, number-of-nodes
	sstr_elem >> tag0 >> tag1 >> nb_nodes;

	auto & data0 =
	file.mesh_accessor.template getData<UInt>("tag_0", akantu_type);
	data0.push_back(tag0);

	auto & data1 =
	file.mesh_accessor.template getData<UInt>("tag_1", akantu_type);
	data1.push_back(tag1);
	}

	Vector<UInt> local_connect(node_per_element);
	for (UInt j = 0; j < node_per_element; ++j) {
	UInt node_index;
	sstr_elem >> node_index;

	AKANTU_DEBUG_ASSERT(node_index <= file.last_node_number,
	"Node number not in range : line "
	<< file.current_line);

	node_index -= file.first_node_number;
	local_connect(read_order[j]) = node_index;
	}
	connectivity.push_back(local_connect);
	}
	};

	readers["$NodeData"] = [&](const std::string &) {
	/* $NodeData
	number-of-string-tags
	< "string-tag" >
	…
	number-of-real-tags
	< real-tag >
	…
	number-of-integer-tags
	< integer-tag >
	…
	node-number value …
	…
	$EndNodeData */

	auto read_data_tags = [&](auto x) {
	UInt number_of_tags{0};
	file.read_line(number_of_tags);
	std::vector<decltype(x)> tags(number_of_tags);

	for (auto && tag : tags) {
	file.read_line(tag);
	}
	return tags;
	};

	auto && string_tags = read_data_tags(std::string{});
	auto && real_tags [[gnu::unused]] = read_data_tags(double{});
	auto && int_tags = read_data_tags(int{});

	auto && data = file.mesh.template registerNodalData<double>(
	trim(string_tags[0], '"'), int_tags[1]);
	data.resize(file.mesh.getNbNodes(), 0.);
	for (auto n [[gnu::unused]] : arange(int_tags[2])) {
	auto && sstr = file.get_line();

	int node;
	double value;
	sstr >> node;
	for (auto c : arange(int_tags[1])) {
	sstr >> value;
	data(node - file.first_node_number, c) = value;
	}
	}
	};

	readers["$Periodic"] = [&](const std::string &) {
	UInt nb_periodic_entities;
	file.read_line(nb_periodic_entities);

	file.mesh_accessor.getNodesFlags().resize(file.mesh.getNbNodes(),
	NodeFlag::_normal);

	for (UInt p = 0; p < nb_periodic_entities; ++p) {
	// dimension slave-tag master-tag
	UInt dimension;
	file.read_line(dimension);

	// transformation
	file.get_line();

	// nb nodes
	UInt nb_nodes;
	file.read_line(nb_nodes);

	for (UInt n = 0; n < nb_nodes; ++n) {
	// slave master
	auto && sstr = file.get_line();

	// The info in the mesh seem inconsistent so they are ignored for now.
	continue;

	if (dimension == file.mesh.getSpatialDimension() - 1) {
	UInt slave, master;

	sstr >> slave;
	sstr >> master;
	file.mesh_accessor.addPeriodicSlave(slave, master);
	}
	}
	}

	// mesh_accessor.markMeshPeriodic();
	};
	}

	/* -------------------------------------------------------------------------- */
	template <typename File, typename Readers>
	void MeshIOMSH::populateReaders4(File & file, Readers & readers) {
	static std::map<int, std::string> entity_type{
	{0, "points"},
	{1, "curve"},
	{2, "surface"},
	{3, "volume"},
	};

	readers["$Entities"] = [&](const std::string &) {
	size_t num_entity[4];
	file.read_line(num_entity[0], num_entity[1], num_entity[2], num_entity[3]);

	for (auto entity_dim : arange(4)) {
	for (auto _ [[gnu::unused]] : arange(num_entity[entity_dim])) {
	auto && sstr = file.get_line();

	int tag;
	double min_x, min_y, min_z, max_x, max_y, max_z;
	size_t num_physical_tags;
	sstr >> tag >> min_x >> min_y >> min_z;

	if (entity_dim > 0 or file.version < 4.1) {
	sstr >> max_x >> max_y >> max_z;
	}

	sstr >> num_physical_tags;

	for (auto _ [[gnu::unused]] : arange(num_physical_tags)) {
	int phys_tag;
	sstr >> phys_tag;

	std::string physical_name;
	if (this->physical_names.find(phys_tag) ==
	this->physical_names.end()) {
	physical_name = "msh_block_" + std::to_string(phys_tag);
	} else {
	physical_name = this->physical_names[phys_tag];
	}

	if (not file.mesh.elementGroupExists(physical_name)) {
	file.mesh.createElementGroup(physical_name, entity_dim);
	} else {
	file.mesh.getElementGroup(physical_name).addDimension(entity_dim);
	}
	file.entity_tag_to_physical_tags.insert(
	std::make_pair(tag, phys_tag));
	}
	}
	}
	};

	readers["$Nodes"] = [&](const std::string &) {
	size_t num_blocks, num_nodes;
	if (file.version >= 4.1) {
	file.read_line(num_blocks, num_nodes, file.first_node_number,
	file.last_node_number);
	} else {
	file.read_line(num_blocks, num_nodes);
	}
	auto & nodes = file.mesh_accessor.getNodes();
	nodes.reserve(num_nodes);
	file.mesh_accessor.setNbGlobalNodes(num_nodes);

	if (num_nodes > std::numeric_limits<UInt>::max()) {
	AKANTU_EXCEPTION(
	"There are more nodes in this files than the index type of akantu "
	"can handle, consider recompiling with a bigger index type");
	}

	for (auto block [[gnu::unused]] : arange(num_blocks)) {
	int entity_dim, entity_tag, parametric;
	size_t num_nodes_in_block;
	size_t node_id{0};

	// auto & grp = file.mesh.createNodeGroup("msh_" + entity_type[entity_dim] +
	// "_block_" + std::to_string(block));

	if (file.version >= 4.1) {
	file.read_line(entity_dim, entity_tag, parametric, num_nodes_in_block);
	if (parametric) {
	AKANTU_EXCEPTION(
	"Akantu does not support parametric nodes in msh files");
	}
	for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
	size_t tag;
	file.read_line(tag);
	file.node_tags[tag] = node_id;
	//grp.add(node_id);
	++node_id;
	}

	for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
	Vector<double> pos(3);
	file.read_line(pos(_x), pos(_y), pos(_z));
	nodes.push_back(pos);
	}
	} else {
	file.read_line(entity_tag, entity_dim, parametric, num_nodes_in_block);

	for (auto _ [[gnu::unused]] : arange(num_nodes_in_block)) {
	size_t tag;
	Vector<double> pos(3);
	file.read_line(tag, pos(_x), pos(_y), pos(_z));

	if (file.version < 4.1) {
	file.first_node_number = std::min(file.first_node_number, tag);
	file.last_node_number = std::max(file.last_node_number, tag);
	}

	nodes.push_back(pos);
	file.node_tags[tag] = node_id;
	//grp.add(node_id);
	++node_id;
	}
	}
	}
	};

	readers["$Elements"] = [&](const std::string &) {
	size_t num_blocks, num_elements;
	file.read_line(num_blocks, num_elements);

	for (auto block [[gnu::unused]] : arange(num_blocks)) {
	int entity_dim, entity_tag, element_type;
	size_t num_elements_in_block;

	if (file.version >= 4.1) {
	file.read_line(entity_dim, entity_tag, element_type,
	num_elements_in_block);
	} else {
	file.read_line(entity_tag, entity_dim, element_type,
	num_elements_in_block);
	}

	/// get the connectivity vector depending on the element type
	auto && akantu_type =
	this->_msh_to_akantu_element_types[(MSHElementType)element_type];

	if (akantu_type == _not_defined) {
	AKANTU_DEBUG_WARNING("Unsuported element kind " << element_type
	<< " at line "
	<< file.current_line);
	continue;
	}

	Element elem{akantu_type, 0, _not_ghost};

	auto & connectivity = file.mesh_accessor.getConnectivity(akantu_type);
	Vector<UInt> local_connect(connectivity.getNbComponent());
	auto && read_order = this->_read_order[akantu_type];

	for (auto _ [[gnu::unused]] : arange(num_elements_in_block)) {
	auto && sstr_elem = file.get_line();
	size_t elem_tag;
	sstr_elem >> elem_tag;
	for (auto && c : arange(connectivity.getNbComponent())) {
	size_t node_tag;
	sstr_elem >> node_tag;

	AKANTU_DEBUG_ASSERT(node_tag <= file.last_node_number,
	"Node number not in range : line "
	<< file.current_line);

	node_tag = file.node_tags[node_tag];
	local_connect(read_order[c]) = node_tag;
	}
	connectivity.push_back(local_connect);
	elem.element = connectivity.size() - 1;
	file.element_tags[elem_tag] = elem;

	auto range = file.entity_tag_to_physical_tags.equal_range(entity_tag);
	for (auto it = range.first; it != range.second; ++it) {
	file.mesh.getElementGroup(this->physical_names[it->second]).add(elem);
	}
	}
	}
	};
	}

	/* -------------------------------------------------------------------------- */
	void MeshIOMSH::read(const std::string & filename, Mesh & mesh) {

	File file(filename, mesh);

	std::map<std::string, std::function<void(const std::string &)>> readers;

	readers["$MeshFormat"] = [&](const std::string &) {
	auto && sstr = file.get_line();

	int format;
	sstr >> file.version >> format;

	if (format != 0)
	AKANTU_ERROR("This reader can only read ASCII files.");

	if (file.version > 2) {
	int size_of_size_t;
	sstr >> size_of_size_t;
	if (size_of_size_t > int(sizeof(UInt))) {
	AKANTU_DEBUG_WARNING("The size of the indexes in akantu are to small "
	"to read this file (akantu "
	<< sizeof(UInt) << " vs. msh file "
	<< size_of_size_t << ")");
	}
	}

	if (file.version < 4) {
	this->populateReaders2(file, readers);
	} else {
	this->populateReaders4(file, readers);
	}
	};

	readers["$PhysicalNames"] = [&](const std::string &) {
	file.has_physical_names = true;
	int num_of_phys_names;
	file.read_line(num_of_phys_names); /// the format line

	for (auto k [[gnu::unused]] : arange(num_of_phys_names)) {
	int phys_name_id;
	int phys_dim;
	std::string phys_name;
	file.read_line(phys_dim, phys_name_id, std::quoted(phys_name));

	this->physical_names[phys_name_id] = phys_name;
	}
	};

	readers["Unsupported"] = [&](const std::string & _block) {
	std::string block = _block.substr(1);
	AKANTU_DEBUG_WARNING("Unsupported block_kind " << block << " at line "
	<< file.current_line);
	auto && end_block = "$End" + block;
	while (file.line != end_block) {
	file.get_line();
	}
	};

	while (file.good()) {
	std::string block;
	file.read_line(block);

	auto && it = readers.find(block);

	if (it != readers.end()) {
	it->second(block);

	std::string end_block;
	file.read_line(end_block);
	block = block.substr(1);

	if (end_block != "$End" + block) {
	AKANTU_EXCEPTION("The reader failed to properly read the block "
	<< block << ". Expected a $End" << block << " at line "
	<< file.current_line);
	}
	} else if (block.size() != 0) {
	readers["Unsupported"](block);
	}
	}

	// mesh.updateTypesOffsets(_not_ghost);

	if (file.version <= 2) {
	this->constructPhysicalNames("tag_0", mesh);
	if (file.has_physical_names)
	mesh.createGroupsFromMeshData<std::string>("physical_names");
	}

	MeshUtils::fillElementToSubElementsData(mesh);
	}

	/* -------------------------------------------------------------------------- */
	void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) {
	std::ofstream outfile;
	const Array<Real> & nodes = mesh.getNodes();

	outfile.open(filename.c_str());

	outfile << "$MeshFormat"
	<< "\n";
	outfile << "2.1 0 8"
	<< "\n";
	outfile << "$EndMeshFormat"
	<< "\n";

	outfile << std::setprecision(std::numeric_limits<Real>::digits10);
	outfile << "$Nodes"
	<< "\n";

	outfile << nodes.size() << "\n";

	outfile << std::uppercase;
	for (UInt i = 0; i < nodes.size(); ++i) {
	Int offset = i * nodes.getNbComponent();
	outfile << i + 1;
	for (UInt j = 0; j < nodes.getNbComponent(); ++j) {
	outfile << " " << nodes.storage()[offset + j];
	}

	for (UInt p = nodes.getNbComponent(); p < 3; ++p)
	outfile << " " << 0.;
	outfile << "\n";
	;
	}
	outfile << std::nouppercase;
	outfile << "$EndNodes"
	<< "\n";

	outfile << "$Elements"
	<< "\n";
	Int nb_elements = 0;
	for (auto && type :
	mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) {
	const Array<UInt> & connectivity = mesh.getConnectivity(type, _not_ghost);
	nb_elements += connectivity.size();
	}
	outfile << nb_elements << "\n";

	UInt element_idx = 1;
	for (auto && type :
	mesh.elementTypes(_all_dimensions, _not_ghost, _ek_not_defined)) {
	const auto & connectivity = mesh.getConnectivity(type, _not_ghost);

	UInt * tag[2] = {nullptr, nullptr};
	if (mesh.hasData<UInt>("tag_0", type, _not_ghost)) {
	const auto & data_tag_0 = mesh.getData<UInt>("tag_0", type, _not_ghost);
	tag[0] = data_tag_0.storage();
	}

	if (mesh.hasData<UInt>("tag_1", type, _not_ghost)) {
	const auto & data_tag_1 = mesh.getData<UInt>("tag_1", type, _not_ghost);
	tag[1] = data_tag_1.storage();
	}

	for (UInt i = 0; i < connectivity.size(); ++i) {
	UInt offset = i * connectivity.getNbComponent();
	outfile << element_idx << " " << _akantu_to_msh_element_types[type]
	<< " 2";

	/// \todo write the real data in the file
	for (UInt t = 0; t < 2; ++t)
	if (tag[t])
	outfile << " " << tag[t][i];
	else
	outfile << " 0";

	for (UInt j = 0; j < connectivity.getNbComponent(); ++j) {
	outfile << " " << connectivity.storage()[offset + j] + 1;
	}
	outfile << "\n";
	element_idx++;
	}
	}

	outfile << "$EndElements"
	<< "\n";

	if (mesh.hasData(MeshDataType::_nodal)) {
	auto && tags = mesh.getTagNames();
	for (auto && tag : tags) {
	if (mesh.getTypeCode(tag, MeshDataType::_nodal) != _tc_real)
	continue;

	auto && data = mesh.getNodalData<double>(tag);
	outfile << "$NodeData"
	<< "\n";
	outfile << "1"
	<< "\n";
	outfile << "\"" << tag << "\"\n";
	outfile << "1\n0.0"
	<< "\n";
	outfile << "3\n0"
	<< "\n";
	outfile << data.getNbComponent() << "\n";
	outfile << data.size() << "\n";
	for (auto && d : enumerate(make_view(data, data.getNbComponent()))) {
	outfile << std::get<0>(d) + 1;
	for (auto && v : std::get<1>(d)) {
	outfile << " " << v;
	}
	outfile << "\n";
	}
	outfile << "$EndNodeData"
	<< "\n";
	}
	}

	outfile.close();
	}

	/* --------------------------------------------------------------------------
	*/

	} // namespace akantu

mesh_io_msh.ccNo OneTemporaryActions

File Metadata

mesh_io_msh.ccView Options

Event Timeline

mesh_io_msh.cc
No OneTemporary
Actions

mesh_io_msh.cc
View Options