preCICE v3.1.1
Loading...
Searching...
No Matches
ExportVTK.cpp
Go to the documentation of this file.
1#include "ExportVTK.hpp"
2#include <Eigen/Core>
3#include <filesystem>
4#include <fstream>
5#include <iomanip>
6#include <memory>
7#include "io/Export.hpp"
9#include "mesh/Data.hpp"
10#include "mesh/Edge.hpp"
11#include "mesh/Mesh.hpp"
13#include "mesh/Triangle.hpp"
14#include "mesh/Vertex.hpp"
15#include "utils/IntraComm.hpp"
16#include "utils/assertion.hpp"
17
18namespace precice::io {
19
21 const std::string &name,
22 const std::string &location,
23 const mesh::Mesh & mesh)
24{
25 PRECICE_TRACE(name, location, mesh.getName());
27 PRECICE_ASSERT(!utils::IntraComm::isParallel(), "ExportVTK only supports serial participants.");
28
29 namespace fs = std::filesystem;
30 fs::path outfile(location);
31 if (not location.empty())
33 outfile = outfile / fs::path(name + ".vtk");
34 std::ofstream outstream(outfile.string(), std::ios::trunc);
35 PRECICE_CHECK(outstream, "VTK export failed to open destination file \"{}\"", outfile.generic_string());
36
37 initializeWriting(outstream);
38 writeHeader(outstream);
39 exportMesh(outstream, mesh);
40 exportData(outstream, mesh);
41 exportGradient(outstream, mesh);
42 outstream.close();
43}
44
46 std::ofstream & outFile,
47 const mesh::Mesh &mesh)
48{
49 PRECICE_TRACE(mesh.getName());
50
51 // Plot vertices
52 outFile << "POINTS " << mesh.nVertices() << " double \n\n";
53 for (const mesh::Vertex &vertex : mesh.vertices()) {
54 writeVertex(vertex.getCoords(), outFile);
55 }
56 outFile << '\n';
57
58 // Plot edges
59 if (mesh.getDimensions() == 2) {
60 outFile << "CELLS " << mesh.edges().size() << ' ' << mesh.edges().size() * 3 << "\n\n";
61 for (auto const &edge : mesh.edges()) {
62 int internalIndices[2];
63 internalIndices[0] = edge.vertex(0).getID();
64 internalIndices[1] = edge.vertex(1).getID();
65 writeLine(internalIndices, outFile);
66 }
67 outFile << "\nCELL_TYPES " << mesh.edges().size() << "\n\n";
68 for (size_t i = 0; i < mesh.edges().size(); ++i) {
69 outFile << "3\n";
70 }
71 }
72
73 // Plot triangles
74 if (mesh.getDimensions() == 3) {
75 size_t sizeTetrahedra = mesh.tetrahedra().size();
76 size_t sizeTriangles = mesh.triangles().size();
77 size_t sizeEdges = mesh.edges().size();
78 size_t sizeElements = sizeTriangles + sizeEdges + sizeTetrahedra;
79
80 outFile << "CELLS " << sizeElements << ' '
81 << sizeTetrahedra * 5 + sizeTriangles * 4 + sizeEdges * 3 << "\n\n";
82 for (auto const &tetra : mesh.tetrahedra()) {
83 int internalIndices[4];
84 internalIndices[0] = tetra.vertex(0).getID();
85 internalIndices[1] = tetra.vertex(1).getID();
86 internalIndices[2] = tetra.vertex(2).getID();
87 internalIndices[3] = tetra.vertex(3).getID();
88 writeTetrahedron(internalIndices, outFile);
89 }
90 for (auto const &triangle : mesh.triangles()) {
91 int internalIndices[3];
92 internalIndices[0] = triangle.vertex(0).getID();
93 internalIndices[1] = triangle.vertex(1).getID();
94 internalIndices[2] = triangle.vertex(2).getID();
95 writeTriangle(internalIndices, outFile);
96 }
97 for (auto const &edge : mesh.edges()) {
98 int internalIndices[2];
99 internalIndices[0] = edge.vertex(0).getID();
100 internalIndices[1] = edge.vertex(1).getID();
101 writeLine(internalIndices, outFile);
102 }
103
104 outFile << "\nCELL_TYPES " << sizeElements << "\n\n";
105 // See VTK reference for CELL_TYPES
106 for (size_t i = 0; i < sizeTetrahedra; i++) {
107 outFile << "10\n";
108 }
109 for (size_t i = 0; i < sizeTriangles; i++) {
110 outFile << "5\n";
111 }
112 for (size_t i = 0; i < sizeEdges; ++i) {
113 outFile << "3\n";
114 }
115 }
116
117 outFile << '\n';
118}
119
121 std::ofstream & outFile,
122 const mesh::Mesh &mesh)
123{
124 outFile << "POINT_DATA " << mesh.nVertices() << "\n\n";
125
126 outFile << "SCALARS Rank unsigned_int\n";
127 outFile << "LOOKUP_TABLE default\n";
129 outFile << "\n\n";
130
131 for (const mesh::PtrData &data : mesh.data()) { // Plot vertex data
132 Eigen::VectorXd &values = data->values();
133 if (data->getDimensions() > 1) {
134 Eigen::VectorXd viewTemp(data->getDimensions());
135 outFile << "VECTORS " << data->getName() << " double\n";
136 for (const mesh::Vertex &vertex : mesh.vertices()) {
137 int offset = vertex.getID() * data->getDimensions();
138 for (int i = 0; i < data->getDimensions(); i++) {
139 viewTemp[i] = values(offset + i);
140 }
141 int i = 0;
142 for (; i < data->getDimensions(); i++) {
143 outFile << viewTemp[i] << ' ';
144 }
145 if (i < 3) {
146 outFile << '0';
147 }
148 outFile << '\n';
149 }
150 outFile << '\n';
151 } else if (data->getDimensions() == 1) {
152 outFile << "SCALARS " << data->getName() << " double\n";
153 outFile << "LOOKUP_TABLE default\n";
154 for (const mesh::Vertex &vertex : mesh.vertices()) {
155 outFile << values(vertex.getID()) << '\n';
156 }
157 outFile << '\n';
158 }
159 }
160}
161
163{
164 const int spaceDim = mesh.getDimensions();
165 for (const mesh::PtrData &data : mesh.data()) {
166 if (data->hasGradient()) { // Check whether this data has gradient
167 auto &gradients = data->gradients();
168 if (data->getDimensions() == 1) { // Scalar data, create a vector <dataname>_gradient
169 outFile << "VECTORS " << data->getName() << "_gradient"
170 << " double\n";
171 for (int i = 0; i < gradients.cols(); i++) { // Loop over vertices
172 int j = 0; // Dimension counter
173 for (; j < gradients.rows(); j++) { // Loop over space directions
174 outFile << gradients.coeff(j, i) << " ";
175 }
176 if (j < 3) { // If 2D data add additional zero as third component
177 outFile << '0';
178 }
179 outFile << "\n";
180 }
181 } else { // Vector data, write n vector for n dimension <dataname>_(dx/dy/dz)
182 outFile << "VECTORS " << data->getName() << "_dx"
183 << " double\n";
184 for (int i = 0; i < gradients.cols(); i += spaceDim) { // Loop over vertices
185 int j = 0;
186 for (; j < gradients.rows(); j++) { // Loop over components
187 outFile << gradients.coeff(j, i) << " ";
188 }
189 if (j < 3) { // If 2D data add additional zero as third component
190 outFile << '0';
191 }
192 outFile << "\n";
193 }
194 outFile << "\n";
195
196 outFile << "VECTORS " << data->getName() << "_dy"
197 << " double\n";
198 for (int i = 1; i < gradients.cols(); i += spaceDim) { // Loop over vertices
199 int j = 0;
200 for (; j < gradients.rows(); j++) { // Loop over components
201 outFile << gradients.coeff(j, i) << " ";
202 }
203 if (j < 3) { // If 2D data add additional zero as third component
204 outFile << '0';
205 }
206 outFile << "\n";
207 }
208 outFile << "\n";
209
210 if (spaceDim == 3) { // dz is only for 3D data
211 outFile << "VECTORS " << data->getName() << "_dz"
212 << " double\n";
213 for (int i = 2; i < gradients.cols(); i += spaceDim) { // Loop over vertices
214 int j = 0;
215 for (; j < gradients.rows(); j++) { // Loop over components
216 outFile << gradients.coeff(j, i) << " ";
217 }
218 if (j < 3) { // If 2D data add additional zero as third component
219 outFile << '0';
220 }
221 outFile << "\n";
222 }
223 }
224 }
225 outFile << '\n';
226 }
227 }
228}
229
231 std::ofstream &filestream)
232{
233 // size_t pos = fullFilename.rfind(".vtk");
234 // if ((pos == std::string::npos) || (pos != fullFilename.size()-4)){
235 // fullFilename += ".vtk";
236 // }
237 filestream.setf(std::ios::showpoint);
238 filestream.setf(std::ios::scientific);
240}
241
243 std::ostream &outFile)
244{
245 outFile << "# vtk DataFile Version 2.0\n\n"
246 << "ASCII\n\n"
247 << "DATASET UNSTRUCTURED_GRID\n\n";
248}
249
251 const Eigen::VectorXd &position,
252 std::ostream & outFile)
253{
254 if (position.size() == 2) {
255 outFile << position(0) << " " << position(1) << " " << 0.0 << '\n';
256 } else {
257 PRECICE_ASSERT(position.size() == 3);
258 outFile << position(0) << " " << position(1) << " " << position(2) << '\n';
259 }
260}
261
263 int vertexIndices[3],
264 std::ostream &outFile)
265{
266 outFile << 3 << ' ';
267 for (int i = 0; i < 3; i++) {
268 outFile << vertexIndices[i] << ' ';
269 }
270 outFile << '\n';
271}
272
274 int vertexIndices[4],
275 std::ostream &outFile)
276{
277 outFile << 4 << ' ';
278 for (int i = 0; i < 4; i++) {
279 outFile << vertexIndices[i] << ' ';
280 }
281 outFile << '\n';
282}
283
285 int vertexIndices[2],
286 std::ostream &outFile)
287{
288 outFile << 2 << ' ';
289 for (int i = 0; i < 2; i++) {
290 outFile << vertexIndices[i] << ' ';
291 }
292 outFile << '\n';
293}
294
295} // namespace precice::io
#define PRECICE_TRACE(...)
Definition LogMacros.hpp:95
#define PRECICE_CHECK(check,...)
Definition LogMacros.hpp:35
std::string name
#define PRECICE_ASSERT(...)
Definition assertion.hpp:87
void exportData(std::ofstream &outFile, const mesh::Mesh &mesh)
static void writeVertex(const Eigen::VectorXd &position, std::ostream &outFile)
static void writeHeader(std::ostream &outFile)
static void writeLine(int vertexIndices[2], std::ostream &outFile)
static void initializeWriting(std::ofstream &filestream)
static void writeTetrahedron(int vertexIndices[4], std::ostream &outFile)
void exportMesh(std::ofstream &outFile, const mesh::Mesh &mesh)
Definition ExportVTK.cpp:45
virtual void doExport(const std::string &name, const std::string &location, const mesh::Mesh &mesh)
Perform writing to VTK file.
Definition ExportVTK.cpp:20
static void writeTriangle(int vertexIndices[3], std::ostream &outFile)
void exportGradient(std::ofstream &outFile, const mesh::Mesh &mesh)
Container and creator for meshes.
Definition Mesh.hpp:39
int getDimensions() const
Definition Mesh.cpp:98
VertexContainer & vertices()
Returns modifieable container holding all vertices.
Definition Mesh.cpp:53
const std::string & getName() const
Returns the name of the mesh, as set in the config file.
Definition Mesh.cpp:218
std::size_t nVertices() const
Returns the number of vertices.
Definition Mesh.cpp:63
TetraContainer & tetrahedra()
Returns modifiable container holding all tetrahedra.
Definition Mesh.cpp:93
const DataContainer & data() const
Allows access to all data.
Definition Mesh.cpp:170
TriangleContainer & triangles()
Returns modifiable container holding all triangles.
Definition Mesh.cpp:78
EdgeContainer & edges()
Returns modifiable container holding all edges.
Definition Mesh.cpp:68
Vertex of a mesh.
Definition Vertex.hpp:16
static bool isParallel()
True if this process is running in parallel.
Definition IntraComm.cpp:62
T close(T... args)
T create_directories(T... args)
T empty(T... args)
T fill_n(T... args)
T generic_string(T... args)
provides Import and Export of the coupling mesh and data.
T setf(T... args)
T setprecision(T... args)
T size(T... args)