Homec4science
Diffusion akantu (bugfixes/39-the-python-example-for-phasefield-are-broken-due-to-a-change-to-setintegrationscheme) Browse

rAKA/src/io/mesh_io/mesh_io_msh.cc
1b154aa804debugfixes/39-the-python…

Last Change
Raw File

mesh_io_msh.cc
View Options

/**
* @file mesh_io_msh.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Fri Jun 18 2010
* @date last modification: Fri Jul 04 2014
*
* @brief Read/Write for MSH files generated by gmsh
*
* @section LICENSE
*
* Copyright (©) 2014 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 <fstream>
/* -------------------------------------------------------------------------- */
#include "mesh_io.hh"
#include "mesh_utils.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_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_prism_18 ] = _not_defined;
_msh_to_akantu_element_types[_msh_hexahedron_20 ] = _hexahedron_20;
_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[_kirchhoff_shell] = _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];
UInt * tmp = new UInt[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_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() {
std::map<ElementType, MSHElementType>::iterator it;
for(it = _akantu_to_msh_element_types.begin();
it != _akantu_to_msh_element_types.end(); ++it) {
delete [] _read_order[it->first];
}
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::read(const std::string & filename, Mesh & mesh) {
std::ifstream infile;
infile.open(filename.c_str());
std::string line;
UInt first_node_number = std::numeric_limits<UInt>::max(), last_node_number = 0,
file_format = 1, current_line = 0;
if(!infile.good()) {
AKANTU_DEBUG_ERROR("Cannot open file " << filename);
}
while(infile.good()) {
std::getline(infile, line);
current_line++;
/// read the header
if(line == "$MeshFormat") {
std::getline(infile, line); /// the format line
std::stringstream sstr(line);
std::string version; sstr >> version;
Int format; sstr >> format;
if(format != 0) AKANTU_DEBUG_ERROR("This reader can only read ASCII files.");
std::getline(infile, line); /// the end of block line
current_line += 2;
file_format = 2;
}
/// read the physical names
if(line == "$PhysicalNames") {
std::getline(infile, line); /// the format line
std::stringstream sstr(line);
UInt num_of_phys_names;
sstr >> num_of_phys_names;
for(UInt k(0); k < num_of_phys_names; k++) {
std::getline(infile, line);
std::stringstream sstr_phys_name(line);
UInt phys_name_id;
UInt phys_dim;
sstr_phys_name >> phys_dim >> phys_name_id;
std::size_t b = line.find("\"");
std::size_t e = line.rfind("\"");
std::string phys_name = line.substr(b + 1, e - b -1);
phys_name_map[phys_name_id] = phys_name;
}
}
/// read all nodes
if(line == "$Nodes" || line == "$NOD") {
UInt nb_nodes;
std::getline(infile, line);
std::stringstream sstr(line);
sstr >> nb_nodes;
current_line++;
Array<Real> & nodes = const_cast<Array<Real> &>(mesh.getNodes());
nodes.resize(nb_nodes);
mesh.nb_global_nodes = nb_nodes;
UInt index;
Real 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;
std::getline(infile, line);
std::stringstream sstr_node(line);
sstr_node >> index >> coord[0] >> coord[1] >> coord[2];
current_line++;
first_node_number = std::min(first_node_number,index);
last_node_number = std::max(last_node_number, index);
/// read the coordinates
for(UInt j = 0; j < spatial_dimension; ++j)
nodes.storage()[offset + j] = coord[j];
}
std::getline(infile, line); /// the end of block line
}
/// read all elements
if(line == "$Elements" || line == "$ELM") {
UInt nb_elements;
UInt * read_order = NULL;
std::getline(infile, line);
std::stringstream sstr(line);
sstr >> nb_elements;
current_line++;
Int index;
UInt msh_type;
ElementType akantu_type, akantu_type_old = _not_defined;
Array<UInt> *connectivity = NULL;
UInt node_per_element = 0;
for(UInt i = 0; i < nb_elements; ++i) {
std::getline(infile, line);
std::stringstream sstr_elem(line);
current_line++;
sstr_elem >> index;
sstr_elem >> msh_type;
/// get the connectivity vector depending on the element type
akantu_type = _msh_to_akantu_element_types[(MSHElementType) msh_type];
if(akantu_type == _not_defined) {
AKANTU_DEBUG_WARNING("Unsuported element kind " << msh_type
<< " at line " << current_line);
continue;
}
if(akantu_type != akantu_type_old) {
connectivity = mesh.getConnectivityPointer(akantu_type);
// connectivity->resize(0);
node_per_element = connectivity->getNbComponent();
akantu_type_old = akantu_type;
read_order = _read_order[akantu_type];
}
/// read tags informations
if(file_format == 2) {
UInt nb_tags;
sstr_elem >> nb_tags;
for(UInt j = 0; j < nb_tags; ++j) {
Int tag;
sstr_elem >> tag;
std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << j;
Array<UInt> * data = mesh.getDataPointer<UInt>(sstr_tag_name.str(), akantu_type, _not_ghost);
data->push_back(tag);
}
} else if (file_format == 1) {
Int tag;
sstr_elem >> tag; //reg-phys
std::string tag_name = "tag_0";
Array<UInt> * data = mesh.getDataPointer<UInt>(tag_name, akantu_type, _not_ghost);
data->push_back(tag);
sstr_elem >> tag; //reg-elem
tag_name = "tag_1";
data = mesh.getDataPointer<UInt>(tag_name, akantu_type, _not_ghost);
data->push_back(tag);
sstr_elem >> tag; //number-of-nodes
}
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 <= last_node_number,
"Node number not in range : line " << current_line);
node_index -= first_node_number;
local_connect[read_order[j]] = node_index;
}
connectivity->push_back(local_connect);
}
std::getline(infile, line); /// the end of block line
}
if((line[0] == '$') && (line.find("End") == std::string::npos)) {
AKANTU_DEBUG_WARNING("Unsuported block_kind " << line
<< " at line " << current_line);
}
}
mesh.updateTypesOffsets(_not_ghost);
infile.close();
this->constructPhysicalNames("tag_0",mesh);
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" << std::endl;
outfile << "2.1 0 8" << std::endl;;
outfile << "$EndMeshFormat" << std::endl;;
outfile << std::setprecision(std::numeric_limits<Real>::digits10);
outfile << "$Nodes" << std::endl;;
outfile << nodes.getSize() << std::endl;
outfile << std::uppercase;
for(UInt i = 0; i < nodes.getSize(); ++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 << std::endl;;
}
outfile << std::nouppercase;
outfile << "$EndNodes" << std::endl;;
outfile << "$Elements" << std::endl;;
Mesh::type_iterator it = mesh.firstType(_all_dimensions, _not_ghost, _ek_not_defined);
Mesh::type_iterator end = mesh.lastType(_all_dimensions, _not_ghost, _ek_not_defined);
Int nb_elements = 0;
for(; it != end; ++it) {
const Array<UInt> & connectivity = mesh.getConnectivity(*it, _not_ghost);
nb_elements += connectivity.getSize();
}
outfile << nb_elements << std::endl;
UInt element_idx = 1;
for(it = mesh.firstType(_all_dimensions, _not_ghost, _ek_not_defined); it != end; ++it) {
ElementType type = *it;
const Array<UInt> & connectivity = mesh.getConnectivity(type, _not_ghost);
UInt * tag[2] = {NULL, NULL};
try {
const Array<UInt> & data_tag_0 = mesh.getData<UInt>("tag_0", type, _not_ghost);
tag[0] = data_tag_0.storage();
} catch(...) { tag[0] = NULL; }
try {
const Array<UInt> & data_tag_1 = mesh.getData<UInt>("tag_1", type, _not_ghost);
tag[1] = data_tag_1.storage();
} catch(...) { tag[1] = NULL; }
for(UInt i = 0; i < connectivity.getSize(); ++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 << std::endl;
element_idx++;
}
}
outfile << "$EndElements" << std::endl;;
outfile.close();
}
/* -------------------------------------------------------------------------- */
__END_AKANTU__
c4science · Help