Homec4science
Diffusion akantu (bugfixes/46-bug-in-strain-energy-split-and-missing-constraint-in-phasefield) Browse

rAKA/third-party/iohelper/test/mesh_io_msh.cc
09b079271f45bugfixes/46-bug-in…

Last Change
Raw File

mesh_io_msh.cc
View Options

/**
* @file mesh_io_msh.cc
*
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Thu Oct 11 2012
* @date last modification: Thu Nov 01 2012
*
* @brief Read/Write for MSH files generated by gmsh
*
*
* Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* IOHelper 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.
*
* IOHelper 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 IOHelper. 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 <limits>
#include <iostream>
#include <cstdlib>
#include <sstream>
#include <cassert>
/* -------------------------------------------------------------------------- */
#include "mesh_io_msh.hh"
/* -------------------------------------------------------------------------- */
typedef unsigned int UInt;
typedef int Int;
typedef double Real;
/* -------------------------------------------------------------------------- */
#define MESH_IO_FATAL(x) { \
std::cerr << "Cannot open file " << filename << std::endl; \
exit(EXIT_FAILURE); \
}
/* -------------------------------------------------------------------------- */
/* Methods Implentations */
/* -------------------------------------------------------------------------- */
MeshIOMSH::MeshIOMSH() {
_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 ] = 1;
_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_prism_18 ] = 18;
_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 ] = _not_defined;
// _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_pyramid_14 ] = _not_defined;
// _msh_to_akantu_element_types[_msh_point ] = _not_defined;
// _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[_point ] = _msh_point;
// _akantu_to_msh_element_types[_bernoulli_beam_2] = _msh_segment_2;
std::map<MSHElementType, UInt>::iterator it;
for(it = _msh_nodes_per_elem.begin();
it != _msh_nodes_per_elem.end(); ++it) {
UInt nb_nodes = (*it).second;
UInt * tmp = new UInt[nb_nodes];
for (UInt i = 0; i < nb_nodes; ++i) {
tmp[i] = i;
}
switch((*it).first) {
case _msh_tetrahedron_10:
tmp[8] = 9;
tmp[9] = 8;
break;
default:
//nothing to change
break;
}
_read_order[it->first] = tmp;
}
}
/* -------------------------------------------------------------------------- */
MeshIOMSH::~MeshIOMSH() {
std::map<MSHElementType, UInt>::iterator it;
for(it = _msh_nodes_per_elem.begin();
it != _msh_nodes_per_elem.end(); ++it) {
delete [] _read_order[it->first];
}
}
/* -------------------------------------------------------------------------- */
void MeshIOMSH::read(const std::string & filename,
UInt spatial_dimension,
MSHElementType selected_type,
std::vector<double> & nodes,
std::vector<int> & connectivity){
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()) {
MESH_IO_FATAL("Cannot open file ");
}
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) MESH_IO_FATAL("This reader can only read ASCII files.");
std::getline(infile, line); /// the end of block line
current_line += 2;
file_format = 2;
}
/// read all nodes
if(line == "$Nodes" || line == "$NOD") {
UInt nb_nodes;
std::getline(infile, line);
std::stringstream sstr(line);
sstr >> nb_nodes;
current_line++;
nodes.resize(nb_nodes*spatial_dimension);
UInt index;
Real coord[3];
/// 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[offset + j] = coord[j];
}
std::getline(infile, line); /// the end of block line
}
/// read all elements
if(line == "$Elements" || line == "$ELM") {
std::getline(infile, line);
std::stringstream sstr(line);
UInt nb_elements = 0;
sstr >> nb_elements;
current_line++;
Int index;
UInt msh_type;
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;
if (selected_type != msh_type) continue;
UInt * read_order = _read_order[static_cast<MSHElementType>(msh_type)];
UInt node_per_element = _msh_nodes_per_elem[static_cast<MSHElementType>(msh_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;
}
} else if (file_format == 1) {
Int tag;
sstr_elem >> tag; //reg-phys
std::string tag_name = "tag_0";
sstr_elem >> tag; //reg-elem
tag_name = "tag_1";
}
UInt local_connect[node_per_element];
for(UInt j = 0; j < node_per_element; ++j) {
UInt node_index;
sstr_elem >> node_index;
assert(node_index <= last_node_number);
node_index -= first_node_number;
local_connect[read_order[j]] = node_index;
}
for (UInt i = 0; i < node_per_element; ++i) {
connectivity.push_back(local_connect[i]);
}
}
std::getline(infile, line); /// the end of block line
}
if((line[0] == '$') && (line.find("End") == std::string::npos)) {
std::cerr << "Unsuported block_kind " << line
<< " at line " << current_line << std::endl;
}
}
infile.close();
}
/* -------------------------------------------------------------------------- */
c4science · Help