Commit 46d41374 authored by Stenger, Florian's avatar Stenger, Florian
Browse files

new branch for physicalTags and msh2-format

parent 396a57a0
Pipeline #6174 failed with stage
in 10 minutes and 7 seconds
......@@ -8,4 +8,4 @@ Version: 0.2
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-grid
Suggests: dune-alugrid dune-foamgrid dune-vtk dune-curvedsurfacegrid
Suggests: dune-alugrid dune-foamgrid dune-vtk dune-curvedgrid
#pragma once
#include <iosfwd>
#include <set>
#include <map>
#include <memory>
#include <vector>
......@@ -10,6 +11,8 @@
#include <dune/gmsh4/filereader.hh>
#include <dune/gmsh4/gridcreators/continuousgridcreator.hh> // default GridCreator
#include <dune/grid/io/file/gmshreader.hh> //for old msh2.2 reader
namespace Dune
{
/// File-Reader for GMsh unstructured .msh files
......@@ -42,6 +45,10 @@ namespace Dune
public:
using size_type = SizeType;
Gmsh4Reader ()
: creator_(new GridCreator)
{}
/// Constructor. Creates a new GridCreator with the passed factory
template <class... Args,
std::enable_if_t<std::is_constructible<GridCreator, Args...>::value,int> = 0>
......@@ -59,6 +66,39 @@ namespace Dune
: creator_(std::move(creator))
{}
//convenience-method to keep compatibility with old GmshReader from dune-grid.
static ToUniquePtr<Grid> read (const std::string& fileName, bool verbose,
bool insertBoundarySegments);
//convenience-method to keep compatibility with old GmshReader from dune-grid.
static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
bool verbose = true, bool insertBoundarySegments=true);
//Needed for following read-method. Also declared in GmshReader. Should probably by public
//there are accessible via friend-declaration...
struct DataArg {
std::vector<int> *data_ = nullptr;
DataArg(std::vector<int> &data) : data_(&data) {}
DataArg(const decltype(std::ignore)&) {}
DataArg() = default;
};
struct DataFlagArg : DataArg {
bool flag_ = false;
using DataArg::DataArg;
DataFlagArg(bool flag) : flag_(flag) {}
};
//convenience-method to keep compatibility with old GmshReader from dune-grid.
static void read (Dune::GridFactory<Grid> &factory, const std::string &fileName,
DataFlagArg boundarySegmentData, DataArg elementData, bool verbose=true);
//convenience-method to keep compatibility with old GmshReader from dune-grid.
static void read (Dune::GridFactory<Grid>& factory, const std::string& fileName,
std::vector<int>& boundarySegmentToPhysicalEntity,
std::vector<int>& elementToPhysicalEntity,
bool verbose=true, bool insertBoundarySegments=true);
/// Read the grid from file with `filename` into the GridFactory stored in the GridCreator
/**
* \param filename The name of the input file
......@@ -108,6 +148,8 @@ namespace Dune
#endif
private:
void prepareBoundingEntityTags();
std::vector<int>& getPhysicalTags(int entityDim, int entityTag);
template<class T>
void readValueBinary(std::ifstream& input, T &v);
void readMeshFormat (std::ifstream& input,
......@@ -228,11 +270,15 @@ namespace Dune
// PhysicalNames section
std::vector<PhysicalNamesAttributes> physicalNames_;
// boundary handling
std::set<int> boundingEntityTags;
std::set<int> boundaryEntities;
// Entities section
std::vector<PointAttributes> points_;
std::vector<EntityAttributes> curves_;
std::vector<EntityAttributes> surfaces_;
std::vector<EntityAttributes> volumes_;
std::map<int, PointAttributes> points_;
std::map<int, EntityAttributes> curves_;
std::map<int, EntityAttributes> surfaces_;
std::map<int, EntityAttributes> volumes_;
// PartitionedEntities section
size_type numPartitions_ = 0;
......@@ -267,6 +313,10 @@ namespace Dune
};
// deduction guides
template <class Grid>
Gmsh4Reader ()
-> Gmsh4Reader<Grid, Gmsh4::ContinuousGridCreator<Grid>>;
template <class Grid>
Gmsh4Reader (GridFactory<Grid>&)
-> Gmsh4Reader<Grid, Gmsh4::ContinuousGridCreator<Grid>>;
......
......@@ -9,6 +9,8 @@
#include <dune/gmsh4/utility/filesystem.hh>
#include <dune/gmsh4/utility/string.hh>
#include <dune/grid/io/file/gmshreader.hh> //for old msh2.2 reader
// Helper-function to deal with endianness in binary msh-files.
inline void swapBytes (char *array, int size)
{
......@@ -240,6 +242,7 @@ struct Gmsh4Reader<G,C,S>::ElementAttributes
int entityDim;
int entityTag;
int elementType;
std::vector<int> physicalTags;
std::vector<Element> elements;
};
......@@ -260,6 +263,38 @@ struct Gmsh4Reader<G,C,S>::PeriodicAttributes
};
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::prepareBoundingEntityTags()
{
switch((int)G::dimension){
case 1:
for(auto &curv : curves_)
for(auto tag : curv.second.boundingEntities) boundingEntityTags.insert(tag);
break;
case 2:
for(auto &surface : surfaces_)
for(auto tag : surface.second.boundingEntities) boundingEntityTags.insert(tag);
break;
case 3:
for(auto &volume : volumes_)
for(auto tag : volume.second.boundingEntities) boundingEntityTags.insert(tag);
}
}
template <class G, class C, class S>
std::vector<int>& Gmsh4Reader<G,C,S>::getPhysicalTags(int entityDim, int entityTag)
{
switch(entityDim){
case 0: return points_[entityTag].physicalTags;
case 1: return curves_[entityTag].physicalTags;
case 2: return surfaces_[entityTag].physicalTags;
case 3: return volumes_[entityTag].physicalTags;
}
DUNE_THROW(Dune::NotImplemented, "Only entities of dimension 0, 1, 2 and 3 allowed.");
}
template <class G, class C, class S>
template <class T>
void Gmsh4Reader<G,C,S>::readValueBinary(std::ifstream& input, T &v)
......@@ -270,6 +305,62 @@ void Gmsh4Reader<G,C,S>::readValueBinary(std::ifstream& input, T &v)
}
//convenience-method to keep compatibility with old GmshReader from dune-grid.
template <class G, class C, class S>
ToUniquePtr<G> Gmsh4Reader<G,C,S>::read (const std::string& fileName, bool verbose,
bool insertBoundarySegments)
{
std::cout << "WARNING! You are using a deprecated syntax for the GmshReader! "
"Use \"read(std::string const& filename, bool fillCreator = true)\" instead, "
"especially if you need to read MSH 4.1 files." << std::endl;
return GmshReader<G>::read(fileName, verbose, insertBoundarySegments);
}
//convenience-function to keep compatibility with old GmshReader from dune-grid.
template<class T> static T &discarded(T &&value) { return value; }
//convenience-method to keep compatibility with old GmshReader from dune-grid.
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::read (Dune::GridFactory<G>& factory, const std::string& fileName,
bool verbose, bool insertBoundarySegments)
{
std::cout << "WARNING! You are using a deprecated syntax for the GmshReader! "
"Use \"read(std::string const& filename, bool fillCreator = true)\" instead, "
"especially if you need to read MSH 4.1 files." << std::endl;
GmshReader<G>::read(factory, fileName, discarded(std::vector<int>{}),
discarded(std::vector<int>{}), verbose, insertBoundarySegments);
}
//convenience-method to keep compatibility with old GmshReader from dune-grid.
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::read (Dune::GridFactory<G> &factory, const std::string &fileName,
DataFlagArg boundarySegmentData, DataArg elementData, bool verbose)
{
std::cout << "WARNING! You are using a deprecated syntax for the GmshReader! "
"Use \"read(std::string const& filename, bool fillCreator = true)\" instead, "
"especially if you need to read MSH 4.1 files." << std::endl;
GmshReader<G>::read(factory, fileName, boundarySegmentData, elementData, verbose);
}
//convenience-method to keep compatibility with old GmshReader from dune-grid.
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::read (Dune::GridFactory<G>& factory, const std::string& fileName,
std::vector<int>& boundarySegmentToPhysicalEntity,
std::vector<int>& elementToPhysicalEntity,
bool verbose, bool insertBoundarySegments)
{
std::cout << "WARNING! You are using a deprecated syntax for the GmshReader! "
"Use \"read(std::string const& filename, bool fillCreator = true)\" instead, "
"especially if you need to read MSH 4.1 files." << std::endl;
GmshReader<G>::read(factory, fileName, boundarySegmentToPhysicalEntity,
elementToPhysicalEntity, verbose, insertBoundarySegments);
}
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::read (std::string const& filename, bool fillCreator)
{
......@@ -277,13 +368,15 @@ void Gmsh4Reader<G,C,S>::read (std::string const& filename, bool fillCreator)
if (!Gmsh4::exists(filename))
DUNE_THROW(IOError, "File " << filename << " does not exist!");
clear();
std::ifstream input(filename, std::ios_base::in | std::ios_base::binary);
GMSH4_ASSERT(input.is_open());
std::string ext = Gmsh4::Path(filename).extension().string();
if (ext == ".msh") {
readSerialFileFromStream(input, fillCreator);
pieces_.push_back(filename);
readSerialFileFromStream(input, fillCreator);
} else if (ext == ".pro") {
readParallelFileFromStream(input, comm().rank(), comm().size(), fillCreator);
} else {
......@@ -295,8 +388,6 @@ void Gmsh4Reader<G,C,S>::read (std::string const& filename, bool fillCreator)
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::readSerialFileFromStream (std::ifstream& input, bool fillCreator)
{
clear();
// MeshFormat section
double version = 0.0;
int file_type = 0;
......@@ -320,14 +411,23 @@ void Gmsh4Reader<G,C,S>::readSerialFileFromStream (std::ifstream& input, bool fi
switch (section) {
case Sections::MESH_FORMAT: {
readMeshFormat(input, version, file_type, data_size);
GMSH4_ASSERT_MSG(version >= 4.0 && version < 5.0,
"Can only read gmsh files versions >= 4.0 and < 5.0");
GMSH4_ASSERT_MSG((version >= 4.0 && version < 5.0) || (version >= 2.0 && version <= 2.2),
"Can only read gmsh files versions >= 4.0 and < 5.0 or versions >= 2.0 and <= 2.2");
if( version >= 2.0 && version <= 2.2 ){
std::vector<int> boundarySegmentToPhysicalEntity, elementToPhysicalEntity;
GmshReader<G>::read(creator_->factory(), pieces_.back(),
boundarySegmentToPhysicalEntity, elementToPhysicalEntity,
/*verbose=*/false, /*insertBoundarySegments=*/true);
//translate physicalTags from vector<int> to vector<vector<int>>
creator_->expandTagsVectors(elementToPhysicalEntity, boundarySegmentToPhysicalEntity);
}else{
GMSH4_ASSERT_MSG(file_type == 0 || file_type == 1,
"Invalid file-type: 0 for ASCII mode, 1 for binary mode");
GMSH4_ASSERT_MSG(data_size >= 4 && data_size <= 16,
"Invalid data-size range: should be in {4, 16}");
GMSH4_ASSERT_MSG(file_type != 1 || data_size == sizeof(size_type),
"Invalid data-size: must be sizeof(size_t)");
}
break;
}
case Sections::PHYSICAL_NAMES:
......@@ -345,6 +445,7 @@ void Gmsh4Reader<G,C,S>::readSerialFileFromStream (std::ifstream& input, bool fi
else readNodesBinary(input);
break;
case Sections::ELEMENTS:
prepareBoundingEntityTags();
if(file_type == 0) readElementsAscii(input);
else readElementsBinary(input);
break;
......@@ -376,7 +477,6 @@ void Gmsh4Reader<G,C,S>::readSerialFileFromStream (std::ifstream& input, bool fi
template <class G, class C, class S>
void Gmsh4Reader<G,C,S>::readParallelFileFromStream (std::ifstream& input, int commRank, int commSize, bool fillCreator)
{
clear();
DUNE_THROW(Dune::NotImplemented, "Reading parallel .pro files not yet implemented.");
if (fillCreator)
fillGridCreator();
......@@ -437,7 +537,6 @@ void Gmsh4Reader<G,C,S>::readEntitiesAscii (std::ifstream& input)
}
// points
points_.reserve(numPoints);
for (size_type i = 0; i < numPoints; ++i) {
if (!std::getline(input,line))
break;
......@@ -453,12 +552,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesAscii (std::ifstream& input)
for (size_type j = 0; j < numPhysicalTags; ++j)
stream >> attr.physicalTags[j];
points_.push_back(attr);
points_[attr.tag] = attr;
}
GMSH4_ASSERT(points_.size() == numPoints);
// curves
curves_.reserve(numCurves);
for (size_type i = 0; i < numCurves; ++i) {
if (!std::getline(input,line))
break;
......@@ -481,12 +579,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesAscii (std::ifstream& input)
for (size_type j = 0; j < numBoundingPoints; ++j)
stream >> attr.boundingEntities[j];
curves_.push_back(attr);
curves_[attr.tag] = attr;
}
GMSH4_ASSERT(curves_.size() == numCurves);
// surfaces
surfaces_.reserve(numSurfaces);
for (size_type i = 0; i < numSurfaces; ++i) {
if (!std::getline(input,line))
break;
......@@ -509,12 +606,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesAscii (std::ifstream& input)
for (size_type j = 0; j < numBoundingCurves; ++j)
stream >> attr.boundingEntities[j];
surfaces_.push_back(attr);
surfaces_[attr.tag] = attr;
}
GMSH4_ASSERT(surfaces_.size() == numSurfaces);
// volumes
volumes_.reserve(numVolumes);
for (size_type i = 0; i < numVolumes; ++i) {
if (!std::getline(input,line))
break;
......@@ -537,7 +633,7 @@ void Gmsh4Reader<G,C,S>::readEntitiesAscii (std::ifstream& input)
for (size_type j = 0; j < numBoundingSurfaces; ++j)
stream >> attr.boundingEntities[j];
volumes_.push_back(attr);
volumes_[attr.tag] = attr;
}
GMSH4_ASSERT(volumes_.size() == numVolumes);
}
......@@ -554,7 +650,6 @@ void Gmsh4Reader<G,C,S>::readEntitiesBinary (std::ifstream& input)
readValueBinary(input, numVolumes);
// points
points_.reserve(numPoints);
for (size_type i = 0; i < numPoints; ++i) {
PointAttributes attr;
......@@ -569,12 +664,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesBinary (std::ifstream& input)
for (size_type j = 0; j < numPhysicalTags; ++j)
readValueBinary(input, attr.physicalTags[j]);
points_.push_back(attr);
points_[attr.tag] = attr;
}
GMSH4_ASSERT(points_.size() == numPoints);
// curves
curves_.reserve(numCurves);
for (size_type i = 0; i < numCurves; ++i) {
EntityAttributes attr;
......@@ -598,12 +692,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesBinary (std::ifstream& input)
for (size_type j = 0; j < numBoundingPoints; ++j)
readValueBinary(input, attr.boundingEntities[j]);
curves_.push_back(attr);
curves_[attr.tag] = attr;
}
GMSH4_ASSERT(curves_.size() == numCurves);
// surfaces
surfaces_.reserve(numSurfaces);
for (size_type i = 0; i < numSurfaces; ++i) {
EntityAttributes attr;
......@@ -627,12 +720,11 @@ void Gmsh4Reader<G,C,S>::readEntitiesBinary (std::ifstream& input)
for (size_type j = 0; j < numBoundingCurves; ++j)
readValueBinary(input, attr.boundingEntities[j]);
surfaces_.push_back(attr);
surfaces_[attr.tag] = attr;
}
GMSH4_ASSERT(surfaces_.size() == numSurfaces);
// volumes
volumes_.reserve(numVolumes);
for (size_type i = 0; i < numVolumes; ++i) {
EntityAttributes attr;
......@@ -656,7 +748,7 @@ void Gmsh4Reader<G,C,S>::readEntitiesBinary (std::ifstream& input)
for (size_type j = 0; j < numBoundingSurfaces; ++j)
readValueBinary(input, attr.boundingEntities[j]);
volumes_.push_back(attr);
volumes_[attr.tag] = attr;
}
GMSH4_ASSERT(volumes_.size() == numVolumes);
std::string line;
......@@ -1123,6 +1215,14 @@ void Gmsh4Reader<G,C,S>::readElementsAscii (std::ifstream& input)
std::istringstream stream(line);
stream >> entityBlock.entityDim >> entityBlock.entityTag >> entityBlock.elementType >> numElementsInBlock;
if(entityBlock.entityDim == G::dimension){
entityBlock.physicalTags = getPhysicalTags(entityBlock.entityDim, entityBlock.entityTag);
}else if(entityBlock.entityDim == G::dimension - 1
&& boundingEntityTags.find(entityBlock.entityTag) != boundingEntityTags.end()){
boundaryEntities.insert(entityBlock.entityTag);
entityBlock.physicalTags = getPhysicalTags(entityBlock.entityDim, entityBlock.entityTag);
}
size_type numNodes = elementType_[entityBlock.elementType];
entityBlock.elements.resize(numElementsInBlock);
......@@ -1162,6 +1262,14 @@ void Gmsh4Reader<G,C,S>::readElementsBinary (std::ifstream& input)
readValueBinary(input, entityBlock.elementType);
readValueBinary(input, numElementsInBlock);
if(entityBlock.entityDim == G::dimension){
entityBlock.physicalTags = getPhysicalTags(entityBlock.entityDim, entityBlock.entityTag);
}else if(entityBlock.entityDim == G::dimension - 1
&& boundingEntityTags.find(entityBlock.entityTag) != boundingEntityTags.end()){
boundaryEntities.insert(entityBlock.entityTag);
entityBlock.physicalTags = getPhysicalTags(entityBlock.entityDim, entityBlock.entityTag);
}
size_type numNodes = elementType_[entityBlock.elementType];
entityBlock.elements.resize(numElementsInBlock);
......@@ -1227,7 +1335,6 @@ void Gmsh4Reader<G,C,S>::fillGridCreator (bool insertPieces)
if (!nodes_.empty())
creator_->insertVertices(numNodes_, {minNodeTag_, maxNodeTag_}, nodes_);
if (!elements_.empty()) {
std::set<int> boundaryEntities;
creator_->insertElements(numElements_, {minElementTag_, maxElementTag_}, elements_, boundaryEntities);
}
if (insertPieces)
......
......@@ -26,11 +26,22 @@ namespace Dune
using GlobalCoordinate = typename Grid::template Codim<0>::Entity::Geometry::GlobalCoordinate;
public:
/// Constructor. Creates an internal GridFactory, that will be deleted in destructor
GridCreatorInterface()
: internalFactory_(true), factory_(new GridFactory<Grid>)
{}
/// Constructor. Stores a reference to the passed GridFactory
GridCreatorInterface (GridFactory<Grid>& factory)
: factory_(&factory)
{}
/// Destructor.
virtual ~GridCreatorInterface()
{
if(internalFactory_) delete factory_;
}
/// Insert all points as vertices into the factory
template <class NodeAttributes>
void insertVertices (std::size_t numNodes,
......@@ -110,8 +121,49 @@ namespace Dune
/* do nothing */;
}
/** \brief obtain a vector of physical-tags for each element
*
* In order to determine the physical tags for a particular element use
* its insertion-index as index into elementTags. The insertion-index can
* be obtained e.g. by
* GridFactory::insertionIndex(const typename Codim< 0 >::Entity &entity)
*
* \returns vector containing for each element a vector of physical-tags
*/
std::vector<std::vector<int>>& elementTags()
{
return elementTags_;
}
/** \brief obtain a vector of physical-tags for each boundary-segment
*
* In order to determine the physical tags for a particular boundary-segment use
* its insertion-index as index into boundaryTags. The insertion-index can be
* obtained e.g. by
* GridFactory::insertionIndex(const typename GridType::LeafIntersection &intersection)
*
* \returns vector containing for each boundary-segment a vector of physical-tags
*/
std::vector<std::vector<int>>& boundaryTags()
{
return boundaryTags_;
}
void expandTagsVectors(std::vector<int> elementTags, std::vector<int> boundaryTags)
{
elementTags_.clear();
elementTags_.reserve(elementTags.size());
for(auto tag : elementTags) elementTags_.emplace_back(1, tag);
boundaryTags_.clear();
boundaryTags_.reserve(boundaryTags.size());
for(auto tag : boundaryTags) boundaryTags_.emplace_back(1, tag);
}
protected:
bool internalFactory_ = false;
GridFactory<Grid>* factory_;
std::vector<std::vector<int>> elementTags_;
std::vector<std::vector<int>> boundaryTags_;
};
} // end namespace Gmsh4
......
......@@ -31,6 +31,8 @@ namespace Dune
public:
using Super::Super;
using Super::factory;
using Super::boundaryTags_;
using Super::elementTags_;
template <class NodeAttributes>
void insertVerticesImpl (std::size_t numNodes,
......@@ -64,15 +66,21 @@ namespace Dune
std::vector<std::int64_t> cornerVertices(nodes_.size(), -1);
for (auto const& entityBlock : entityBlocks) {
if (entityBlock.entityDim < Grid::dimension-1)
bool boundaryDimension = (entityBlock.entityDim == Grid::dimension - 1);
if (entityBlock.entityDim != Grid::dimension && (!boundaryDimension
|| boundaryEntities.find(entityBlock.entityTag) == boundaryEntities.end()))
continue;
auto type = Gmsh4::to_geometry(entityBlock.elementType);
Gmsh4::CellType cell{type};
if (entityBlock.entityDim == Grid::dimension) { //element
if(boundaryDimension){
auto refElem = referenceElement<double,Grid::dimension - 1>(cell.type());
connectivity.resize(refElem.size(Grid::dimension - 1));
}else{
auto refElem = referenceElement<double,Grid::dimension>(cell.type());
connectivity.resize(refElem.size(Grid::dimension));
}
for (auto const& element : entityBlock.elements) {
GMSH4_ASSERT(element.nodes.size() >= connectivity.size());
......@@ -86,7 +94,12 @@ namespace Dune
connectivity[cell.permutation(j)] = vertex;
}
if(boundaryDimension){
factory().insertBoundarySegment(connectivity);
boundaryTags_.push_back(entityBlock.physicalTags);
}else{
factory().insertElement(cell.type(), connectivity);
elementTags_.push_back(entityBlock.physicalTags);
}
}
}
......@@ -94,7 +107,7 @@ namespace Dune
}
private:
/// All point coordinates inclusing the higher-order lagrange points
/// All point coordinates including the higher-order lagrange points
Nodes nodes_;
std::vector<std::size_t> vertexMap_;
......
......@@ -61,6 +61,8 @@ namespace Dune
public:
using Super::Super;
using Super::factory;
using Super::boundaryTags_;
using Super::elementTags_;
/// Implementation of the interface function `insertVertices()`
template <class NodeAttributes>
......@@ -101,15 +103,21 @@ namespace Dune
std::vector<std::int64_t> cornerVertices(nodes_.size(), -1);
for (auto const& entityBlock : entityBlocks) {
if (entityBlock.entityDim < dim - 1)
bool boundaryDimension = (entityBlock.entityDim == GridType::dimension - 1);
if (entityBlock.entityDim != GridType::dimension && (!boundaryDimension
|| boundaryEntities.find(entityBlock.entityTag) == boundaryEntities.end()))
continue;
auto type = Gmsh4::to_geometry(entityBlock.elementType);
Gmsh4::CellType cell{type};
if (entityBlock.entityDim == dim) { //element
auto refElem = referenceElement<double, dim>(type);
if(boundaryDimension){
auto refElem = referenceElement<double,dim - 1>(cell.type());
connectivity.resize(refElem.size(dim - 1));
}else{
auto refElem = referenceElement<double,dim>(cell.type());
connectivity.resize(refElem.size(dim));
}