DUNE PDELab (git)

gmshwriter.hh
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_GRID_IO_FILE_GMSHWRITER_HH
6#define DUNE_GRID_IO_FILE_GMSHWRITER_HH
7
8#include <fstream>
9#include <iostream>
10#include <iomanip>
11#include <string>
12#include <vector>
13
15#include <dune/geometry/type.hh>
16#include <dune/geometry/referenceelements.hh>
19
20namespace Dune {
21
35 template <class GridView>
37 {
38 private:
39 const GridView gv;
40 int precision;
41
42 static const unsigned int dim = GridView::dimension;
43 static const unsigned int dimWorld = GridView::dimensionworld;
44 static_assert( (dimWorld <= 3), "GmshWriter requires dimWorld <= 3." );
45
47 template<typename Entity>
48 std::size_t nodeIndexFromEntity(const Entity& entity, int i) const {
49 return gv.indexSet().subIndex(entity, i, dim)+1;
50 }
51
55 static std::size_t translateDuneToGmshType(const GeometryType& type) {
56 std::size_t element_type;
57
58 if (type.isLine())
59 element_type = 1;
60 else if (type.isTriangle())
61 element_type = 2;
62 else if (type.isQuadrilateral())
63 element_type = 3;
64 else if (type.isTetrahedron())
65 element_type = 4;
66 else if (type.isHexahedron())
67 element_type = 5;
68 else if (type.isPrism())
69 element_type = 6;
70 else if (type.isPyramid())
71 element_type = 7;
72 else if (type.isVertex())
73 element_type = 15;
74 else
75 DUNE_THROW(Dune::IOError, "GeometryType " << type << " is not supported by gmsh.");
76
77 return element_type;
78 }
79
94 void outputElements(std::ofstream& file, const std::vector<int>& physicalEntities, const std::vector<int>& physicalBoundaries) const {
96 std::size_t counter(1);
97 for (const auto& entity : elements(gv)) {
98 // Check whether the type is compatible. If not, close file and rethrow exception.
99 try {
100 std::size_t element_type = translateDuneToGmshType(entity.type());
101
102 file << counter << " " << element_type;
103 // If present, set the first tag to the physical entity
104 if (!physicalEntities.empty())
105 file << " " << 1 << " " << physicalEntities[elementMapper.index(entity)];
106 else
107 file << " " << 0; // "0" for "I do not use any tags."
108
109 // Output list of nodes.
110 // 3, 5 and 7 got different vertex numbering compared to Dune
111 if (3 == element_type)
112 file << " "
113 << nodeIndexFromEntity(entity, 0) << " " << nodeIndexFromEntity(entity, 1) << " "
114 << nodeIndexFromEntity(entity, 3) << " " << nodeIndexFromEntity(entity, 2);
115 else if (5 == element_type)
116 file << " "
117 << nodeIndexFromEntity(entity, 0) << " " << nodeIndexFromEntity(entity, 1) << " "
118 << nodeIndexFromEntity(entity, 3) << " " << nodeIndexFromEntity(entity, 2) << " "
119 << nodeIndexFromEntity(entity, 4) << " " << nodeIndexFromEntity(entity, 5) << " "
120 << nodeIndexFromEntity(entity, 7) << " " << nodeIndexFromEntity(entity, 6);
121 else if (7 == element_type)
122 file << " "
123 << nodeIndexFromEntity(entity, 0) << " " << nodeIndexFromEntity(entity, 1) << " "
124 << nodeIndexFromEntity(entity, 3) << " " << nodeIndexFromEntity(entity, 2) << " "
125 << nodeIndexFromEntity(entity, 4);
126 else {
127 for (int k = 0; k < entity.geometry().corners(); ++k)
128 file << " " << nodeIndexFromEntity(entity, k);
129 }
130 ++counter;
131
132 file << std::endl;
133
134 // Write boundaries
135 if (!physicalBoundaries.empty()) {
136 auto refElement = referenceElement<typename GridView::ctype,dim>(entity.type());
137 for(const auto& intersection : intersections(gv, entity)) {
138 if(intersection.boundary()) {
139 const auto faceLocalIndex(intersection.indexInInside());
140 file << counter << " " << translateDuneToGmshType(intersection.type())
141 << " " << 1 << " " << physicalBoundaries[intersection.boundarySegmentIndex()];
142 for (int k = 0; k < intersection.geometry().corners(); ++k)
143 {
144 const auto vtxLocalIndex(refElement.subEntity(faceLocalIndex, 1, k, dim));
145 file << " " << nodeIndexFromEntity(entity, vtxLocalIndex);
146 }
147 ++counter;
148 file << std::endl;
149 }
150 }
151 }
152
153 } catch(Exception& e) {
154 file.close();
155 throw;
156 }
157 }
158 }
159
160
167 void outputNodes(std::ofstream& file) const {
168 for (const auto& vertex : vertices(gv)) {
169 const auto globalCoord = vertex.geometry().center();
170 const auto nodeIndex = gv.indexSet().index(vertex)+1; // Start counting indices by "1".
171
172 if (1 == dimWorld)
173 file << nodeIndex << " " << globalCoord[0] << " " << 0 << " " << 0 << std::endl;
174 else if (2 == dimWorld)
175 file << nodeIndex << " " << globalCoord[0] << " " << globalCoord[1] << " " << 0 << std::endl;
176 else // (3 == dimWorld)
177 file << nodeIndex << " " << globalCoord[0] << " " << globalCoord[1] << " " << globalCoord[2] << std::endl;
178 }
179 }
180
181 public:
187 GmshWriter(const GridView& gridView, int numDigits=6) : gv(gridView), precision(numDigits) {}
188
193 void setPrecision(int numDigits) {
194 precision = numDigits;
195 }
196
218 void write(const std::string& fileName,
219 const std::vector<int>& physicalEntities=std::vector<int>(),
220 const std::vector<int>& physicalBoundaries=std::vector<int>()) const {
221 // Open file
222 std::ofstream file(fileName.c_str());
223 if (!file.is_open())
224 DUNE_THROW(Dune::IOError, "Could not open " << fileName << " with write access.");
225
226 // Set precision
227 file << std::setprecision( precision );
228
229 // Output Header
230 file << "$MeshFormat" << std::endl
231 << "2.0 0 " << sizeof(double) << std::endl // "2.0" for "version 2.0", "0" for ASCII
232 << "$EndMeshFormat" << std::endl;
233
234 // Output Nodes
235 file << "$Nodes" << std::endl
236 << gv.size(dim) << std::endl;
237
238 outputNodes(file);
239
240 file << "$EndNodes" << std::endl;
241
242 // Output Elements;
243 int boundariesSize(0);
244 if(!physicalBoundaries.empty())
245 for(const auto& entity : elements(gv))
246 for(const auto& intersection : intersections(gv, entity))
247 if(intersection.boundary())
248 ++boundariesSize;
249
250 file << "$Elements" << std::endl
251 << gv.size(0) + boundariesSize<< std::endl;
252
253 outputElements(file, physicalEntities, physicalBoundaries);
254
255 file << "$EndElements" << std::endl;
256 }
257
258 };
259
260} // namespace Dune
261
262#endif // DUNE_GRID_IO_FILE_GMSHWRITER_HH
Wrapper class for entities.
Definition: entity.hh:66
Geometry geometry() const
obtain geometric realization of the entity
Definition: entity.hh:141
GeometryType type() const
Return the name of the reference element. The type can be used to access the Dune::ReferenceElement.
Definition: entity.hh:146
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Unique label for each type of entities that can occur in DUNE grids.
Definition: type.hh:114
constexpr bool isPyramid() const
Return true if entity is a pyramid.
Definition: type.hh:304
constexpr bool isTetrahedron() const
Return true if entity is a tetrahedron.
Definition: type.hh:299
constexpr bool isPrism() const
Return true if entity is a prism.
Definition: type.hh:309
constexpr bool isVertex() const
Return true if entity is a vertex.
Definition: type.hh:279
constexpr bool isTriangle() const
Return true if entity is a triangle.
Definition: type.hh:289
constexpr bool isLine() const
Return true if entity is a line segment.
Definition: type.hh:284
constexpr bool isQuadrilateral() const
Return true if entity is a quadrilateral.
Definition: type.hh:294
constexpr bool isHexahedron() const
Return true if entity is a hexahedron.
Definition: type.hh:314
Write Gmsh mesh file.
Definition: gmshwriter.hh:37
GmshWriter(const GridView &gridView, int numDigits=6)
Constructor expecting GridView of Grid to be written.
Definition: gmshwriter.hh:187
void setPrecision(int numDigits)
Set the number of digits to be used when writing the vertices. By default is 6.
Definition: gmshwriter.hh:193
void write(const std::string &fileName, const std::vector< int > &physicalEntities=std::vector< int >(), const std::vector< int > &physicalBoundaries=std::vector< int >()) const
Write given grid in Gmsh 2.0 compatible ASCII file.
Definition: gmshwriter.hh:218
Grid view abstract base class.
Definition: gridview.hh:66
Default exception class for I/O errors.
Definition: exceptions.hh:231
Index index(const EntityType &e) const
Map entity to starting index in array for dof block.
Definition: mcmgmapper.hh:171
Different resources needed by all grid implementations.
A few common exception classes.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
const IndexSet & indexSet() const
obtain the index set
Definition: gridview.hh:177
int size(int codim) const
obtain number of entities in a given codimension
Definition: gridview.hh:183
static constexpr int dimension
The dimension of the grid.
Definition: gridview.hh:134
static constexpr int dimensionworld
The dimension of the world the grid lives in.
Definition: gridview.hh:137
constexpr GeometryType vertex
GeometryType representing a vertex.
Definition: type.hh:492
MCMGLayout mcmgElementLayout()
layout for elements (codim-0 entities)
Definition: mcmgmapper.hh:97
Mapper for multiple codim and multiple geometry types.
Dune namespace.
Definition: alignedallocator.hh:13
A unique label for each type of element that can occur in a grid.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)