3 #ifndef DUNE_GRID_IO_FILE_GMSHWRITER_HH
4 #define DUNE_GRID_IO_FILE_GMSHWRITER_HH
12 #include <dune/common/exceptions.hh>
13 #include <dune/geometry/type.hh>
14 #include <dune/geometry/referenceelements.hh>
33 template <
class Gr
idView>
42 static_assert( (dimWorld <= 3),
"GmshWriter requires dimWorld <= 3." );
45 template<
typename Entity>
46 std::size_t nodeIndexFromEntity(
const Entity& entity,
int i)
const {
47 return gv.
indexSet().subIndex(entity, i, dim)+1;
53 static std::size_t translateDuneToGmshType(
const GeometryType& type) {
54 std::size_t element_type;
58 else if (type.isTriangle())
60 else if (type.isQuadrilateral())
62 else if (type.isTetrahedron())
64 else if (type.isHexahedron())
66 else if (type.isPrism())
68 else if (type.isPyramid())
70 else if (type.isVertex())
73 DUNE_THROW(Dune::IOError,
"GeometryType " << type <<
" is not supported by gmsh.");
92 void outputElements(std::ofstream& file,
const std::vector<int>& physicalEntities,
const std::vector<int>& physicalBoundaries)
const {
94 std::size_t counter(1);
95 for (
const auto& entity : elements(gv)) {
98 std::size_t element_type = translateDuneToGmshType(entity.
type());
100 file << counter <<
" " << element_type;
102 if (!physicalEntities.empty())
103 file <<
" " << 1 <<
" " << physicalEntities[elementMapper.
index(entity)];
109 if (3 == element_type)
111 << nodeIndexFromEntity(entity, 0) <<
" " << nodeIndexFromEntity(entity, 1) <<
" "
112 << nodeIndexFromEntity(entity, 3) <<
" " << nodeIndexFromEntity(entity, 2);
113 else if (5 == element_type)
115 << nodeIndexFromEntity(entity, 0) <<
" " << nodeIndexFromEntity(entity, 1) <<
" "
116 << nodeIndexFromEntity(entity, 3) <<
" " << nodeIndexFromEntity(entity, 2) <<
" "
117 << nodeIndexFromEntity(entity, 4) <<
" " << nodeIndexFromEntity(entity, 5) <<
" "
118 << nodeIndexFromEntity(entity, 7) <<
" " << nodeIndexFromEntity(entity, 6);
119 else if (7 == element_type)
121 << nodeIndexFromEntity(entity, 0) <<
" " << nodeIndexFromEntity(entity, 1) <<
" "
122 << nodeIndexFromEntity(entity, 3) <<
" " << nodeIndexFromEntity(entity, 2) <<
" "
123 << nodeIndexFromEntity(entity, 4);
125 for (
int k = 0; k < entity.
geometry().corners(); ++k)
126 file <<
" " << nodeIndexFromEntity(entity, k);
133 if (!physicalBoundaries.empty()) {
134 const auto& refElement(ReferenceElements<typename GridView::ctype,dim>::general(entity.
type()));
135 for(
const auto& intersection : intersections(gv, entity)) {
136 if(intersection.boundary()) {
137 const auto faceLocalIndex(intersection.indexInInside());
138 file << counter <<
" " << translateDuneToGmshType(intersection.type())
139 <<
" " << 1 <<
" " << physicalBoundaries[intersection.boundarySegmentIndex()];
140 for (
int k = 0; k < intersection.geometry().corners(); ++k)
142 const auto vtxLocalIndex(refElement.subEntity(faceLocalIndex, 1, k, dim));
143 file <<
" " << nodeIndexFromEntity(entity, vtxLocalIndex);
151 }
catch(Exception& e) {
165 void outputNodes(std::ofstream& file)
const {
166 for (
const auto&
vertex : vertices(gv)) {
167 const auto globalCoord =
vertex.geometry().center();
171 file << nodeIndex <<
" " << globalCoord[0] <<
" " << 0 <<
" " << 0 << std::endl;
172 else if (2 == dimWorld)
173 file << nodeIndex <<
" " << globalCoord[0] <<
" " << globalCoord[1] <<
" " << 0 << std::endl;
175 file << nodeIndex <<
" " << globalCoord[0] <<
" " << globalCoord[1] <<
" " << globalCoord[2] << std::endl;
192 precision = numDigits;
216 void write(
const std::string& fileName,
217 const std::vector<int>& physicalEntities=std::vector<int>(),
218 const std::vector<int>& physicalBoundaries=std::vector<int>())
const {
220 std::ofstream file(fileName.c_str());
222 DUNE_THROW(Dune::IOError,
"Could not open " << fileName <<
" with write access.");
225 file << std::setprecision( precision );
228 file <<
"$MeshFormat" << std::endl
229 <<
"2.0 0 " <<
sizeof(double) << std::endl
230 <<
"$EndMeshFormat" << std::endl;
233 file <<
"$Nodes" << std::endl
234 << gv.
size(dim) << std::endl;
238 file <<
"$EndNodes" << std::endl;
241 int boundariesSize(0);
242 if(!physicalBoundaries.empty())
243 for(
const auto& entity : elements(gv))
244 for(
const auto& intersection : intersections(gv, entity))
245 if(intersection.boundary())
248 file <<
"$Elements" << std::endl
249 << gv.
size(0) + boundariesSize<< std::endl;
251 outputElements(file, physicalEntities, physicalBoundaries);
253 file <<
"$EndElements" << std::endl;
260 #endif // DUNE_GRID_IO_FILE_GMSHWRITER_HH