MFEM  v4.5.2
Finite element discretization library
reflector.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 // ------------------------------------------------------------------------
13 // Reflector Miniapp: Reflect a mesh about a plane
14 // ------------------------------------------------------------------------
15 //
16 // This miniapp reflects a 3D mesh about a plane defined by a point and a
17 // normal vector. Element and boundary element attributes are copied from the
18 // corresponding elements in the original mesh, except for boundary elements on
19 // the plane of reflection.
20 //
21 // Compile with: make reflector
22 //
23 // Sample runs: reflector -m ../../data/pipe-nurbs.mesh -n '0 0 1'
24 // reflector -m ../../data/fichera.mesh -o '1 0 0' -n '1 0 0'
25 
26 
27 #include "mfem.hpp"
28 #include <fstream>
29 #include <iostream>
30 #include <set>
31 #include <array>
32 
33 using namespace std;
34 using namespace mfem;
35 
36 
37 void ReflectPoint(Vector & p, Vector const& origin, Vector const& normal)
38 {
39  Vector diff(3);
40  Vector proj(3);
41  subtract(p, origin, diff);
42  const double ip = diff * normal;
43 
44  diff = normal;
45  diff *= -2.0 * ip;
46 
47  add(p, diff, p);
48 }
49 
50 class ReflectedCoefficient : public VectorCoefficient
51 {
52 private:
54  const Vector origin, normal;
55  Mesh *meshOrig;
56 
57  // Map from reflected to original mesh elements
58  std::vector<int> *r2o;
59 
60  std::vector<std::vector<int>> *perm;
61 
62 public:
63  ReflectedCoefficient(VectorCoefficient &A, Vector const& origin_,
64  Vector const& normal_, std::vector<int> *r2o_,
65  Mesh *mesh, std::vector<std::vector<int>> *refPerm) :
66  VectorCoefficient(3), a(&A), origin(origin_), normal(normal_),
67  meshOrig(mesh), r2o(r2o_), perm(refPerm)
68  { }
69 
70  virtual void Eval(Vector &V, ElementTransformation &T,
71  const IntegrationPoint &ip);
72 
73  using VectorCoefficient::Eval;
74 };
75 
77  const IntegrationPoint &ip)
78 {
79  double x[3];
80  Vector transip(x, 3);
81 
82  T.Transform(ip, transip);
83 
84  const int elem = T.ElementNo;
85  const bool reflected = (*r2o)[elem] < 0;
86  const int originalElem = reflected ? -1 - (*r2o)[elem] : (*r2o)[elem];
87 
88  ElementTransformation *T_orig = meshOrig->GetElementTransformation(
89  originalElem);
90 
91  Vector rp(transip);
92 
93  if (reflected)
94  {
95  ReflectPoint(rp, origin, normal);
96  }
97 
98  IntegrationPoint ipo;
99 
100  a->Eval(V, *T_orig, ip);
101 
102  if (reflected)
103  {
104  // Map from ip in reflected elements to ip in initial elements, in mesh
105  // `reflected`. y = Ax + b in reference spaces, where A has 9 entries,
106  // and b has 3, totaling 12 unknowns. The number of data points is
107  // 3 * 8 = 24, so it is overdetermined. We use 4 points out of 8, (0,0,0),
108  // (1,0,0), (0,1,0), (0,0,1). For x=(0,0,0), b = y, and the other choices
109  // give the columns of A.
110 
111  // Permutation p is such that hex_reflected[i] = hex_init[p[i]]
112  const std::vector<int>& p = (*perm)[elem];
113 
114  // ip is on reflected hex. We map from the reflected hex to the initial
115  // hex, in reference space. Thus we use y = Ax + b, where x is in the
116  // reflected reference space, and y is in the initial hex reference space.
117 
118  const IntegrationRule *ir = Geometries.GetVertices(Geometry::CUBE);
119 
120  Vector b(3);
121  b[0] = (*ir)[p[0]].x;
122  b[1] = (*ir)[p[0]].y;
123  b[2] = (*ir)[p[0]].z;
124 
125  DenseMatrix A(3);
126 
127  // Vertex 1 is x=(1,0,0), so Ax is the first column of A.
128  A(0,0) = (*ir)[p[1]].x - b[0];
129  A(1,0) = (*ir)[p[1]].y - b[1];
130  A(2,0) = (*ir)[p[1]].z - b[2];
131 
132  // Vertex 3 is x=(0,1,0), so Ax is the second column of A.
133  A(0,1) = (*ir)[p[3]].x - b[0];
134  A(1,1) = (*ir)[p[3]].y - b[1];
135  A(2,1) = (*ir)[p[3]].z - b[2];
136 
137  // Vertex 4 is x=(0,0,1), so Ax is the third column of A.
138  A(0,2) = (*ir)[p[4]].x - b[0];
139  A(1,2) = (*ir)[p[4]].y - b[1];
140  A(2,2) = (*ir)[p[4]].z - b[2];
141 
142  Vector r(3);
143  Vector y(3);
144 
145  r[0] = ip.x;
146  r[1] = ip.y;
147  r[2] = ip.z;
148 
149  A.Mult(r, y);
150  y += b;
151 
152  ipo.x = y[0];
153  ipo.y = y[1];
154  ipo.z = y[2];
155 
156  a->Eval(V, *T_orig, ipo);
157  }
158 
159  if (reflected)
160  {
161  ReflectPoint(V, origin, normal);
162  }
163 }
164 
165 // Find perm such that h1[i] = h2[perm[i]]
166 void GetHexPermutation(Array<int> const& h1, Array<int> const& h2,
167  std::vector<int> & perm)
168 {
169  std::map<int, int> h2inv;
170  const int n = perm.size();
171 
172  for (int i=0; i<n; ++i)
173  {
174  h2inv[h2[i]] = i;
175  }
176 
177  for (int i=0; i<n; ++i)
178  {
179  perm[i] = h2inv[h1[i]];
180  }
181 }
182 
183 // This class facilitates constructing a hexahedral mesh one element at a time,
184 // using AddElement. The hexahedron input to AddElement is specified by vertices
185 // without requiring consistent global orientations. The mesh can be constructed
186 // simply by calling AddVertex for all vertices and then AddElement for all
187 // hexahedra. The mesh is not owned by this class and should be deleted outside
188 // this class.
189 class HexMeshBuilder
190 {
191 public:
192  HexMeshBuilder(int nv, int ne) : v2f(nv)
193  {
194  mesh = new Mesh(3, nv, ne);
195  f2v =
196  {
197  { {0,3,2,1},
198  {4,5,6,7},
199  {0,1,5,4},
200  {3,7,6,2},
201  {0,4,7,3},
202  {1,2,6,5}
203  }
204  };
205 
206  e2v =
207  {
208  { {0,1},
209  {1,2},
210  {2,3},
211  {0,3},
212  {4,5},
213  {5,6},
214  {6,7},
215  {4,7},
216  {0,4},
217  {1,5},
218  {2,6},
219  {3,7}
220  }
221  };
222  }
223 
224  /** @brief Add a single vertex to the mesh, specified by 3 coordinates. */
225  int AddVertex(const double *coords) const
226  {
227  return mesh->AddVertex(coords);
228  }
229 
230  /** @brief Add a single hexahedral element to the mesh, specified by 8
231  vertices. */
232  /** If reorder is false, the ordering in @vertices is used, so it must be
233  known in advance to have consistent orientations. Otherwise, a new
234  ordering will be found to ensure consistent orientations in the mesh.
235  */
236  int AddElement(Array<int> const& vertices, const bool reorder);
237 
238  Mesh *mesh;
239  std::vector<std::vector<int>> refPerm;
240 
241 private:
242  std::vector<std::vector<int>> faces;
243  std::vector<std::set<int>> v2f;
244  std::vector<std::vector<int>> f2e;
245  std::array<std::array<int, 4>, 6> f2v;
246  std::array<std::array<int, 2>, 12> e2v;
247 
248  int FindFourthVertexOnFace(Array<int> const& hex,
249  std::vector<int> const& v3) const;
250 
251  void ReorderHex(Array<int> & hex) const;
252  void ReverseHexX(Array<int> & hex) const;
253  void ReverseHexY(Array<int> & hex) const;
254  void ReverseHexZ(Array<int> & hex) const;
255  void ReverseHexFace(Array<int> & hex, const int face) const;
256  int FindHexFace(Array<int> const& hex, std::vector<int> const& face) const;
257 
258  bool ReorderHex_faceOrientations(Array<int> & hex) const;
259  void SaveHexFaces(const int elem, Array<int> const& hex);
260 };
261 
262 int HexMeshBuilder::AddElement(Array<int> const& vertices, const bool reorder)
263 {
264  MFEM_ASSERT(vertices.Size() == 8, "Hexahedron must have 8 vertices");
265 
266  Array<int> rvert(vertices);
267  if (reorder)
268  {
269  ReorderHex(rvert); // First reorder to set (0,0,0) and (1,1,1) vertices.
270 
271  // Now reorder to get consistent face orientations.
272  bool reordered = true;
273  int iter = 0;
274  do
275  {
276  reordered = ReorderHex_faceOrientations(rvert);
277  iter++;
278  MFEM_VERIFY(iter < 5, "");
279  }
280  while (reordered);
281 
282  std::vector<int> perm_e(8);
283  GetHexPermutation(rvert, vertices, perm_e);
284  refPerm.push_back(perm_e);
285  }
286  else
287  {
288  refPerm.push_back(std::vector<int> {0, 1, 2, 3, 4, 5, 6, 7});
289  }
290 
291  SaveHexFaces(mesh->GetNE(), rvert);
292 
293  Element * nel = mesh->NewElement(Geometry::Type::CUBE);
294 
295  nel->SetVertices(rvert);
296  return mesh->AddElement(nel);
297 }
298 
299 int HexMeshBuilder::FindFourthVertexOnFace(Array<int> const& hex,
300  std::vector<int> const& v3) const
301 {
302  int f0 = -1;
303  for (int f=0; f<6; ++f)
304  {
305  bool all3found = true;
306 
307  for (int i=0; i<3; ++i)
308  {
309  // Check whether v3[i] is in face f
310  bool found = false;
311 
312  for (int j=0; j<4; ++j)
313  {
314  if (hex[f2v[f][j]] == v3[i])
315  {
316  found = true;
317  }
318  }
319 
320  if (!found)
321  {
322  all3found = false;
323  break;
324  }
325  }
326 
327  if (all3found)
328  {
329  MFEM_ASSERT(f0 == -1, "");
330  f0 = f;
331  }
332  }
333 
334  MFEM_VERIFY(f0 >= 0, "");
335 
336  // Find the vertex of f0 not in v3
337  int v = -1;
338 
339  for (int j=0; j<4; ++j)
340  {
341  bool found = false;
342  for (int i=0; i<3; ++i)
343  {
344  // Check whether v3[i] is in face f
345  if (hex[f2v[f0][j]] == v3[i])
346  {
347  found = true;
348  }
349  }
350 
351  if (!found)
352  {
353  MFEM_ASSERT(v == -1, "");
354  v = hex[f2v[f0][j]];
355  }
356  }
357 
358  MFEM_VERIFY(v >= 0, "");
359 
360  return v;
361 }
362 
363 void HexMeshBuilder::ReorderHex(Array<int> & hex) const
364 {
365  MFEM_VERIFY(hex.Size() == 8, "hex");
366 
367  Array<int> h(hex);
368 
369  const int v0 = hex.Min();
370 
371  std::map<int, int> v2hex0;
372 
373  for (int i=0; i<hex.Size(); ++i)
374  {
375  v2hex0[hex[i]] = i;
376  }
377 
378  // Find the 3 vertices sharing an edge with v0.
379  std::vector<int> v0e;
380  for (int e=0; e<12; ++e)
381  {
382  if (v0 == hex[e2v[e][0]] || v0 == hex[e2v[e][1]])
383  {
384  v0e.push_back(e);
385  }
386  }
387 
388  MFEM_VERIFY(v0e.size() == 3, "");
389 
390  std::vector<int> v0n; // Neighbors of v0
391  for (auto e : v0e)
392  {
393  if (v0 == hex[e2v[e][0]])
394  {
395  v0n.push_back(hex[e2v[e][1]]);
396  }
397  else
398  {
399  v0n.push_back(hex[e2v[e][0]]);
400  }
401  }
402 
403  MFEM_VERIFY(v0n.size() == 3, "");
404 
405  sort(v0n.begin(), v0n.end());
406 
407  h[0] = v0;
408  h[1] = v0n[0];
409  h[3] = v0n[1];
410  h[4] = v0n[2];
411 
412  // Set h[2] by finding the face containing h[0], h[1], h[3]
413  std::vector<int> v3(3);
414  v3[0] = h[0];
415  v3[1] = h[1];
416  v3[2] = h[3];
417  h[2] = FindFourthVertexOnFace(hex, v3);
418 
419  // Set h[5] based on h[0], h[1], h[4]
420  v3[2] = h[4];
421  h[5] = FindFourthVertexOnFace(hex, v3);
422 
423  // Set h[7] based on h[0], h[3], h[4]
424  v3[1] = h[3];
425  h[7] = FindFourthVertexOnFace(hex, v3);
426 
427  // Set h[6] based on h[1], h[2], h[5]
428  v3[0] = h[1];
429  v3[1] = h[2];
430  v3[2] = h[5];
431  h[6] = FindFourthVertexOnFace(hex, v3);
432 
433  hex = h;
434 }
435 
436 void HexMeshBuilder::ReverseHexZ(Array<int> & hex) const
437 {
438  // faces {0,1,2,3} and {4,5,6,7} are reversed to become
439  // {0,3,2,1} and {4,7,6,5}
440  // This is accomplished by swapping vertices 1 and 3, and vertices 5 and 7.
441  int s = hex[1];
442  hex[1] = hex[3];
443  hex[3] = s;
444 
445  s = hex[5];
446  hex[5] = hex[7];
447  hex[7] = s;
448 }
449 
450 void HexMeshBuilder::ReverseHexY(Array<int> & hex) const
451 {
452  // faces {0,1,5,4} and {3,2,6,7} are reversed to become
453  // {0,4,5,1} and {3,7,6,2}
454  // This is accomplished by swapping vertices 1 and 4, and vertices 2 and 7.
455  int s = hex[1];
456  hex[1] = hex[4];
457  hex[4] = s;
458 
459  s = hex[2];
460  hex[2] = hex[7];
461  hex[7] = s;
462 }
463 
464 void HexMeshBuilder::ReverseHexX(Array<int> & hex) const
465 {
466  // faces {0,3,7,4} and {1,2,6,5} are reversed to become
467  // {0,4,7,3} and {1,5,6,2}
468  // This is accomplished by swapping vertices 3 and 4, and vertices 2 and 5.
469  int s = hex[4];
470  hex[4] = hex[3];
471  hex[3] = s;
472 
473  s = hex[5];
474  hex[5] = hex[2];
475  hex[2] = s;
476 }
477 
478 // Reverse face orientations without changing reference vertices 0 or 6.
479 void HexMeshBuilder::ReverseHexFace(Array<int> & hex, const int face) const
480 {
481  const int f = 2 * (face / 2); // f is in {0, 2, 4}
482 
483  switch (f)
484  {
485  case 0:
486  ReverseHexZ(hex);
487  break;
488  case 2:
489  ReverseHexY(hex);
490  break;
491  default: // case 4
492  ReverseHexX(hex);
493  }
494 }
495 
496 int HexMeshBuilder::FindHexFace(Array<int> const& hex,
497  std::vector<int> const& face) const
498 {
499  int localFace = -1;
500  for (int f=0; f<6; ++f)
501  {
502  std::vector<int> fv(4);
503  for (int i=0; i<4; ++i)
504  {
505  fv[i] = hex[f2v[f][i]];
506  }
507 
508  sort(fv.begin(), fv.end());
509 
510  if (fv == face)
511  {
512  MFEM_VERIFY(localFace == -1, "");
513  localFace = f;
514  }
515  }
516 
517  MFEM_VERIFY(localFace >= 0, "");
518 
519  return localFace;
520 }
521 
522 bool HexMeshBuilder::ReorderHex_faceOrientations(Array<int> & hex) const
523 {
524  std::vector<int> localFacesFound, globalFacesFound;
525  for (int f=0; f<6; ++f)
526  {
527  std::vector<int> fv(4);
528  for (int i=0; i<4; ++i)
529  {
530  fv[i] = hex[f2v[f][i]];
531  }
532 
533  sort(fv.begin(), fv.end());
534 
535  const int vmin = fv[0];
536  int globalFace = -1;
537  for (auto gf : v2f[vmin])
538  {
539  if (fv == faces[gf])
540  {
541  globalFace = gf;
542  }
543  }
544 
545  if (globalFace >= 0)
546  {
547  globalFacesFound.push_back(globalFace);
548  localFacesFound.push_back(f);
549  }
550  }
551 
552  const int numFoundFaces = globalFacesFound.size();
553 
554  for (int ff=0; ff<numFoundFaces; ++ff)
555  {
556  const int globalFace = globalFacesFound[ff];
557  const int localFace = localFacesFound[ff];
558 
559  MFEM_VERIFY(f2e[globalFace].size() == 1, "");
560  const int neighborElem = f2e[globalFace][0];
561 
562  Array<int> neighborElemVert;
563  mesh->GetElementVertices(neighborElem, neighborElemVert);
564  const int neighborLocalFace = FindHexFace(neighborElemVert, faces[globalFace]);
565 
566  std::vector<int> fv(4);
567  std::vector<int> nv(4);
568  for (int i=0; i<4; ++i)
569  {
570  fv[i] = hex[f2v[localFace][i]];
571  nv[i] = neighborElemVert[f2v[neighborLocalFace][i]];
572  }
573 
574  // As in Mesh::GetQuadOrientation, check whether fv and nv are oriented
575  // in the same direction.
576 
577  int id0;
578  for (id0 = 0; id0 < 4; id0++)
579  {
580  if (fv[id0] == nv[0])
581  {
582  break;
583  }
584  }
585 
586  MFEM_VERIFY(id0 < 4, "");
587 
588  bool same = (fv[(id0+1) % 4] == nv[1]);
589  if (same)
590  {
591  // Orientation should not be the same, so reverse the orientation of
592  // face localFace and its opposite face in hex.
593  ReverseHexFace(hex, localFace);
594  return true;
595  }
596  }
597 
598  return false;
599 }
600 
601 void HexMeshBuilder::SaveHexFaces(const int elem, Array<int> const& hex)
602 {
603  for (int f=0; f<6; ++f)
604  {
605  std::vector<int> fv(4);
606  for (int i=0; i<4; ++i)
607  {
608  fv[i] = hex[f2v[f][i]];
609  }
610 
611  sort(fv.begin(), fv.end());
612 
613  const int vmin = fv[0];
614  int globalFace = -1;
615  for (auto gf : v2f[vmin])
616  {
617  if (fv == faces[gf])
618  {
619  globalFace = gf;
620  }
621  }
622 
623  if (globalFace == -1)
624  {
625  // Face not found, so add it.
626  faces.push_back(fv);
627  globalFace = faces.size() - 1;
628 
629  std::vector<int> firstElem = {elem};
630  f2e.push_back(firstElem);
631  }
632  else
633  {
634  // Face found, so add elem to f2e
635  MFEM_VERIFY(f2e[globalFace].size() == 1 &&
636  f2e[globalFace][0] != elem, "");
637  f2e[globalFace].push_back(elem);
638  }
639 
640  MFEM_VERIFY(faces.size() == f2e.size(), "");
641  v2f[vmin].insert(globalFace);
642  }
643 }
644 
645 double GetElementEdgeMin(Mesh const& mesh, const int elem)
646 {
647  Array<int> edges, cor;
648  mesh.GetElementEdges(elem, edges, cor);
649 
650  double diam = -1.0;
651  for (auto e : edges)
652  {
653  Array<int> vert;
654  mesh.GetEdgeVertices(e, vert);
655  const double *v0 = mesh.GetVertex(vert[0]);
656  const double *v1 = mesh.GetVertex(vert[1]);
657 
658  double L = 0.0;
659  for (int i=0; i<3; ++i)
660  {
661  L += (v0[i] - v1[i]) * (v0[i] - v1[i]);
662  }
663 
664  L = sqrt(L);
665 
666  if (diam < 0.0 || L < diam) { diam = L; }
667  }
668 
669  return diam;
670 }
671 
672 void FindElementsTouchingPlane(Mesh const& mesh, Vector const& origin,
673  Vector const& normal, std::vector<int> & el)
674 {
675  const double relTol = 1.0e-6;
676  Vector diff(3);
677 
678  for (int e=0; e<mesh.GetNE(); ++e)
679  {
680  const double diam = GetElementEdgeMin(mesh, e);
681  Array<int> vert;
682  mesh.GetElementVertices(e, vert);
683 
684  bool onplane = false;
685  for (auto v : vert)
686  {
687  const double *vcrd = mesh.GetVertex(v);
688  for (int i=0; i<3; ++i)
689  {
690  diff[i] = vcrd[i] - origin[i];
691  }
692 
693  if (std::abs(diff * normal) < relTol * diam)
694  {
695  onplane = true;
696  }
697  }
698 
699  if (onplane) { el.push_back(e); }
700  }
701 }
702 
703 // Order the elements in layers, starting at the plane of reflection.
704 bool GetMeshElementOrder(Mesh const& mesh, Vector const& origin,
705  Vector const& normal, std::vector<int> & elOrder)
706 {
707  const int ne = mesh.GetNE();
708  elOrder.assign(ne, -1);
709 
710  std::vector<bool> elementMarked;
711 
712  elementMarked.assign(ne, false);
713 
714  std::vector<int> layer;
715  FindElementsTouchingPlane(mesh, origin, normal, layer);
716 
717  if (layer.size() == 0)
718  {
719  // If the mesh does not touch the plane, any ordering will work.
720  for (int i=0; i<ne; ++i)
721  {
722  elOrder[i] = i;
723  }
724 
725  return false;
726  }
727 
728  int cnt = 0;
729  while (cnt < ne)
730  {
731  for (auto e : layer)
732  {
733  elOrder[cnt] = e;
734  cnt++;
735  elementMarked[e] = true;
736  }
737 
738  if (cnt == ne) { break; }
739 
740  std::set<int> layerNext;
741  for (auto e : layer)
742  {
743  Array<int> nghb = mesh.FindFaceNeighbors(e);
744  for (auto n : nghb)
745  {
746  if (!elementMarked[n]) { layerNext.insert(n); }
747  }
748  }
749 
750  MFEM_VERIFY(layerNext.size() > 0, "");
751 
752  layer.clear();
753  layer.reserve(layerNext.size());
754  for (auto e : layerNext)
755  {
756  layer.push_back(e);
757  }
758  }
759 
760  MFEM_VERIFY(cnt == ne, "");
761 
762  return true;
763 }
764 
765 Mesh* ReflectHighOrderMesh(Mesh & mesh, Vector origin, Vector normal)
766 {
767  MFEM_VERIFY(mesh.Dimension() == 3, "Only 3D meshes can be reflected");
768 
769  // Find the minimum edge length, to use for a relative tolerance.
770  double minLength = 0.0;
771  for (int i=0; i<mesh.GetNE(); i++)
772  {
773  Array<int> vert;
774  mesh.GetEdgeVertices(i, vert);
775  const Vector v0(mesh.GetVertex(vert[0]), 3);
776  const Vector v1(mesh.GetVertex(vert[1]), 3);
777  Vector diff(3);
778  subtract(v0, v1, diff);
779  const double length = diff.Norml2();
780  if (i == 0 || length < minLength)
781  {
782  minLength = length;
783  }
784  }
785 
786  const double relTol = 1.0e-6;
787 
788  // Find vertices in reflection plane.
789  std::set<int> planeVertices;
790  for (int i=0; i<mesh.GetNV(); i++)
791  {
792  Vector v(mesh.GetVertex(i), 3);
793  Vector diff(3);
794  subtract(v, origin, diff);
795  const double ip = diff * normal;
796  if (fabs(ip) < relTol * minLength)
797  {
798  planeVertices.insert(i);
799  }
800  }
801 
802  const int nv = mesh.GetNV();
803  const int ne = mesh.GetNE();
804 
805  std::vector<int> r2o;
806 
807  const int nv_reflected = (2*nv) - planeVertices.size();
808 
809  HexMeshBuilder builder(nv_reflected, 2*ne);
810 
811  r2o.assign(2*ne, -2-ne); // Initialize to invalid value.
812 
813  std::vector<int> v2r;
814  v2r.assign(mesh.GetNV(), -1);
815 
816  // Copy vertices
817  for (int v=0; v<mesh.GetNV(); v++)
818  {
819  builder.AddVertex(mesh.GetVertex(v));
820  }
821 
822  for (int v=0; v<mesh.GetNV(); v++)
823  {
824  // Check whether vertex v is in the plane
825  if (planeVertices.find(v) == planeVertices.end())
826  {
827  // For vertices not in plane, reflect and add.
828  Vector vr(3);
829  for (int i=0; i<3; ++i)
830  {
831  vr[i] = mesh.GetVertex(v)[i];
832  }
833 
834  ReflectPoint(vr, origin, normal);
835 
836  v2r[v] = builder.AddVertex(vr.GetData());
837  }
838  }
839 
840  std::vector<int> elOrder;
841  const bool onPlane = GetMeshElementOrder(mesh, origin, normal, elOrder);
842 
843  for (int eidx=0; eidx<mesh.GetNE(); eidx++)
844  {
845  const int e = elOrder[eidx];
846 
847  // Copy the original element
848  Array<int> elvert;
849  mesh.GetElementVertices(e, elvert);
850 
851  MFEM_VERIFY(elvert.Size() == 8, "Only hexahedral elements are supported");
852 
853  const int copiedElem = builder.AddElement(elvert, false);
854  r2o[copiedElem] = e;
855 
856  // Add the new reflected element
857  Array<int> rvert(elvert.Size());
858  for (int i=0; i<elvert.Size(); ++i)
859  {
860  const int v = elvert[i];
861  rvert[i] = (v2r[v] == -1) ? v : v2r[v];
862  }
863 
864  const int newElem = builder.AddElement(rvert, onPlane);
865  r2o[newElem] = -1 - e;
866  }
867 
868  Mesh *reflected = builder.mesh;
869 
870  // Set attributes
871  MFEM_VERIFY((int) r2o.size() == reflected->GetNE(), "");
872  for (int i = 0; i < (int) r2o.size(); ++i)
873  {
874  const int e = (r2o[i] >= 0) ? r2o[i] : -1 - r2o[i];
875  reflected->SetAttribute(i, mesh.GetAttribute(e));
876  }
877 
878  // In order to set boundary attributes, first set a map from original mesh
879  // boundary elements to reflected mesh boundary elements, by using the vertex
880  // map v2r. Note that for v < mesh.GetNV(), vertex v of `mesh` coincides with
881  // vertex v of `reflected`, and if that vertex is not in the reflection
882  // plane, v2r[v] >= mesh.GetNV() is the index of the vertex in `reflected`
883  // that is its reflection.
884 
885  // Identify each quadrilateral boundary element with the unique pair of
886  // vertices (v1, v2) such that v1 is the minimum vertex index in the
887  // quadrilateral, and v2 is diagonally opposite v1.
888 
889  std::map<std::pair<int, int>, int> mapBE;
890  for (int i=0; i<mesh.GetNBE(); ++i)
891  {
892  const Element *be = mesh.GetBdrElement(i);
893  Array<int> v;
894  be->GetVertices(v);
895  MFEM_VERIFY(v.Size() == 4, "Boundary elements must be quadrilateral");
896 
897  const int v1 = v.Min();
898  int v1i = -1;
899  for (int j=0; j<v.Size(); ++j)
900  {
901  if (v[j] == v1)
902  {
903  v1i = j;
904  }
905  }
906 
907  const int v2 = v[(v1i + 2) % 4];
908 
909  mapBE[std::pair<int, int>(v1, v2)] = i;
910 
911  // Find the indices of vertices in `reflected` of the reflected quadrilateral.
912  Array<int> rv(4);
913  int rv1 = -1; // Find the minimum reflected vertex index.
914  int rv1i = -1;
915  bool inPlane = true;
916  for (int j=0; j<v.Size(); ++j)
917  {
918  rv[j] = (v2r[v[j]] == -1) ? v[j] : v2r[v[j]];
919 
920  if (v2r[v[j]] != -1)
921  {
922  inPlane = false;
923  }
924 
925  if (rv1 == -1 || rv[j] < rv1)
926  {
927  rv1 = rv[j];
928  rv1i = j;
929  }
930  }
931 
932  // Note that in-plane boundary elements are skipped.
933  if (!inPlane)
934  {
935  const int rv2 = rv[(rv1i + 2) % 4];
936 
937  mapBE[std::pair<int, int>(rv1, rv2)] = i;
938 
939  mfem::Swap(rv[0], rv[2]); // Fix the orientation
940 
941  const Geometry::Type orig_geom = mesh.GetBdrElementBaseGeometry(i);
942  Element *rbe = reflected->NewElement(orig_geom);
943  rbe->SetVertices(v);
944  reflected->AddBdrElement(rbe);
945 
946  rbe = reflected->NewElement(orig_geom);
947  rbe->SetVertices(rv);
948  reflected->AddBdrElement(rbe);
949  }
950  }
951 
952  for (int i=0; i<reflected->GetNBE(); ++i)
953  {
954  Element *be = reflected->GetBdrElement(i);
955  Array<int> rv;
956  be->GetVertices(rv);
957  MFEM_VERIFY(rv.Size() == 4, "Boundary elements must be quadrilateral");
958 
959  // Reflected boundary element i is identified with vertices
960 
961  const int v1 = rv.Min();
962  int v1i = -1;
963  for (int j=0; j<rv.Size(); ++j)
964  {
965  if (rv[j] == v1)
966  {
967  v1i = j;
968  }
969  }
970 
971  const int v2 = rv[(v1i + 2) % 4];
972 
973  const int originalBE = mapBE[std::pair<int, int>(v1, v2)];
974  const int originalAttribute = mesh.GetBdrAttribute(originalBE);
975  reflected->SetBdrAttribute(i, originalAttribute);
976  }
977 
978  reflected->FinalizeTopology();
979  reflected->Finalize();
980  reflected->RemoveUnusedVertices();
981 
982  if (mesh.GetNodes())
983  {
984  // Extract Nodes GridFunction and determine its type
985  const GridFunction * Nodes = mesh.GetNodes();
986  const FiniteElementSpace * fes = Nodes->FESpace();
987 
988  Ordering::Type ordering = fes->GetOrdering();
989  int order = fes->FEColl()->GetOrder();
990  int sdim = mesh.SpaceDimension();
991  bool discont =
992  dynamic_cast<const L2_FECollection*>(fes->FEColl()) != NULL;
993 
994  // Set curvature of the same type as original mesh
995  reflected->SetCurvature(order, discont, sdim, ordering);
996 
997  GridFunction * reflected_nodes = reflected->GetNodes();
998  GridFunction newReflectedNodes(*reflected_nodes);
999 
1000  VectorGridFunctionCoefficient nodesCoef(Nodes);
1001 
1002  ReflectedCoefficient rc(nodesCoef, origin, normal, &r2o, &mesh,
1003  &builder.refPerm);
1004 
1005  newReflectedNodes.ProjectCoefficient(rc);
1006  *reflected_nodes = newReflectedNodes;
1007  }
1008 
1009  return reflected;
1010 }
1011 
1012 int main(int argc, char *argv[])
1013 {
1014  // Parse command-line options.
1015  const char *mesh_file = "../../data/pipe-nurbs.mesh";
1016  bool visualization = 1;
1017  Vector normal(3);
1018  Vector origin(3);
1019 
1020  normal = 0.0;
1021  normal[2] = 1.0;
1022  origin = 0.0;
1023 
1024  OptionsParser args(argc, argv);
1025  args.AddOption(&mesh_file, "-m", "--mesh",
1026  "Mesh file to use.");
1027  args.AddOption(&normal, "-n", "--normal",
1028  "Normal vector of plane.");
1029  args.AddOption(&origin, "-o", "--origin",
1030  "A point in the plane.");
1031  args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
1032  "--no-visualization",
1033  "Enable or disable GLVis visualization.");
1034  args.Parse();
1035  if (!args.Good())
1036  {
1037  args.PrintUsage(cout);
1038  return 1;
1039  }
1040  args.PrintOptions(cout);
1041 
1042  MFEM_VERIFY(fabs(normal.Norml2() - 1.0) < 1.0e-14, "");
1043 
1044  Mesh mesh(mesh_file, 0, 0);
1045 
1046  Mesh *reflected = ReflectHighOrderMesh(mesh, origin, normal);
1047 
1048  // Save the final mesh
1049  ofstream mesh_ofs("reflected.mesh");
1050  mesh_ofs.precision(8);
1051  reflected->Print(mesh_ofs);
1052 
1053  if (visualization)
1054  {
1055  // GLVis server to visualize to
1056  char vishost[] = "localhost";
1057  int visport = 19916;
1058  socketstream sol_sock(vishost, visport);
1059  sol_sock.precision(8);
1060  sol_sock << "mesh\n" << *reflected << flush;
1061  }
1062 
1063  delete reflected;
1064 
1065  return 0;
1066 }
T Min() const
Find the minimal element in the array, using the comparison operator < for class T.
Definition: array.cpp:85
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void GetElementEdges(int i, Array< int > &edges, Array< int > &cor) const
Return the indices and the orientations of all edges of element i.
Definition: mesh.cpp:5945
Base class for vector Coefficients that optionally depend on time and space.
virtual void GetVertices(Array< int > &v) const =0
Returns element&#39;s vertices.
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:324
int Dimension() const
Definition: mesh.hpp:1047
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:454
int size
Size of the array.
Definition: array.hpp:51
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:218
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void SetVertices(const int *ind)
Set the indices the element according to the input.
Definition: element.cpp:17
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:150
STL namespace.
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:265
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:939
Geometry Geometries
Definition: fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:314
Element * NewElement(int geom)
Definition: mesh.cpp:3854
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:933
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1544
bool GetMeshElementOrder(Mesh const &mesh, Vector const &origin, Vector const &normal, std::vector< int > &elOrder)
Definition: reflector.cpp:704
void FindElementsTouchingPlane(Mesh const &mesh, Vector const &origin, Vector const &normal, std::vector< int > &el)
Definition: reflector.cpp:672
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
constexpr char vishost[]
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
double b
Definition: lissajous.cpp:42
constexpr int visport
Type
Ordering methods:
Definition: fespace.hpp:33
virtual void SetCurvature(int order, bool discont=false, int space_dim=-1, int ordering=1)
Set the curvature of the mesh nodes using the given polynomial degree.
Definition: mesh.cpp:5356
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:6029
virtual void Eval(DenseMatrix &M, ElementTransformation &T, const IntegrationRule &ir)
Evaluate the vector coefficient in the element described by T at all points of ir, storing the result in M.
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:1053
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1148
Array< int > FindFaceNeighbors(const int elem) const
Returns the sorted, unique indices of elements sharing a face with element elem, including elem...
Definition: mesh.cpp:6184
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
void GetHexPermutation(Array< int > const &h1, Array< int > const &h2, std::vector< int > &perm)
Definition: reflector.cpp:166
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1550
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:635
void ReflectPoint(Vector &p, Vector const &origin, Vector const &normal)
Definition: reflector.cpp:37
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:11742
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:461
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set &#39;var&#39; to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:2902
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
Definition: mesh.hpp:1553
int SpaceDimension() const
Definition: mesh.hpp:1048
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
double a
Definition: lissajous.cpp:41
Class for integration point with weight.
Definition: intrules.hpp:25
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:2995
double GetElementEdgeMin(Mesh const &mesh, const int elem)
Definition: reflector.cpp:645
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2396
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:804
int AddBdrElement(Element *elem)
Definition: mesh.cpp:1828
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1547
Mesh * ReflectHighOrderMesh(Mesh &mesh, Vector origin, Vector normal)
Definition: reflector.cpp:765
Vector data type.
Definition: vector.hpp:60
virtual void Print(std::ostream &os=mfem::out) const
Definition: mesh.hpp:1749
Vector coefficient defined by a vector GridFunction.
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:7949
RefCoord s[3]
const Element * GetBdrElement(int i) const
Definition: mesh.hpp:1085
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:1111
Abstract data type element.
Definition: element.hpp:28
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320
int main(int argc, char *argv[])
Definition: reflector.cpp:1012
double f(const Vector &p)