Forked from
ogs / ogs
22349 commits behind the upstream repository.
-
Marc Walther authored
now reads physical surface id, before read plane surface id
Marc Walther authorednow reads physical surface id, before read plane surface id
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
GMSHInterface.cpp 14.15 KiB
/**
* \file
* \author Thomas Fischer
* \date 2010-04-29
* \brief Implementation of the GMSHInterface class.
*
* \copyright
* Copyright (c) 2012-2015, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
* @file GMSHInterface.cpp
* @date 2010-04-29
* @author Thomas Fischer
*/
#include <fstream>
#include <vector>
// ThirdParty/logog
#include "logog/include/logog.hpp"
// BaseLib
#include "BaseLib/BuildInfo.h"
#include "FileTools.h"
#include "StringTools.h"
// FileIO
#include "GMSHInterface.h"
#include "GmshIO/GMSHAdaptiveMeshDensity.h"
#include "GmshIO/GMSHFixedMeshDensity.h"
#include "GmshIO/GMSHNoMeshDensity.h"
// GeoLib
#include "Point.h"
#include "Polygon.h"
#include "Polyline.h"
#include "PolylineWithSegmentMarker.h"
#include "QuadTree.h"
// MSH
#include "Elements/Line.h"
#include "Elements/Hex.h"
#include "Elements/Prism.h"
#include "Elements/Pyramid.h"
#include "Elements/Quad.h"
#include "Elements/Tet.h"
#include "Elements/Tri.h"
#include "Mesh.h"
#include "MeshLib/Node.h"
namespace FileIO
{
GMSHInterface::GMSHInterface(GeoLib::GEOObjects & geo_objs,
bool /*include_stations_as_constraints*/,
GMSH::MeshDensityAlgorithm mesh_density_algorithm,
double param1,
double param2,
std::size_t param3,
std::vector<std::string>& selected_geometries) :
_n_lines(0), _n_plane_sfc(0), _geo_objs(geo_objs), _selected_geometries(selected_geometries)
{
switch (mesh_density_algorithm) {
case GMSH::MeshDensityAlgorithm::NoMeshDensity:
_mesh_density_strategy = new GMSH::GMSHNoMeshDensity;
break;
case GMSH::MeshDensityAlgorithm::FixedMeshDensity:
_mesh_density_strategy = new GMSH::GMSHFixedMeshDensity(param1);
break; case GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity:
_mesh_density_strategy = new GMSH::GMSHAdaptiveMeshDensity(param1, param2, param3);
break;
}
}
int GMSHInterface::writeGeoFile(GeoLib::GEOObjects &geo_objects, std::string const& file_name)
{
std::vector<std::string> names;
geo_objects.getGeometryNames(names);
if (names.empty())
{
ERR ("No geometry information available.");
return 1;
}
bool const multiple_geometries = (names.size() > 1);
std::string merge_name("MergedGeometry");
if (multiple_geometries)
{
if (geo_objects.mergeGeometries (names, merge_name) != 1)
return 2;
names.clear();
names.push_back(merge_name);
}
else
merge_name = names[0];
// default parameters for GMSH interface
double param1(0.5); // mesh density scaling on normal points
double param2(0.05); // mesh density scaling on station points
std::size_t param3(2); // points per leaf
GMSHInterface gmsh_io(geo_objects, true, FileIO::GMSH::MeshDensityAlgorithm::AdaptiveMeshDensity, param1, param2, param3, names);
int const writer_return_val = gmsh_io.writeToFile(file_name);
if (multiple_geometries)
{
geo_objects.removeSurfaceVec(merge_name);
geo_objects.removePolylineVec(merge_name);
geo_objects.removePointVec(merge_name);
}
return (writer_return_val==1) ? 0 : 3;
}
bool GMSHInterface::isGMSHMeshFile(const std::string& fname)
{
std::ifstream input(fname.c_str());
if (!input) {
ERR("GMSHInterface::isGMSHMeshFile(): Could not open file %s.", fname.c_str());
return false;
}
std::string header_first_line;
input >> header_first_line;
if (header_first_line.find("$MeshFormat") != std::string::npos) {
// read version
std::string version;
getline(input, version);
getline(input, version);
INFO("GMSHInterface::isGMSHMeshFile(): Found GMSH mesh file version: %s.",
version.c_str());
input.close();
return true;
}
return false;
}
MeshLib::Mesh* GMSHInterface::readGMSHMesh(std::string const& fname)
{
std::string line;
std::ifstream in(fname.c_str(), std::ios::in);
if (!in.is_open())
{
WARN ("GMSHInterface::readGMSHMesh() - Could not open file %s.", fname.c_str());
return nullptr;
}
getline(in, line); // $MeshFormat keyword
if (line.find("$MeshFormat") == std::string::npos)
{
in.close();
WARN ("No GMSH file format recognized.");
return nullptr;
}
getline(in, line); // version-number file-type data-size
if (line.substr(0,3).compare("2.2") != 0) {
WARN("Wrong gmsh file format version.");
return nullptr;
}
if (line.substr(4,1).compare("0") != 0) {
WARN("Currently reading gmsh binary file type is not supported.");
return nullptr;
}
getline(in, line); //$EndMeshFormat
std::vector<MeshLib::Node*> nodes;
std::vector<MeshLib::Element*> elements;
std::vector<int> materials;
std::map<unsigned, unsigned> id_map;
while (line.find("$EndElements") == std::string::npos)
{
// Node data
getline(in, line); //$Nodes Keywords
if (line.find("$Nodes") != std::string::npos)
{
std::size_t n_nodes(0);
long id;
double x, y, z;
in >> n_nodes >> std::ws;
nodes.resize(n_nodes);
for (std::size_t i = 0; i < n_nodes; i++) {
in >> id >> x >> y >> z >> std::ws;
id_map.insert(std::map<unsigned, unsigned>::value_type(id, i));
nodes[i] = new MeshLib::Node(x,y,z,id);
}
getline(in, line); // End Node keyword $EndNodes
}
// Element data
if (line.find("$Elements") != std::string::npos)
{
std::size_t n_elements(0);
if (! (in >> n_elements >> std::ws)) { // number-of-elements
ERR("Read GMSH mesh does not contain any elements");
}
elements.reserve(n_elements);
materials.reserve(n_elements);
for (std::size_t i = 0; i < n_elements; i++)
{
MeshLib::Element* elem(nullptr);
int mat_id(0);
std::tie(elem, mat_id) = readElement(in, nodes, id_map);
if (elem) {
elements.push_back(elem); materials.push_back(mat_id);
}
}
getline(in, line); // END keyword
}
if (line.find("PhysicalNames") != std::string::npos)
{
std::size_t n_lines(0);
in >> n_lines >> std::ws; // number-of-lines
for (std::size_t i = 0; i < n_lines; i++)
getline(in, line);
getline(in, line); // END keyword
}
}
in.close();
if (elements.empty()) {
for (auto it(nodes.begin()); it != nodes.end(); ++it) {
delete *it;
}
return nullptr;
}
MeshLib::Mesh * mesh(new MeshLib::Mesh(
BaseLib::extractBaseNameWithoutExtension(fname), nodes, elements));
boost::optional<MeshLib::PropertyVector<int> &> opt_material_ids(
mesh->getProperties().createNewPropertyVector<int>(
"MaterialIDs", MeshLib::MeshItemType::Cell, 1)
);
if (!opt_material_ids) {
WARN("Could not create PropertyVector for MaterialIDs in Mesh.");
} else {
MeshLib::PropertyVector<int> & material_ids(opt_material_ids.get());
material_ids.insert(material_ids.end(), materials.cbegin(),
materials.cend());
}
INFO("\t... finished.");
INFO("Nr. Nodes: %d.", nodes.size());
INFO("Nr. Elements: %d.", elements.size());
return mesh;
}
void GMSHInterface::readNodeIDs(std::ifstream &in,
unsigned n_nodes,
std::vector<unsigned> &node_ids,
std::map<unsigned, unsigned> const& id_map)
{
unsigned idx;
for (unsigned i = 0; i < n_nodes; i++)
{
in >> idx;
node_ids.push_back(id_map.at(idx));
}
}
std::pair<MeshLib::Element*, int>
GMSHInterface::readElement(std::ifstream &in,
std::vector<MeshLib::Node*> const& nodes,
std::map<unsigned, unsigned> const& id_map)
{
unsigned idx, type, n_tags, dummy;
int mat_id;
std::vector<unsigned> node_ids;
std::vector<MeshLib::Node*> elem_nodes;
in >> idx >> type >> n_tags >> mat_id >> dummy;
// skip tags
for (std::size_t j = 2; j < n_tags; j++) in >> dummy;
switch (type)
{
case 1: {
readNodeIDs(in, 2, node_ids, id_map);
// edge_nodes array will be deleted from Line object
MeshLib::Node** edge_nodes = new MeshLib::Node*[2];
edge_nodes[0] = nodes[node_ids[0]];
edge_nodes[1] = nodes[node_ids[1]];
return std::make_pair(new MeshLib::Line(edge_nodes), 0);
}
case 2: {
readNodeIDs(in, 3, node_ids, id_map);
MeshLib::Node** tri_nodes = new MeshLib::Node*[3];
tri_nodes[0] = nodes[node_ids[2]];
tri_nodes[1] = nodes[node_ids[1]];
tri_nodes[2] = nodes[node_ids[0]];
return std::make_pair(new MeshLib::Tri(tri_nodes), mat_id);
}
case 3: {
readNodeIDs(in, 4, node_ids, id_map);
MeshLib::Node** quad_nodes = new MeshLib::Node*[4];
for (unsigned k(0); k < 4; k++)
quad_nodes[k] = nodes[node_ids[k]];
return std::make_pair(new MeshLib::Quad(quad_nodes), mat_id);
}
case 4: {
readNodeIDs(in, 4, node_ids, id_map);
MeshLib::Node** tet_nodes = new MeshLib::Node*[5];
for (unsigned k(0); k < 4; k++)
tet_nodes[k] = nodes[node_ids[k]];
return std::make_pair(new MeshLib::Tet(tet_nodes), mat_id);
}
case 5: {
readNodeIDs(in, 8, node_ids, id_map);
MeshLib::Node** hex_nodes = new MeshLib::Node*[8];
for (unsigned k(0); k < 8; k++)
hex_nodes[k] = nodes[node_ids[k]];
return std::make_pair(new MeshLib::Hex(hex_nodes), mat_id);
}
case 6: {
readNodeIDs(in, 6, node_ids, id_map);
MeshLib::Node** prism_nodes = new MeshLib::Node*[6];
for (unsigned k(0); k < 6; k++)
prism_nodes[k] = nodes[node_ids[k]];
return std::make_pair(new MeshLib::Prism(prism_nodes), mat_id);
}
case 7: {
readNodeIDs(in, 5, node_ids, id_map);
MeshLib::Node** pyramid_nodes = new MeshLib::Node*[5];
for (unsigned k(0); k < 5; k++)
pyramid_nodes[k] = nodes[node_ids[k]];
return std::make_pair(new MeshLib::Pyramid(pyramid_nodes), mat_id);
}
case 15:
in >> dummy; // skip rest of line
break;
default:
WARN("GMSHInterface::readGMSHMesh(): Unknown element type %d.", type);
}
return std::make_pair(nullptr, -1);
}
bool GMSHInterface::write()
{
_out << "// GMSH input file created by OpenGeoSys " << BaseLib::BuildInfo::git_describe;
#ifdef BUILD_TIMESTAMP
_out << " built on " << BaseLib::BuildInfo::build_timestamp;#endif
_out << "\n\n";
writeGMSHInputFile(_out);
return true;
}
void GMSHInterface::writeGMSHInputFile(std::ostream& out)
{
DBUG("GMSHInterface::writeGMSHInputFile(): get data from GEOObjects.");
if (_selected_geometries.empty())
return;
bool remove_geometry(false);
// *** get and merge data from _geo_objs
if (_selected_geometries.size() > 1) {
_gmsh_geo_name = "GMSHGeometry";
remove_geometry = true;
_geo_objs.mergeGeometries(_selected_geometries, _gmsh_geo_name);
} else {
_gmsh_geo_name = _selected_geometries[0];
remove_geometry = false;
}
std::vector<GeoLib::Point*> * merged_pnts(const_cast<std::vector<GeoLib::Point*> *>(_geo_objs.getPointVec(_gmsh_geo_name)));
if (! merged_pnts) {
ERR("GMSHInterface::writeGMSHInputFile(): Did not found any points.");
return;
} else {
const std::size_t n_pnts(merged_pnts->size());
for (std::size_t k(0); k<n_pnts; k++) {
(*((*merged_pnts)[k]))[2] = 0.0;
}
}
std::vector<GeoLib::Polyline*> const* merged_plys(_geo_objs.getPolylineVec(_gmsh_geo_name));
DBUG("GMSHInterface::writeGMSHInputFile(): \t ok.");
// *** compute topological hierarchy of polygons
if (merged_plys) {
for (std::vector<GeoLib::Polyline*>::const_iterator it(merged_plys->begin());
it!=merged_plys->end(); ++it) {
if ((*it)->isClosed()) {
_polygon_tree_list.push_back(new GMSH::GMSHPolygonTree(new GeoLib::Polygon(*(*it), true), NULL, _geo_objs, _gmsh_geo_name, _mesh_density_strategy));
}
}
DBUG("GMSHInterface::writeGMSHInputFile(): Compute topological hierarchy - detected %d polygons.", _polygon_tree_list.size());
GeoLib::createPolygonTrees<GMSH::GMSHPolygonTree>(_polygon_tree_list);
DBUG("GMSHInterface::writeGMSHInputFile(): Compute topological hierarchy - calculated %d polygon trees.", _polygon_tree_list.size());
} else {
return;
}
// *** insert stations and polylines (except polygons) in the appropriate object of
// class GMSHPolygonTree
// *** insert stations
const std::size_t n_geo_names(_selected_geometries.size());
for (std::size_t j(0); j < n_geo_names; j++) {
const std::vector<GeoLib::Point*>* stations (_geo_objs.getStationVec(_selected_geometries[j]));
if (stations) {
const std::size_t n_stations(stations->size());
for (std::size_t k(0); k < n_stations; k++) {
bool found(false);
for (std::list<GMSH::GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
it != _polygon_tree_list.end() && !found; ++it) {
if ((*it)->insertStation((*stations)[k])) {
found = true;
}
}
}
} }
// *** insert polylines
const std::size_t n_plys(merged_plys->size());
for (std::size_t k(0); k<n_plys; k++) {
if (! (*merged_plys)[k]->isClosed()) {
for (std::list<GMSH::GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
it != _polygon_tree_list.end(); ++it) {
(*it)->insertPolyline(new GeoLib::PolylineWithSegmentMarker(*(*merged_plys)[k]));
}
}
}
// *** init mesh density strategies
for (std::list<GMSH::GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
it != _polygon_tree_list.end(); ++it) {
(*it)->initMeshDensityStrategy();
}
// *** create GMSH data structures
const std::size_t n_merged_pnts(merged_pnts->size());
_gmsh_pnts.resize(n_merged_pnts);
for (std::size_t k(0); k<n_merged_pnts; k++) {
_gmsh_pnts[k] = NULL;
}
for (std::list<GMSH::GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
it != _polygon_tree_list.end(); ++it) {
(*it)->createGMSHPoints(_gmsh_pnts);
}
// *** finally write data :-)
writePoints(out);
std::size_t pnt_id_offset(_gmsh_pnts.size());
for (std::list<GMSH::GMSHPolygonTree*>::iterator it(_polygon_tree_list.begin());
it != _polygon_tree_list.end(); ++it) {
(*it)->writeLineLoop(_n_lines, _n_plane_sfc, out);
(*it)->writeSubPolygonsAsLineConstraints(_n_lines, _n_plane_sfc-1, out);
(*it)->writeLineConstraints(_n_lines, _n_plane_sfc-1, out);
(*it)->writeStations(pnt_id_offset, _n_plane_sfc-1, out);
(*it)->writeAdditionalPointData(pnt_id_offset, _n_plane_sfc-1, out);
}
if (remove_geometry) {
_geo_objs.removeSurfaceVec(_gmsh_geo_name);
_geo_objs.removePolylineVec(_gmsh_geo_name);
_geo_objs.removePointVec(_gmsh_geo_name);
}
}
void GMSHInterface::writePoints(std::ostream& out) const
{
const std::size_t n_gmsh_pnts(_gmsh_pnts.size());
for (std::size_t k(0); k<n_gmsh_pnts; k++) {
if (_gmsh_pnts[k]) {
out << *(_gmsh_pnts[k]) << "\n";
}
}
}
} // end namespace FileIO