swh:1:snp:6088ab52ef49920e01e3f334cdf4d5d6c8a822b9
Tip revision: 718b674cc0af413d7a07a561c7533c9988097d93 authored by Lars Bilke on 23 March 2021, 10:20:21 UTC
[cmake] Simplified (MSVC) folder setup.
[cmake] Simplified (MSVC) folder setup.
Tip revision: 718b674
TetGenInterface.cpp
/**
* \file
* \author Thomas Fischer
* \date 2011-09-12
* \brief Implementation of the TetGenInterface class.
*
* \copyright
* Copyright (c) 2012-2021, OpenGeoSys Community (http://www.opengeosys.org)
* Distributed under a Modified BSD License.
* See accompanying file LICENSE.txt or
* http://www.opengeosys.org/project/license
*
*/
#include "TetGenInterface.h"
#include <cstddef>
#include <string>
#include <fstream>
#include "BaseLib/Logging.h"
#include "BaseLib/FileTools.h"
#include "BaseLib/StringTools.h"
#include "GeoLib/Triangle.h"
#include "MeshLib/Mesh.h"
#include "MeshLib/Node.h"
#include "MeshLib/Elements/Element.h"
#include "MeshLib/Elements/Tet.h"
#include "MeshLib/MeshInformation.h"
namespace FileIO
{
TetGenInterface::TetGenInterface() = default;
bool TetGenInterface::readTetGenGeometry (std::string const& geo_fname,
GeoLib::GEOObjects &geo_objects)
{
std::ifstream poly_stream (geo_fname.c_str());
if (!poly_stream)
{
ERR("TetGenInterface::readTetGenGeometry() failed to open {:s}",
geo_fname);
return false;
}
std::string ext (BaseLib::getFileExtension(geo_fname));
if (ext != ".smesh")
{
ERR ("TetGenInterface::readTetGenGeometry() - unknown file type (only *.smesh is supported).");
return false;
}
std::vector<MeshLib::Node*> nodes;
if (!readNodesFromStream (poly_stream, nodes))
{
// remove nodes read until now
for (auto& node : nodes)
{
delete node;
}
return false;
}
const std::size_t nNodes (nodes.size());
auto points = std::make_unique<std::vector<GeoLib::Point*>>();
points->reserve(nNodes);
for (std::size_t k(0); k<nNodes; ++k)
{
points->push_back(new GeoLib::Point(*(nodes[k]), nodes[k]->getID()));
delete nodes[k];
}
std::string geo_name (BaseLib::extractBaseNameWithoutExtension(geo_fname));
geo_objects.addPointVec(std::move(points), geo_name);
const std::vector<std::size_t> &id_map (geo_objects.getPointVecObj(geo_name)->getIDMap());
auto surfaces = std::make_unique<std::vector<GeoLib::Surface*>>();
if (!parseSmeshFacets(poly_stream, *surfaces, *geo_objects.getPointVec(geo_name), id_map))
{
// remove surfaces read until now but keep the points
for (std::size_t k = 0; k < surfaces->size(); k++)
{
delete (*surfaces)[k];
}
}
geo_objects.addSurfaceVec(std::move(surfaces), geo_name);
return true;
}
std::size_t TetGenInterface::getNFacets(std::ifstream &input)
{
std::string line;
while (!input.fail())
{
getline (input, line);
if (input.fail())
{
ERR("TetGenInterface::getNFacets(): Error reading number of facets.");
return 0;
}
BaseLib::simplify(line);
if (line.empty() || line.compare(0, 1, "#") == 0)
{
continue;
}
const std::list<std::string> fields = BaseLib::splitString(line, ' ');
auto it = fields.begin();
const auto nFacets(BaseLib::str2number<std::size_t>(*it));
if (fields.size() > 1)
{
_boundary_markers = BaseLib::str2number<std::size_t>(*(++it)) != 0;
}
return nFacets;
}
return 0;
}
bool TetGenInterface::parseSmeshFacets(std::ifstream &input,
std::vector<GeoLib::Surface*> &surfaces,
const std::vector<GeoLib::Point*> &points,
const std::vector<std::size_t> &pnt_id_map)
{
const std::size_t nFacets (this->getNFacets(input));
std::string line;
surfaces.reserve(nFacets);
std::list<std::string>::const_iterator it;
const unsigned offset = (_zero_based_idx) ? 0 : 1;
std::vector<std::size_t> idx_map;
std::size_t k(0);
while (k<nFacets && !input.fail())
{
getline (input, line);
if (input.fail())
{
ERR("TetGenInterface::parseFacets(): Error reading facet {:d}.", k);
return false;
}
BaseLib::simplify(line);
if (line.empty() || line.compare(0,1,"#") == 0)
{
continue;
}
// read facets
const std::list<std::string> point_fields = BaseLib::splitString(line, ' ');
it = point_fields.begin();
const auto nPoints = BaseLib::str2number<std::size_t>(*it);
if (nPoints != 3)
{
ERR ("Smesh-files are currently only supported for triangle meshes.");
return false;
}
std::vector<std::size_t> point_ids;
const std::size_t point_field_size = (_boundary_markers) ? nPoints+1 : nPoints;
if (point_fields.size() > point_field_size)
{
for (std::size_t j(0); j < nPoints; ++j)
{
point_ids.push_back(
pnt_id_map[BaseLib::str2number<std::size_t>(*(++it)) -
offset]);
}
const std::size_t sfc_marker = (_boundary_markers) ? BaseLib::str2number<std::size_t>(*(++it)) : 0;
const std::size_t idx = std::find(idx_map.begin(), idx_map.end(), sfc_marker) - idx_map.begin();
if (idx >= surfaces.size())
{
idx_map.push_back(sfc_marker);
surfaces.push_back(new GeoLib::Surface(points));
}
surfaces[idx]->addTriangle(point_ids[0], point_ids[1], point_ids[2]);
}
else
{
ERR("TetGenInterface::parseFacets(): Error reading points for "
"facet {:d}.",
k);
return false;
}
++k;
}
// here the poly-file potentially defines a number of points to mark holes within the volumes defined by the facets, these are ignored for now
// here the poly-file potentially defines a number of region attributes, these are ignored for now
std::size_t nTotalTriangles (0);
for (auto& surface : surfaces)
{
nTotalTriangles += surface->getNumberOfTriangles();
}
if (nTotalTriangles == nFacets)
{
return true;
}
ERR("TetGenInterface::parseFacets(): Number of expected total triangles "
"({:d}) does not match number of found triangles ({:d}).",
surfaces.size(),
nTotalTriangles);
return false;
}
MeshLib::Mesh* TetGenInterface::readTetGenMesh (std::string const& nodes_fname,
std::string const& ele_fname)
{
std::ifstream ins_nodes (nodes_fname.c_str());
std::ifstream ins_ele (ele_fname.c_str());
if (!ins_nodes || !ins_ele)
{
if (!ins_nodes)
{
ERR("TetGenInterface::readTetGenMesh failed to open {:s}",
nodes_fname);
}
if (!ins_ele)
{
ERR("TetGenInterface::readTetGenMesh failed to open {:s}",
ele_fname);
}
return nullptr;
}
std::vector<MeshLib::Node*> nodes;
if (!readNodesFromStream (ins_nodes, nodes)) {
// remove nodes read until now
for (auto & node : nodes) {
delete node;
}
return nullptr;
}
std::vector<MeshLib::Element*> elements;
std::vector<int> materials;
if (!readElementsFromStream (ins_ele, elements, materials, nodes)) {
BaseLib::cleanupVectorElements(nodes, elements);
return nullptr;
}
MeshLib::Properties properties;
// Transmit material values if there is any material value != 0
if (std::any_of(materials.cbegin(), materials.cend(), [](int m){ return m != 0; }))
{
auto* const mat_props = properties.createNewPropertyVector<int>(
"MaterialIDs", MeshLib::MeshItemType::Cell);
mat_props->reserve(elements.size());
std::copy(materials.cbegin(),
materials.cend(),
std::back_inserter(*mat_props));
}
const std::string mesh_name (BaseLib::extractBaseNameWithoutExtension(nodes_fname));
return new MeshLib::Mesh(mesh_name, nodes, elements, properties);
}
bool TetGenInterface::readNodesFromStream (std::ifstream &ins,
std::vector<MeshLib::Node*> &nodes)
{
std::string line;
getline (ins, line);
std::size_t n_nodes;
std::size_t dim;
std::size_t n_attributes;
bool boundary_markers;
while (!ins.fail())
{
BaseLib::simplify(line);
if (line.empty() || line.compare(0,1,"#") == 0)
{
// this line is a comment - skip
getline (ins, line);
continue;
}
// read header line
bool header_okay = parseNodesFileHeader(line, n_nodes, dim, n_attributes, boundary_markers);
if (!header_okay)
{
return false;
}
if (!parseNodes(ins, nodes, n_nodes, dim))
{
return false;
}
return true;
}
return false;
}
bool TetGenInterface::parseNodesFileHeader(std::string const& line,
std::size_t& n_nodes,
std::size_t& dim,
std::size_t& n_attributes,
bool& boundary_markers)
{
std::list<std::string> pnt_header = BaseLib::splitString(line, ' ');
if (pnt_header.empty())
{
ERR("TetGenInterface::parseNodesFileHeader(): could not read number of nodes specified in header.");
return false;
}
auto it = pnt_header.begin();
n_nodes = BaseLib::str2number<std::size_t> (*it);
dim = (pnt_header.size()==1) ? 3 : BaseLib::str2number<std::size_t> (*(++it));
if (pnt_header.size()<4)
{
n_attributes = 0;
boundary_markers = false;
return true;
}
n_attributes = BaseLib::str2number<std::size_t> (*(++it));
boundary_markers = *(++it) == "1";
return true;
}
bool TetGenInterface::parseNodes(std::ifstream &ins,
std::vector<MeshLib::Node*> &nodes,
std::size_t n_nodes,
std::size_t dim)
{
std::string line;
nodes.reserve(n_nodes);
std::size_t k(0);
while (k < n_nodes && !ins.fail())
{
std::vector<double> coordinates(dim);
std::getline(ins, line);
if (ins.fail())
{
ERR("TetGenInterface::parseNodes(): Error reading node {:d}.", k);
return false;
}
std::size_t id;
std::size_t pos_end = 0;
std::size_t pos_beg = line.find_first_not_of(' ', pos_end);
pos_end = line.find_first_of(" \n", pos_beg);
if (line.empty() || pos_beg==pos_end || line.compare(pos_beg,1,"#") == 0)
{
continue;
}
if (pos_beg != std::string::npos && pos_end != std::string::npos) {
id = BaseLib::str2number<std::size_t> (line.substr(pos_beg, pos_end - pos_beg));
if (k == 0 && id == 0)
{
_zero_based_idx = true;
}
} else {
ERR("TetGenInterface::parseNodes(): Error reading ID of node {:d}.",
k);
return false;
}
// read coordinates
const unsigned offset = (_zero_based_idx) ? 0 : 1;
for (std::size_t i(0); i < dim; i++) {
pos_beg = line.find_first_not_of(' ', pos_end);
pos_end = line.find_first_of(" \n", pos_beg);
if (pos_end == std::string::npos)
{
pos_end = line.size();
}
if (pos_beg != std::string::npos)
{
coordinates[i] = BaseLib::str2number<double>(line.substr(pos_beg, pos_end-pos_beg));
}
else
{
ERR("TetGenInterface::parseNodes(): error reading coordinate "
"{:d} of node {:d}.",
i,
k);
return false;
}
}
nodes.push_back(new MeshLib::Node(coordinates.data(), id-offset));
// read attributes and boundary markers ... - at the moment we do not use this information
++k;
}
return true;
}
bool TetGenInterface::readElementsFromStream(
std::ifstream& ins,
std::vector<MeshLib::Element*>& elements,
std::vector<int>& materials,
const std::vector<MeshLib::Node*>& nodes) const
{
std::string line;
std::getline(ins, line);
std::size_t n_tets;
std::size_t n_nodes_per_tet;
bool region_attributes;
while (!ins.fail())
{
BaseLib::simplify(line);
if (line.empty() || line.compare(0, 1, "#") == 0)
{
// this line is a comment - skip
std::getline(ins, line);
continue;
}
// read header line
bool header_okay = parseElementsFileHeader(
line, n_tets, n_nodes_per_tet, region_attributes);
if (!header_okay)
{
return false;
}
if (!parseElements(ins,
elements,
materials,
nodes,
n_tets,
n_nodes_per_tet,
region_attributes))
{
return false;
}
return true;
}
return false;
}
bool TetGenInterface::parseElementsFileHeader(std::string& line,
std::size_t& n_tets,
std::size_t& n_nodes_per_tet,
bool& region_attribute)
{
std::size_t pos_beg;
std::size_t pos_end;
// number of tetrahedras
pos_beg = line.find_first_not_of (' ');
pos_end = line.find_first_of(' ', pos_beg);
if (pos_beg != std::string::npos && pos_end != std::string::npos)
{
n_tets = BaseLib::str2number<std::size_t> (line.substr(pos_beg, pos_end - pos_beg));
}
else
{
ERR("TetGenInterface::parseElementsFileHeader(): Could not read number of tetrahedra specified in header.");
return false;
}
// nodes per tet - either 4 or 10
pos_beg = line.find_first_not_of (" \t", pos_end);
pos_end = line.find_first_of(" \t", pos_beg);
n_nodes_per_tet = BaseLib::str2number<std::size_t> (line.substr(pos_beg, pos_end - pos_beg));
// region attribute at tetrahedra?
pos_beg = line.find_first_not_of (" \t", pos_end);
pos_end = line.find_first_of(" \t\n", pos_beg);
if (pos_end == std::string::npos)
{
pos_end = line.size();
}
region_attribute = line.substr(pos_beg, pos_end - pos_beg) == "1";
return true;
}
bool TetGenInterface::parseElements(std::ifstream& ins,
std::vector<MeshLib::Element*> &elements,
std::vector<int> &materials,
const std::vector<MeshLib::Node*> &nodes,
std::size_t n_tets,
std::size_t n_nodes_per_tet,
bool region_attribute) const
{
std::string line;
std::vector<std::size_t> ids(n_nodes_per_tet);
elements.reserve(n_tets);
materials.reserve(n_tets);
const unsigned offset = (_zero_based_idx) ? 0 : 1;
for (std::size_t k(0); k < n_tets && !ins.fail(); k++)
{
getline (ins, line);
if (ins.fail())
{
ERR("TetGenInterface::parseElements(): Error reading tetrahedron "
"{:d}.",
k);
return false;
}
std::size_t pos_end = 0;
std::size_t pos_beg = line.find_first_not_of(' ', pos_end);
pos_end = line.find_first_of(" \n", pos_beg);
if (line.empty() || pos_beg==pos_end || line.compare(pos_beg,1,"#") == 0)
{
k--;
continue;
}
if (pos_beg == std::string::npos || pos_end == std::string::npos)
{
ERR("TetGenInterface::parseElements(): Error reading id of "
"tetrahedron {:d}.",
k);
return false;
}
// read node ids
for (std::size_t i(0); i < n_nodes_per_tet; i++)
{
pos_beg = line.find_first_not_of(' ', pos_end);
pos_end = line.find_first_of(' ', pos_beg);
if (pos_end == std::string::npos)
{
pos_end = line.size();
}
if (pos_beg != std::string::npos && pos_end != std::string::npos)
{
ids[i] = BaseLib::str2number<std::size_t>(line.substr(pos_beg, pos_end - pos_beg)) - offset;
}
else
{
ERR("TetGenInterface::parseElements(): Error reading node {:d} "
"of tetrahedron {:d}.",
i,
k);
return false;
}
}
// read region attribute - this is something like material group
int region (0);
if (region_attribute) {
pos_beg = line.find_first_not_of(' ', pos_end);
pos_end = line.find_first_of(' ', pos_beg);
if (pos_end == std::string::npos)
{
pos_end = line.size();
}
if (pos_beg != std::string::npos && pos_end != std::string::npos)
{
region = BaseLib::str2number<int> (line.substr(pos_beg, pos_end - pos_beg));
}
else
{
ERR("TetGenInterface::parseElements(): Error reading region "
"attribute of tetrahedron {:d}.",
k);
return false;
}
}
// insert new element into vector
auto** tet_nodes = new MeshLib::Node*[4];
for (int n = 0; n < 4; n++)
{
tet_nodes[n] = nodes[ids[n]];
}
elements.push_back (new MeshLib::Tet(tet_nodes));
materials.push_back(region);
}
return true;
}
bool TetGenInterface::writeTetGenSmesh(
const std::string& file_name,
const GeoLib::GEOObjects& geo_objects,
const std::string& geo_name,
const std::vector<GeoLib::Point>& attribute_points)
{
std::vector<GeoLib::Point*> const*const points = geo_objects.getPointVec(geo_name);
std::vector<GeoLib::Surface*> const*const surfaces = geo_objects.getSurfaceVec(geo_name);
if (points==nullptr)
{
ERR("Geometry {:s} not found.", geo_name);
return false;
}
if (surfaces==nullptr)
{
WARN("No surfaces found for geometry {:s}. Writing points only.",
geo_name);
}
std::ofstream out( file_name.c_str(), std::ios::out );
out.precision(std::numeric_limits<double>::digits10);
// the points header
const std::size_t nPoints (points->size());
out << nPoints << " 3\n";
// the point list
for (std::size_t i = 0; i < nPoints; ++i)
{
out << i << " " << (*(*points)[i])[0] << " " << (*(*points)[i])[1]
<< " " << (*(*points)[i])[2] << "\n";
}
// the surfaces header
const std::size_t nSurfaces = (surfaces) ? surfaces->size() : 0;
std::size_t nTotalTriangles (0);
for (std::size_t i = 0; i < nSurfaces; ++i)
{
nTotalTriangles += (*surfaces)[i]->getNumberOfTriangles();
}
out << nTotalTriangles << " 1\n";
for (std::size_t i=0; i<nSurfaces; ++i)
{
const std::size_t nTriangles ((*surfaces)[i]->getNumberOfTriangles());
const std::size_t marker (i+1); // must NOT be 0!
// the poly list
for (std::size_t j=0; j<nTriangles; ++j)
{
const GeoLib::Triangle &tri = *(*(*surfaces)[i])[j];
out << "3 " << tri[0] << " " << tri[1] << " " << tri[2] << " " << marker << "\n";
}
}
out << "0\n"; // the polygon holes list
// the region attributes list
if (attribute_points.empty())
{
out << "0\n";
}
else
{
const std::size_t nAttributePoints (attribute_points.size());
out << nAttributePoints << "\n";
for (std::size_t i = 0; i < nAttributePoints; ++i)
{
out << i + 1 << " " << attribute_points[i][0] << " "
<< attribute_points[i][1] << " " << attribute_points[i][2]
<< " " << 10 * attribute_points[i].getID() << "\n";
}
}
INFO(
"TetGenInterface::writeTetGenSmesh() - {:d} points and {:d} surfaces "
"successfully written.",
nPoints,
nSurfaces);
out.close();
return true;
}
bool TetGenInterface::writeTetGenSmesh(const std::string &file_name,
const MeshLib::Mesh &mesh,
std::vector<MeshLib::Node> &attribute_points) const
{
if (mesh.getDimension() == 1)
{
return false;
}
const std::vector<MeshLib::Node*> &nodes = mesh.getNodes();
std::ofstream out( file_name.c_str(), std::ios::out );
out.precision(std::numeric_limits<double>::digits10);
// the points header
const std::size_t nPoints (nodes.size());
out << nPoints << " 3\n";
// the point list
for (std::size_t i = 0; i < nPoints; ++i)
{
out << i << " " << (*nodes[i])[0] << " " << (*nodes[i])[1] << " "
<< (*nodes[i])[2] << "\n";
}
if (mesh.getDimension() == 2)
{
write2dElements(out, mesh);
}
else
{
write3dElements(out, mesh, attribute_points);
}
out << "0\n"; // the polygon holes list
// the region attributes list
if (attribute_points.empty())
{
out << "0\n";
}
else
{
const std::size_t nAttributePoints (attribute_points.size());
out << nAttributePoints << "\n";
for (std::size_t i = 0; i < nAttributePoints; ++i)
{
out << i + 1 << " " << attribute_points[i][0] << " "
<< attribute_points[i][1] << " " << attribute_points[i][2]
<< " " << 10 * attribute_points[i].getID() << "\n";
}
}
INFO(
"TetGenInterface::writeTetGenPoly() - {:d} points and {:d} surfaces "
"successfully written.",
nPoints,
mesh.getNumberOfElements());
out.close();
return true;
}
void TetGenInterface::write2dElements(std::ofstream &out,
const MeshLib::Mesh &mesh) const
{
// the surfaces header
auto const& types = MeshLib::MeshInformation::getNumberOfElementTypes(mesh);
std::size_t const n_tri =
(types.find(MeshLib::MeshElemType::TRIANGLE) != types.end())
? types.at(MeshLib::MeshElemType::TRIANGLE) : 0;
std::size_t const n_quad =
(types.find(MeshLib::MeshElemType::QUAD) != types.end())
? types.at(MeshLib::MeshElemType::QUAD) : 0;
unsigned const nTotalTriangles = n_tri + 2 * n_quad;
out << nTotalTriangles << " 1\n";
const std::vector<MeshLib::Element*> &elements = mesh.getElements();
MeshLib::PropertyVector<int> const*const mat_ids =
mesh.getProperties().existsPropertyVector<int>("MaterialIDs")
? mesh.getProperties().getPropertyVector<int>("MaterialIDs")
: nullptr;
const std::size_t nElements (elements.size());
unsigned element_count(0);
for (std::size_t i=0; i<nElements; ++i)
{
std::string const matId = mat_ids ? std::to_string((*mat_ids)[i]) : "";
this->writeElementToFacets(out, *elements[i], element_count, matId);
}
}
void TetGenInterface::write3dElements(std::ofstream &out,
const MeshLib::Mesh &mesh,
std::vector<MeshLib::Node> &attribute_points) const
{
const std::vector<MeshLib::Element*> &elements = mesh.getElements();
const std::size_t nElements (elements.size());
if (!attribute_points.empty())
{
attribute_points.clear();
}
// get position where number of facets need to be written and figure out worst case of chars that are needed
const std::streamoff before_elems_pos (out.tellp());
const unsigned n_spaces (static_cast<unsigned>(std::floor(log(nElements*8))) + 1);
out << std::string(n_spaces, ' ') << " 1\n";
auto const* const materialIds =
mesh.getProperties().getPropertyVector<int>("MaterialIDs");
unsigned element_count(0);
for (std::size_t i=0; i<nElements; ++i)
{
if (elements[i]->getDimension() < 3)
{
continue;
}
const unsigned nFaces (elements[i]->getNumberOfNeighbors());
std::string const mat_id_str =
materialIds ? std::to_string((*materialIds)[i]) : "";
for (std::size_t j=0; j<nFaces; ++j)
{
MeshLib::Element const*const neighbor ( elements[i]->getNeighbor(j) );
if (neighbor && materialIds &&
(*materialIds)[i] <= (*materialIds)[neighbor->getID()])
{
continue;
}
std::unique_ptr<MeshLib::Element const> const face (elements[i]->getFace(j));
this->writeElementToFacets(out, *face, element_count, mat_id_str);
}
if (materialIds)
{
attribute_points.emplace_back(
MeshLib::getCenterOfGravity(*elements[i]).getCoords(),
(*materialIds)[i]);
}
}
// add number of facets at correct position and jump back
const std::streamoff after_elems_pos (out.tellp());
out.seekp(before_elems_pos);
out << element_count;
out.seekp(after_elems_pos);
}
void TetGenInterface::writeElementToFacets(std::ofstream& out,
const MeshLib::Element& element,
unsigned& element_count,
std::string const& matId)
{
element_count++;
if (element.getGeomType() == MeshLib::MeshElemType::TRIANGLE)
{
out << "3 " << element.getNodeIndex(0) << " " << element.getNodeIndex(1) << " " << element.getNodeIndex(2) << " " << matId << " # " << element_count << "\n";
}
else if (element.getGeomType() == MeshLib::MeshElemType::QUAD)
{
out << "3 " << element.getNodeIndex(0) << " " << element.getNodeIndex(1) << " " << element.getNodeIndex(2) << " " << matId << " # " << element_count << "\n";
element_count++;
out << "3 " << element.getNodeIndex(0) << " " << element.getNodeIndex(2) << " " << element.getNodeIndex(3) << " " << matId << " # " << element_count << "\n";
}
}
} // end namespace FileIO