MFEM  v4.2.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
vtk.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "vtk.hpp"
13 #include "../general/binaryio.hpp"
14 #ifdef MFEM_USE_ZLIB
15 #include <zlib.h>
16 #endif
17 
18 namespace mfem
19 {
20 
21 int BarycentricToVTKTriangle(int *b, int ref)
22 {
23  // Cf. https://git.io/JvW8f
24  int max = ref;
25  int min = 0;
26  int bmin = std::min(std::min(b[0], b[1]), b[2]);
27  int idx = 0;
28 
29  // scope into the correct triangle
30  while (bmin > min)
31  {
32  idx += 3*ref;
33  max -= 2;
34  ++min;
35  ref -= 3;
36  }
37  for (int d=0; d<3; ++d)
38  {
39  if (b[(d+2)%3] == max)
40  {
41  // we are on a vertex
42  return idx;
43  }
44  ++idx;
45  }
46  for (int d=0; d<3; ++d)
47  {
48  if (b[(d+1)%3] == min)
49  {
50  // we are on an edge
51  return idx + b[d] - (min + 1);
52  }
53  idx += max - (min + 1);
54  }
55  return idx;
56 }
57 
58 int BarycentricToVTKTetra(int *b, int ref)
59 {
60  // Cf. https://git.io/JvW8c
61  int idx = 0;
62 
63  int max = ref;
64  int min = 0;
65 
66  int bmin = std::min(std::min(std::min(b[0], b[1]), b[2]), b[3]);
67 
68  // scope into the correct tetra
69  while (bmin > min)
70  {
71  idx += 2*(ref*ref + 1);
72  max -= 3;
73  min++;
74  ref -= 4;
75  }
76 
77  // When a linearized tetra vertex is cast into barycentric coordinates, one of
78  // its coordinates is maximal and the other three are minimal. These are the
79  // indices of the maximal barycentric coordinate for each vertex.
80  static const int VertexMaxCoords[4] = {3,0,1,2};
81  // Each linearized tetra edge holds two barycentric tetra coordinates constant
82  // and varies the other two. These are the coordinates that are held constant
83  // for each edge.
84  static const int EdgeMinCoords[6][2] = {{1,2},{2,3},{0,2}, {0,1},{1,3},{0,3}};
85  // The coordinate that increments when traversing an edge (i.e. the coordinate
86  // of the nonzero component of the second vertex of the edge).
87  static const int EdgeCountingCoord[6] = {0,1,3,2,2,2};
88  // When describing a linearized tetra face, there is a mapping between the
89  // four-component barycentric tetra system and the three-component barycentric
90  // triangle system. These are the constant indices within the four-component
91  // system for each face (e.g. face 0 holds barycentric tetra coordinate 1
92  // constant).
93  static const int FaceMinCoord[4] = {1,3,0,2};
94  // When describing a linearized tetra face, there is a mapping between the
95  // four-component barycentric tetra system and the three-component barycentric
96  // triangle system. These are the relevant indices within the four-component
97  // system for each face (e.g. face 0 varies across the barycentric tetra
98  // coordinates 0, 2 and 3).
99  static const int FaceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}};
100 
101 
102  for (int vertex = 0; vertex < 4; vertex++)
103  {
104  if (b[VertexMaxCoords[vertex]] == max)
105  {
106  // we are on a vertex
107  return idx;
108  }
109  idx++;
110  }
111 
112  for (int edge = 0; edge < 6; edge++)
113  {
114  if (b[EdgeMinCoords[edge][0]] == min && b[EdgeMinCoords[edge][1]] == min)
115  {
116  // we are on an edge
117  return idx + b[EdgeCountingCoord[edge]] - (min + 1);
118  }
119  idx += max - (min + 1);
120  }
121 
122  for (int face = 0; face < 4; face++)
123  {
124  if (b[FaceMinCoord[face]] == min)
125  {
126  // we are on a face
127  int projectedb[3];
128  for (int i = 0; i < 3; i++)
129  {
130  projectedb[i] = b[FaceBCoords[face][i]] - min;
131  }
132  // we must subtract the indices of the face's vertices and edges, which
133  // total to 3*ref
134  return (idx + BarycentricToVTKTriangle(projectedb, ref) - 3*ref);
135  }
136  idx += (ref+1)*(ref+2)/2 - 3*ref;
137  }
138  return idx;
139 }
140 
141 int VTKTriangleDOFOffset(int ref, int i, int j)
142 {
143  return i + ref*(j - 1) - (j*(j + 1))/2;
144 }
145 
146 int CartesianToVTKPrism(int i, int j, int k, int ref)
147 {
148  // Cf. https://git.io/JvW0M
149  int om1 = ref - 1;
150  int ibdr = (i == 0);
151  int jbdr = (j == 0);
152  int ijbdr = (i + j == ref);
153  int kbdr = (k == 0 || k == ref);
154  // How many boundaries do we lie on at once?
155  int nbdr = ibdr + jbdr + ijbdr + kbdr;
156 
157  // Return an invalid index given invalid coordinates
158  if (i < 0 || i > ref || j < 0 || j > ref || i + j > ref || k < 0 || k > ref)
159  {
160  MFEM_ABORT("Invalid index")
161  }
162 
163  if (nbdr == 3) // Vertex DOF
164  {
165  // ijk is a corner node. Return the proper index (somewhere in [0,5]):
166  return (ibdr && jbdr ? 0 : (jbdr && ijbdr ? 1 : 2)) + (k ? 3 : 0);
167  }
168 
169  int offset = 6;
170  if (nbdr == 2) // Edge DOF
171  {
172  if (!kbdr)
173  {
174  // Must be on a vertical edge and 2 of {ibdr, jbdr, ijbdr} are true
175  offset += om1*6;
176  return offset + (k-1)
177  + ((ibdr && jbdr) ? 0 : (jbdr && ijbdr ? 1 : 2))*om1;
178  }
179  else
180  {
181  // Must be on a horizontal edge and kbdr plus 1 of {ibdr, jbdr, ijbdr} is true
182  // Skip past first 3 edges if we are on the top (k = ref) face:
183  offset += (k == ref ? 3*om1 : 0);
184  if (jbdr)
185  {
186  return offset + i - 1;
187  }
188  offset += om1; // Skip the i-axis edge
189  if (ijbdr)
190  {
191  return offset + j - 1;
192  }
193  offset += om1; // Skip the ij-axis edge
194  // if (ibdr)
195  return offset + (ref - j - 1);
196  }
197  }
198 
199  offset += 9*om1; // Skip all the edges
200 
201  // Number of points on a triangular face (but not on edge/corner):
202  int ntfdof = (om1 - 1)*om1/2;
203  int nqfdof = om1*om1;
204  if (nbdr == 1) // Face DOF
205  {
206  if (kbdr)
207  {
208  // We are on a triangular face.
209  if (k > 0)
210  {
211  offset += ntfdof;
212  }
213  return offset + VTKTriangleDOFOffset(ref, i, j);
214  }
215  // Not a k-normal face, so skip them:
216  offset += 2*ntfdof;
217 
218  // Face is quadrilateral (ref - 1) x (ref - 1)
219  // First face is i-normal, then ij-normal, then j-normal
220  if (jbdr) // On i-normal face
221  {
222  return offset + (i - 1) + om1*(k - 1);
223  }
224  offset += nqfdof; // Skip i-normal face
225  if (ijbdr) // on ij-normal face
226  {
227  return offset + (ref - i - 1) + om1*(k - 1);
228  }
229  offset += nqfdof; // Skip ij-normal face
230  return offset + j - 1 + om1*(k - 1);
231  }
232 
233  // Skip all face DOF
234  offset += 2*ntfdof + 3*nqfdof;
235 
236  // nbdr == 0: Body DOF
237  return offset + VTKTriangleDOFOffset(ref, i, j) + ntfdof*(k - 1);
238  // (i - 1) + (ref-1)*((j - 1) + (ref - 1)*(k - 1)));
239 }
240 
241 int CartesianToVTKTensor(int idx_in, int ref, Geometry::Type geom)
242 {
243  int n = ref + 1;
244  switch (geom)
245  {
246  case Geometry::POINT:
247  return idx_in;
248  case Geometry::SEGMENT:
249  if (idx_in == 0 || idx_in == ref)
250  {
251  return idx_in ? 1 : 0;
252  }
253  return idx_in + 1;
254  case Geometry::SQUARE:
255  {
256  // Cf: https://git.io/JvZLT
257  int i = idx_in % n;
258  int j = idx_in / n;
259  // Do we lie on any of the edges
260  bool ibdr = (i == 0 || i == ref);
261  bool jbdr = (j == 0 || j == ref);
262  if (ibdr && jbdr) // Vertex DOF
263  {
264  return (i ? (j ? 2 : 1) : (j ? 3 : 0));
265  }
266  int offset = 4;
267  if (jbdr) // Edge DOF on j==0 or j==ref
268  {
269  return (i - 1) + (j ? ref - 1 + ref - 1 : 0) + offset;
270  }
271  else if (ibdr) // Edge DOF on i==0 or i==ref
272  {
273  return (j - 1) + (i ? ref - 1 : 2 * (ref - 1) + ref - 1) + offset;
274  }
275  else // Interior DOF
276  {
277  offset += 2 * (ref - 1 + ref - 1);
278  return offset + (i - 1) + (ref - 1) * ((j - 1));
279  }
280  }
281  case Geometry::CUBE:
282  {
283  // Cf: https://git.io/JvZLe
284  int i = idx_in % n;
285  int j = (idx_in / n) % n;
286  int k = idx_in / (n*n);
287  bool ibdr = (i == 0 || i == ref);
288  bool jbdr = (j == 0 || j == ref);
289  bool kbdr = (k == 0 || k == ref);
290  // How many boundaries do we lie on at once?
291  int nbdr = (ibdr ? 1 : 0) + (jbdr ? 1 : 0) + (kbdr ? 1 : 0);
292  if (nbdr == 3) // Vertex DOF
293  {
294  // ijk is a corner node. Return the proper index (in [0,7])
295  return (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
296  }
297 
298  int offset = 8;
299  if (nbdr == 2) // Edge DOF
300  {
301  if (!ibdr)
302  {
303  // On i axis
304  return (i - 1) +
305  (j ? ref - 1 + ref - 1 : 0) +
306  (k ? 2*(ref - 1 + ref - 1) : 0) +
307  offset;
308  }
309  if (!jbdr)
310  {
311  // On j axis
312  return (j - 1) +
313  (i ? ref - 1 : 2*(ref - 1) + ref - 1) +
314  (k ? 2*(ref - 1 + ref - 1) : 0) +
315  offset;
316  }
317  // !kbdr, On k axis
318  offset += 4*(ref - 1) + 4*(ref - 1);
319  return (k - 1) + (ref - 1)*(i ? (j ? 3 : 1) : (j ? 2 : 0))
320  + offset;
321  }
322 
323  offset += 4*(ref - 1 + ref - 1 + ref - 1);
324  if (nbdr == 1) // Face DOF
325  {
326  if (ibdr) // On i-normal face
327  {
328  return (j - 1) + ((ref - 1)*(k - 1))
329  + (i ? (ref - 1)*(ref - 1) : 0) + offset;
330  }
331  offset += 2*(ref - 1)*(ref - 1);
332  if (jbdr) // On j-normal face
333  {
334  return (i - 1)
335  + ((ref - 1)*(k - 1))
336  + (j ? (ref - 1)*(ref - 1) : 0) + offset;
337  }
338  offset += 2*(ref - 1)*(ref - 1);
339  // kbdr, On k-normal face
340  return (i - 1) + ((ref - 1)*(j - 1))
341  + (k ? (ref - 1)*(ref - 1) : 0) + offset;
342  }
343 
344  // nbdr == 0: Interior DOF
345  offset += 2*((ref - 1)*(ref - 1) +
346  (ref - 1)*(ref - 1) +
347  (ref - 1)*(ref - 1));
348  return offset + (i - 1) + (ref - 1)*((j - 1) + (ref - 1)*(k - 1));
349  }
350  default:
351  MFEM_ABORT("CartesianToVTKOrderingTensor only supports tensor"
352  " geometries.");
353  return -1;
354  }
355 }
356 
358 {
359 
360  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
361  int nnodes = RefG->RefPts.GetNPoints();
362  con.SetSize(nnodes);
363  if (geom == Geometry::TRIANGLE)
364  {
365  int b[3];
366  int idx = 0;
367  for (b[1]=0; b[1]<=ref; ++b[1])
368  {
369  for (b[0]=0; b[0]<=ref-b[1]; ++b[0])
370  {
371  b[2] = ref - b[0] - b[1];
372  con[BarycentricToVTKTriangle(b, ref)] = idx++;
373  }
374  }
375  }
376  else if (geom == Geometry::TETRAHEDRON)
377  {
378  int idx = 0;
379  int b[4];
380  for (int k=0; k<=ref; k++)
381  {
382  for (int j=0; j<=k; j++)
383  {
384  for (int i=0; i<=j; i++)
385  {
386  b[0] = k-j;
387  b[1] = i;
388  b[2] = j-i;
389  b[3] = ref-b[0]-b[1]-b[2];
390  con[BarycentricToVTKTetra(b, ref)] = idx++;
391  }
392  }
393  }
394  }
395  else if (geom == Geometry::PRISM)
396  {
397  int idx = 0;
398  for (int k=0; k<=ref; k++)
399  {
400  for (int j=0; j<=ref; j++)
401  {
402  for (int i=0; i<=ref-j; i++)
403  {
404  con[CartesianToVTKPrism(i, j, k, ref)] = idx++;
405  }
406  }
407  }
408  }
409  else
410  {
411  for (int idx=0; idx<nnodes; ++idx)
412  {
413  con[CartesianToVTKTensor(idx, ref, geom)] = idx;
414  }
415  }
416 }
417 
418 void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes,
419  uint32_t nbytes, int compression_level)
420 {
421  if (compression_level == 0)
422  {
423  // First write size of buffer (as uint32_t), encoded with base 64
424  bin_io::WriteBase64(out, &nbytes, sizeof(nbytes));
425  // Then write all the bytes in the buffer, encoded with base 64
426  bin_io::WriteBase64(out, bytes, nbytes);
427  }
428  else
429  {
430 #ifdef MFEM_USE_ZLIB
431  MFEM_ASSERT(compression_level >= -1 && compression_level <= 9,
432  "Compression level must be between -1 and 9 (inclusive).");
433  uLongf buf_sz = compressBound(nbytes);
434  std::vector<unsigned char> buf(buf_sz);
435  compress2(buf.data(), &buf_sz, static_cast<const Bytef *>(bytes), nbytes,
436  compression_level);
437 
438  // Write the header
439  std::vector<uint32_t> header(4);
440  header[0] = 1; // number of blocks
441  header[1] = nbytes; // uncompressed size
442  header[2] = 0; // size of partial block
443  header[3] = buf_sz; // compressed size
444  bin_io::WriteBase64(out, header.data(), header.size()*sizeof(uint32_t));
445  // Write the compressed data
446  bin_io::WriteBase64(out, buf.data(), buf_sz);
447 #else
448  MFEM_ABORT("MFEM must be compiled with ZLib support to output "
449  "compressed binary data.")
450 #endif
451  }
452 }
453 
455 {
456  int16_t x16 = 1;
457  int8_t *x8 = reinterpret_cast<int8_t *>(&x16);
458  return !*x8;
459 }
460 
461 const char *VTKByteOrder()
462 {
463  if (IsBigEndian())
464  {
465  return "BigEndian";
466  }
467  else
468  {
469  return "LittleEndian";
470  }
471 
472 }
473 
474 } // namespace mfem
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:245
int CartesianToVTKPrism(int i, int j, int k, int ref)
Definition: vtk.cpp:146
const Geometry::Type geom
Definition: ex1.cpp:40
bool IsBigEndian()
Definition: vtk.cpp:454
int BarycentricToVTKTetra(int *b, int ref)
Definition: vtk.cpp:58
int VTKTriangleDOFOffset(int ref, int i, int j)
Definition: vtk.cpp:141
int BarycentricToVTKTriangle(int *b, int ref)
Definition: vtk.cpp:21
double b
Definition: lissajous.cpp:42
void WriteBase64(std::ostream &out, const void *bytes, size_t nbytes)
Definition: binaryio.cpp:25
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1520
IntegrationRule RefPts
Definition: geom.hpp:243
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:942
int CartesianToVTKTensor(int idx_in, int ref, Geometry::Type geom)
Definition: vtk.cpp:241
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:654
const char * VTKByteOrder()
Definition: vtk.cpp:461
void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes, uint32_t nbytes, int compression_level)
Definition: vtk.cpp:418
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Definition: vtk.cpp:357
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66