MFEM  v4.3.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-2021, 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 
22 {
23  POINT, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM
24 };
25 
27 {
28  POINT, QUADRATIC_SEGMENT, QUADRATIC_TRIANGLE, BIQUADRATIC_SQUARE,
29  QUADRATIC_TETRAHEDRON, TRIQUADRATIC_CUBE, BIQUADRATIC_QUADRATIC_PRISM
30 };
31 
33 {
34  POINT, LAGRANGE_SEGMENT, LAGRANGE_TRIANGLE, LAGRANGE_SQUARE,
35  LAGRANGE_TETRAHEDRON, LAGRANGE_CUBE, LAGRANGE_PRISM
36 };
37 
38 const int VTKGeometry::PrismMap[6] = {0, 2, 1, 3, 5, 4};
39 
41 {
42  NULL, NULL, NULL, NULL, NULL, NULL, VTKGeometry::PrismMap
43 };
44 
46 {
47  switch (vtk_geom)
48  {
49  case POINT:
50  return Geometry::POINT;
51  case SEGMENT:
52  case QUADRATIC_SEGMENT:
53  case LAGRANGE_SEGMENT:
54  return Geometry::SEGMENT;
55  case TRIANGLE:
56  case QUADRATIC_TRIANGLE:
57  case LAGRANGE_TRIANGLE:
58  return Geometry::TRIANGLE;
59  case SQUARE:
60  case BIQUADRATIC_SQUARE:
61  case LAGRANGE_SQUARE:
62  return Geometry::SQUARE;
63  case TETRAHEDRON:
66  return Geometry::TETRAHEDRON;
67  case CUBE:
68  case TRIQUADRATIC_CUBE:
69  case LAGRANGE_CUBE:
70  return Geometry::CUBE;
71  case PRISM:
73  case LAGRANGE_PRISM:
74  return Geometry::PRISM;
75  default:
76  return Geometry::INVALID;
77  }
78 }
79 
80 bool VTKGeometry::IsLagrange(int vtk_geom)
81 {
82  return vtk_geom >= LAGRANGE_SEGMENT && vtk_geom <= LAGRANGE_PRISM;
83 }
84 
85 bool VTKGeometry::IsQuadratic(int vtk_geom)
86 {
87  return vtk_geom >= QUADRATIC_SEGMENT
88  && vtk_geom <= BIQUADRATIC_QUADRATIC_PRISM;
89 }
90 
91 int VTKGeometry::GetOrder(int vtk_geom, int npoints)
92 {
93  if (IsQuadratic(vtk_geom))
94  {
95  return 2;
96  }
97  else if (IsLagrange(vtk_geom))
98  {
99  switch (vtk_geom)
100  {
101  case LAGRANGE_SEGMENT:
102  return npoints - 1;
103  case LAGRANGE_TRIANGLE:
104  return (std::sqrt(8*npoints + 1) - 3)/2;
105  case LAGRANGE_SQUARE:
106  return std::round(std::sqrt(npoints)) - 1;
108  switch (npoints)
109  {
110  // Note that for given order, npoints is given by
111  // npoints_order = (order + 1)*(order + 2)*(order + 3)/6,
112  case 4: return 1;
113  case 10: return 2;
114  case 20: return 3;
115  case 35: return 4;
116  case 56: return 5;
117  case 84: return 6;
118  case 120: return 7;
119  case 165: return 8;
120  case 220: return 9;
121  case 286: return 10;
122  default:
123  {
124  constexpr int max_order = 20;
125  int order = 11, npoints_order;
126  for (; order<max_order; ++order)
127  {
128  npoints_order = (order + 1)*(order + 2)*(order + 3)/6;
129  if (npoints_order == npoints) { break; }
130  }
131  MFEM_VERIFY(npoints == npoints_order, "");
132  return order;
133  }
134  }
135  case LAGRANGE_CUBE:
136  return std::round(std::cbrt(npoints)) - 1;
137  case LAGRANGE_PRISM:
138  {
139  const double n = npoints;
140  static const double third = 1.0/3.0;
141  static const double ninth = 1.0/9.0;
142  static const double twentyseventh = 1.0/27.0;
143  const double term =
144  std::cbrt(third*sqrt(third)*sqrt((27.0*n - 2.0)*n) + n
145  - twentyseventh);
146  return std::round(term + ninth / term - 4*third);
147  }
148  }
149  }
150  return 1;
151 }
152 
153 int BarycentricToVTKTriangle(int *b, int ref)
154 {
155  // Cf. https://git.io/JvW8f
156  int max = ref;
157  int min = 0;
158  int bmin = std::min(std::min(b[0], b[1]), b[2]);
159  int idx = 0;
160 
161  // scope into the correct triangle
162  while (bmin > min)
163  {
164  idx += 3*ref;
165  max -= 2;
166  ++min;
167  ref -= 3;
168  }
169  for (int d=0; d<3; ++d)
170  {
171  if (b[(d+2)%3] == max)
172  {
173  // we are on a vertex
174  return idx;
175  }
176  ++idx;
177  }
178  for (int d=0; d<3; ++d)
179  {
180  if (b[(d+1)%3] == min)
181  {
182  // we are on an edge
183  return idx + b[d] - (min + 1);
184  }
185  idx += max - (min + 1);
186  }
187  return idx;
188 }
189 
190 int BarycentricToVTKTetra(int *b, int ref)
191 {
192  // Cf. https://git.io/JvW8c
193  int idx = 0;
194 
195  int max = ref;
196  int min = 0;
197 
198  int bmin = std::min(std::min(std::min(b[0], b[1]), b[2]), b[3]);
199 
200  // scope into the correct tetra
201  while (bmin > min)
202  {
203  idx += 2*(ref*ref + 1);
204  max -= 3;
205  min++;
206  ref -= 4;
207  }
208 
209  // When a linearized tetra vertex is cast into barycentric coordinates, one of
210  // its coordinates is maximal and the other three are minimal. These are the
211  // indices of the maximal barycentric coordinate for each vertex.
212  static const int VertexMaxCoords[4] = {3,0,1,2};
213  // Each linearized tetra edge holds two barycentric tetra coordinates constant
214  // and varies the other two. These are the coordinates that are held constant
215  // for each edge.
216  static const int EdgeMinCoords[6][2] = {{1,2},{2,3},{0,2}, {0,1},{1,3},{0,3}};
217  // The coordinate that increments when traversing an edge (i.e. the coordinate
218  // of the nonzero component of the second vertex of the edge).
219  static const int EdgeCountingCoord[6] = {0,1,3,2,2,2};
220  // When describing a linearized tetra face, there is a mapping between the
221  // four-component barycentric tetra system and the three-component barycentric
222  // triangle system. These are the constant indices within the four-component
223  // system for each face (e.g. face 0 holds barycentric tetra coordinate 1
224  // constant).
225  static const int FaceMinCoord[4] = {1,3,0,2};
226  // When describing a linearized tetra face, there is a mapping between the
227  // four-component barycentric tetra system and the three-component barycentric
228  // triangle system. These are the relevant indices within the four-component
229  // system for each face (e.g. face 0 varies across the barycentric tetra
230  // coordinates 0, 2 and 3).
231  static const int FaceBCoords[4][3] = {{0,2,3}, {2,0,1}, {2,1,3}, {1,0,3}};
232 
233 
234  for (int vertex = 0; vertex < 4; vertex++)
235  {
236  if (b[VertexMaxCoords[vertex]] == max)
237  {
238  // we are on a vertex
239  return idx;
240  }
241  idx++;
242  }
243 
244  for (int edge = 0; edge < 6; edge++)
245  {
246  if (b[EdgeMinCoords[edge][0]] == min && b[EdgeMinCoords[edge][1]] == min)
247  {
248  // we are on an edge
249  return idx + b[EdgeCountingCoord[edge]] - (min + 1);
250  }
251  idx += max - (min + 1);
252  }
253 
254  for (int face = 0; face < 4; face++)
255  {
256  if (b[FaceMinCoord[face]] == min)
257  {
258  // we are on a face
259  int projectedb[3];
260  for (int i = 0; i < 3; i++)
261  {
262  projectedb[i] = b[FaceBCoords[face][i]] - min;
263  }
264  // we must subtract the indices of the face's vertices and edges, which
265  // total to 3*ref
266  return (idx + BarycentricToVTKTriangle(projectedb, ref) - 3*ref);
267  }
268  idx += (ref+1)*(ref+2)/2 - 3*ref;
269  }
270  return idx;
271 }
272 
273 int VTKTriangleDOFOffset(int ref, int i, int j)
274 {
275  return i + ref*(j - 1) - (j*(j + 1))/2;
276 }
277 
278 int CartesianToVTKPrism(int i, int j, int k, int ref)
279 {
280  // Cf. https://git.io/JvW0M
281  int om1 = ref - 1;
282  int ibdr = (i == 0);
283  int jbdr = (j == 0);
284  int ijbdr = (i + j == ref);
285  int kbdr = (k == 0 || k == ref);
286  // How many boundaries do we lie on at once?
287  int nbdr = ibdr + jbdr + ijbdr + kbdr;
288 
289  // Return an invalid index given invalid coordinates
290  if (i < 0 || i > ref || j < 0 || j > ref || i + j > ref || k < 0 || k > ref)
291  {
292  MFEM_ABORT("Invalid index")
293  }
294 
295  if (nbdr == 3) // Vertex DOF
296  {
297  // ijk is a corner node. Return the proper index (somewhere in [0,5]):
298  return (ibdr && jbdr ? 0 : (jbdr && ijbdr ? 1 : 2)) + (k ? 3 : 0);
299  }
300 
301  int offset = 6;
302  if (nbdr == 2) // Edge DOF
303  {
304  if (!kbdr)
305  {
306  // Must be on a vertical edge and 2 of {ibdr, jbdr, ijbdr} are true
307  offset += om1*6;
308  return offset + (k-1)
309  + ((ibdr && jbdr) ? 0 : (jbdr && ijbdr ? 1 : 2))*om1;
310  }
311  else
312  {
313  // Must be on a horizontal edge and kbdr plus 1 of {ibdr, jbdr, ijbdr} is true
314  // Skip past first 3 edges if we are on the top (k = ref) face:
315  offset += (k == ref ? 3*om1 : 0);
316  if (jbdr)
317  {
318  return offset + i - 1;
319  }
320  offset += om1; // Skip the i-axis edge
321  if (ijbdr)
322  {
323  return offset + j - 1;
324  }
325  offset += om1; // Skip the ij-axis edge
326  // if (ibdr)
327  return offset + (ref - j - 1);
328  }
329  }
330 
331  offset += 9*om1; // Skip all the edges
332 
333  // Number of points on a triangular face (but not on edge/corner):
334  int ntfdof = (om1 - 1)*om1/2;
335  int nqfdof = om1*om1;
336  if (nbdr == 1) // Face DOF
337  {
338  if (kbdr)
339  {
340  // We are on a triangular face.
341  if (k > 0)
342  {
343  offset += ntfdof;
344  }
345  return offset + VTKTriangleDOFOffset(ref, i, j);
346  }
347  // Not a k-normal face, so skip them:
348  offset += 2*ntfdof;
349 
350  // Face is quadrilateral (ref - 1) x (ref - 1)
351  // First face is i-normal, then ij-normal, then j-normal
352  if (jbdr) // On i-normal face
353  {
354  return offset + (i - 1) + om1*(k - 1);
355  }
356  offset += nqfdof; // Skip i-normal face
357  if (ijbdr) // on ij-normal face
358  {
359  return offset + (ref - i - 1) + om1*(k - 1);
360  }
361  offset += nqfdof; // Skip ij-normal face
362  return offset + j - 1 + om1*(k - 1);
363  }
364 
365  // Skip all face DOF
366  offset += 2*ntfdof + 3*nqfdof;
367 
368  // nbdr == 0: Body DOF
369  return offset + VTKTriangleDOFOffset(ref, i, j) + ntfdof*(k - 1);
370  // (i - 1) + (ref-1)*((j - 1) + (ref - 1)*(k - 1)));
371 }
372 
373 int CartesianToVTKTensor(int idx_in, int ref, Geometry::Type geom)
374 {
375  int n = ref + 1;
376  switch (geom)
377  {
378  case Geometry::POINT:
379  return idx_in;
380  case Geometry::SEGMENT:
381  if (idx_in == 0 || idx_in == ref)
382  {
383  return idx_in ? 1 : 0;
384  }
385  return idx_in + 1;
386  case Geometry::SQUARE:
387  {
388  // Cf: https://git.io/JvZLT
389  int i = idx_in % n;
390  int j = idx_in / n;
391  // Do we lie on any of the edges
392  bool ibdr = (i == 0 || i == ref);
393  bool jbdr = (j == 0 || j == ref);
394  if (ibdr && jbdr) // Vertex DOF
395  {
396  return (i ? (j ? 2 : 1) : (j ? 3 : 0));
397  }
398  int offset = 4;
399  if (jbdr) // Edge DOF on j==0 or j==ref
400  {
401  return (i - 1) + (j ? ref - 1 + ref - 1 : 0) + offset;
402  }
403  else if (ibdr) // Edge DOF on i==0 or i==ref
404  {
405  return (j - 1) + (i ? ref - 1 : 2 * (ref - 1) + ref - 1) + offset;
406  }
407  else // Interior DOF
408  {
409  offset += 2 * (ref - 1 + ref - 1);
410  return offset + (i - 1) + (ref - 1) * ((j - 1));
411  }
412  }
413  case Geometry::CUBE:
414  {
415  // Cf: https://git.io/JvZLe
416  int i = idx_in % n;
417  int j = (idx_in / n) % n;
418  int k = idx_in / (n*n);
419  bool ibdr = (i == 0 || i == ref);
420  bool jbdr = (j == 0 || j == ref);
421  bool kbdr = (k == 0 || k == ref);
422  // How many boundaries do we lie on at once?
423  int nbdr = (ibdr ? 1 : 0) + (jbdr ? 1 : 0) + (kbdr ? 1 : 0);
424  if (nbdr == 3) // Vertex DOF
425  {
426  // ijk is a corner node. Return the proper index (in [0,7])
427  return (i ? (j ? 2 : 1) : (j ? 3 : 0)) + (k ? 4 : 0);
428  }
429 
430  int offset = 8;
431  if (nbdr == 2) // Edge DOF
432  {
433  if (!ibdr)
434  {
435  // On i axis
436  return (i - 1) +
437  (j ? ref - 1 + ref - 1 : 0) +
438  (k ? 2*(ref - 1 + ref - 1) : 0) +
439  offset;
440  }
441  if (!jbdr)
442  {
443  // On j axis
444  return (j - 1) +
445  (i ? ref - 1 : 2*(ref - 1) + ref - 1) +
446  (k ? 2*(ref - 1 + ref - 1) : 0) +
447  offset;
448  }
449  // !kbdr, On k axis
450  offset += 4*(ref - 1) + 4*(ref - 1);
451  return (k - 1) + (ref - 1)*(i ? (j ? 3 : 1) : (j ? 2 : 0))
452  + offset;
453  }
454 
455  offset += 4*(ref - 1 + ref - 1 + ref - 1);
456  if (nbdr == 1) // Face DOF
457  {
458  if (ibdr) // On i-normal face
459  {
460  return (j - 1) + ((ref - 1)*(k - 1))
461  + (i ? (ref - 1)*(ref - 1) : 0) + offset;
462  }
463  offset += 2*(ref - 1)*(ref - 1);
464  if (jbdr) // On j-normal face
465  {
466  return (i - 1)
467  + ((ref - 1)*(k - 1))
468  + (j ? (ref - 1)*(ref - 1) : 0) + offset;
469  }
470  offset += 2*(ref - 1)*(ref - 1);
471  // kbdr, On k-normal face
472  return (i - 1) + ((ref - 1)*(j - 1))
473  + (k ? (ref - 1)*(ref - 1) : 0) + offset;
474  }
475 
476  // nbdr == 0: Interior DOF
477  offset += 2*((ref - 1)*(ref - 1) +
478  (ref - 1)*(ref - 1) +
479  (ref - 1)*(ref - 1));
480  return offset + (i - 1) + (ref - 1)*((j - 1) + (ref - 1)*(k - 1));
481  }
482  default:
483  MFEM_ABORT("CartesianToVTKOrderingTensor only supports tensor"
484  " geometries.");
485  return -1;
486  }
487 }
488 
490 {
491 
492  RefinedGeometry *RefG = GlobGeometryRefiner.Refine(geom, ref, 1);
493  int nnodes = RefG->RefPts.GetNPoints();
494  con.SetSize(nnodes);
495  if (geom == Geometry::TRIANGLE)
496  {
497  int b[3];
498  int idx = 0;
499  for (b[1]=0; b[1]<=ref; ++b[1])
500  {
501  for (b[0]=0; b[0]<=ref-b[1]; ++b[0])
502  {
503  b[2] = ref - b[0] - b[1];
504  con[BarycentricToVTKTriangle(b, ref)] = idx++;
505  }
506  }
507  }
508  else if (geom == Geometry::TETRAHEDRON)
509  {
510  int idx = 0;
511  int b[4];
512  for (b[2]=0; b[2]<=ref; b[2]++)
513  {
514  for (b[1]=0; b[1]<=ref-b[2]; b[1]++)
515  {
516  for (b[0]=0; b[0]<=ref-b[1]-b[2]; b[0]++)
517  {
518  b[3] = ref-b[0]-b[1]-b[2];
519  con[BarycentricToVTKTetra(b, ref)] = idx++;
520  }
521  }
522  }
523  }
524  else if (geom == Geometry::PRISM)
525  {
526  int idx = 0;
527  for (int k=0; k<=ref; k++)
528  {
529  for (int j=0; j<=ref; j++)
530  {
531  for (int i=0; i<=ref-j; i++)
532  {
533  con[CartesianToVTKPrism(i, j, k, ref)] = idx++;
534  }
535  }
536  }
537  }
538  else
539  {
540  for (int idx=0; idx<nnodes; ++idx)
541  {
542  con[CartesianToVTKTensor(idx, ref, geom)] = idx;
543  }
544  }
545 }
546 
547 void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes,
548  uint32_t nbytes, int compression_level)
549 {
550  if (compression_level == 0)
551  {
552  // First write size of buffer (as uint32_t), encoded with base 64
553  bin_io::WriteBase64(out, &nbytes, sizeof(nbytes));
554  // Then write all the bytes in the buffer, encoded with base 64
555  bin_io::WriteBase64(out, bytes, nbytes);
556  }
557  else
558  {
559 #ifdef MFEM_USE_ZLIB
560  MFEM_ASSERT(compression_level >= -1 && compression_level <= 9,
561  "Compression level must be between -1 and 9 (inclusive).");
562  uLongf buf_sz = compressBound(nbytes);
563  std::vector<unsigned char> buf(buf_sz);
564  compress2(buf.data(), &buf_sz, static_cast<const Bytef *>(bytes), nbytes,
565  compression_level);
566 
567  // Write the header
568  std::vector<uint32_t> header(4);
569  header[0] = 1; // number of blocks
570  header[1] = nbytes; // uncompressed size
571  header[2] = 0; // size of partial block
572  header[3] = buf_sz; // compressed size
573  bin_io::WriteBase64(out, header.data(), header.size()*sizeof(uint32_t));
574  // Write the compressed data
575  bin_io::WriteBase64(out, buf.data(), buf_sz);
576 #else
577  MFEM_ABORT("MFEM must be compiled with ZLib support to output "
578  "compressed binary data.")
579 #endif
580  }
581 }
582 
584 {
585  int16_t x16 = 1;
586  int8_t *x8 = reinterpret_cast<int8_t *>(&x16);
587  return !*x8;
588 }
589 
590 const char *VTKByteOrder()
591 {
592  if (IsBigEndian())
593  {
594  return "BigEndian";
595  }
596  else
597  {
598  return "LittleEndian";
599  }
600 
601 }
602 
603 } // namespace mfem
static const int BIQUADRATIC_SQUARE
Definition: vtk.hpp:35
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
int CartesianToVTKPrism(int i, int j, int k, int ref)
Definition: vtk.cpp:278
static const int HighOrderMap[Geometry::NUM_GEOMETRIES]
Definition: vtk.hpp:53
static const int QUADRATIC_TETRAHEDRON
Definition: vtk.hpp:36
const Geometry::Type geom
Definition: ex1.cpp:40
bool IsBigEndian()
Definition: vtk.cpp:583
static const int SQUARE
Definition: vtk.hpp:28
static bool IsLagrange(int vtk_geom)
Definition: vtk.cpp:80
int BarycentricToVTKTetra(int *b, int ref)
Definition: vtk.cpp:190
static const int LAGRANGE_CUBE
Definition: vtk.hpp:45
static const int CUBE
Definition: vtk.hpp:30
static const int PRISM
Definition: vtk.hpp:31
int VTKTriangleDOFOffset(int ref, int i, int j)
Definition: vtk.cpp:273
static const int QUADRATIC_TRIANGLE
Definition: vtk.hpp:34
int BarycentricToVTKTriangle(int *b, int ref)
Definition: vtk.cpp:153
static const int TRIQUADRATIC_CUBE
Definition: vtk.hpp:37
double b
Definition: lissajous.cpp:42
void WriteBase64(std::ostream &out, const void *bytes, size_t nbytes)
Definition: binaryio.cpp:25
static const int Map[Geometry::NUM_GEOMETRIES]
Definition: vtk.hpp:51
static const int LAGRANGE_SEGMENT
Definition: vtk.hpp:41
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1518
static const int LAGRANGE_SQUARE
Definition: vtk.hpp:43
IntegrationRule RefPts
Definition: geom.hpp:262
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:942
static const int BIQUADRATIC_QUADRATIC_PRISM
Definition: vtk.hpp:39
static const int QuadraticMap[Geometry::NUM_GEOMETRIES]
Definition: vtk.hpp:52
int CartesianToVTKTensor(int idx_in, int ref, Geometry::Type geom)
Definition: vtk.cpp:373
static Geometry::Type GetMFEMGeometry(int vtk_geom)
Definition: vtk.cpp:45
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:674
static const int TRIANGLE
Definition: vtk.hpp:27
static const int LAGRANGE_TETRAHEDRON
Definition: vtk.hpp:44
static const int QUADRATIC_SEGMENT
Definition: vtk.hpp:33
static const int TETRAHEDRON
Definition: vtk.hpp:29
const char * VTKByteOrder()
Definition: vtk.cpp:590
static const int POINT
Definition: vtk.hpp:25
void WriteVTKEncodedCompressed(std::ostream &out, const void *bytes, uint32_t nbytes, int compression_level)
Definition: vtk.cpp:547
static const int LAGRANGE_TRIANGLE
Definition: vtk.hpp:42
static const int PrismMap[6]
Definition: vtk.hpp:48
void CreateVTKElementConnectivity(Array< int > &con, Geometry::Type geom, int ref)
Definition: vtk.cpp:489
static const int SEGMENT
Definition: vtk.hpp:26
static bool IsQuadratic(int vtk_geom)
Definition: vtk.cpp:85
static const int * VertexPermutation[Geometry::NUM_GEOMETRIES]
Definition: vtk.hpp:49
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
static int GetOrder(int vtk_geom, int npoints)
Definition: vtk.cpp:91
static const int LAGRANGE_PRISM
Definition: vtk.hpp:46