MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
reflector.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
33using namespace std;
34using namespace mfem;
35
36
37void 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 real_t ip = diff * normal;
43
44 diff = normal;
45 diff *= -2.0 * ip;
46
47 add(p, diff, p);
48}
49
50class ReflectedCoefficient : public VectorCoefficient
51{
52private:
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
62public:
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
74};
75
77 const IntegrationPoint &ip)
78{
79 real_t 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
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
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]]
166void 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.
189class HexMeshBuilder
190{
191public:
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 real_t *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
241private:
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
262int 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
299int 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
363void 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
436void 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
450void 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
464void 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.
479void 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
496int 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
522bool 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
601void 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
645real_t GetElementEdgeMin(Mesh const& mesh, const int elem)
646{
647 Array<int> edges, cor;
648 mesh.GetElementEdges(elem, edges, cor);
649
650 real_t diam = -1.0;
651 for (auto e : edges)
652 {
653 Array<int> vert;
654 mesh.GetEdgeVertices(e, vert);
655 const real_t *v0 = mesh.GetVertex(vert[0]);
656 const real_t *v1 = mesh.GetVertex(vert[1]);
657
658 real_t 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
672void FindElementsTouchingPlane(Mesh const& mesh, Vector const& origin,
673 Vector const& normal, std::vector<int> & el)
674{
675 const real_t relTol = 1.0e-6;
676 Vector diff(3);
677
678 for (int e=0; e<mesh.GetNE(); ++e)
679 {
680 const real_t 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 real_t *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.
704bool 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
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 real_t 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 real_t length = diff.Norml2();
780 if (i == 0 || length < minLength)
781 {
782 minLength = length;
783 }
784 }
785
786 const real_t 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 real_t ip = diff * normal;
796 if (std::abs(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.GetBdrElementGeometry(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
1012int 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(std::abs(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
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Abstract data type element.
Definition: element.hpp:29
virtual void GetVertices(Array< int > &v) const =0
Get the indices defining the vertices.
virtual void SetVertices(const Array< int > &v)=0
Set the indices defining the vertices.
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:220
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:725
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:727
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
Definition: geom.cpp:293
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:31
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:696
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2346
Class for integration point with weight.
Definition: intrules.hpp:35
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:330
Mesh data type.
Definition: mesh.hpp:56
Element * NewElement(int geom)
Definition: mesh.cpp:4401
int AddBdrElement(Element *elem)
Definition: mesh.cpp:2028
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:7224
Geometry::Type GetBdrElementGeometry(int i) const
Definition: mesh.hpp:1376
int GetAttribute(int i) const
Return the attribute of element i.
Definition: mesh.hpp:1333
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1438
int GetBdrAttribute(int i) const
Return the attribute of boundary element i.
Definition: mesh.hpp:1339
void SetAttribute(int i, int attr)
Set the attribute of element i.
Definition: mesh.hpp:1336
void FinalizeTopology(bool generate_bdr=true)
Finalize the construction of the secondary topology (connectivity) data of a Mesh.
Definition: mesh.cpp:3135
virtual void Print(std::ostream &os=mfem::out, const std::string &comments="") const
Definition: mesh.hpp:2288
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
const Element * GetBdrElement(int i) const
Return pointer to the i'th boundary element object.
Definition: mesh.hpp:1298
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1163
void GetNodes(Vector &node_coord) const
Definition: mesh.cpp:8973
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:1223
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:7069
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:1229
virtual void Finalize(bool refine=false, bool fix_orientation=false)
Finalize the construction of a general Mesh.
Definition: mesh.cpp:3241
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:6985
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:6211
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
Definition: mesh.hpp:1342
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition: mesh.hpp:1265
void RemoveUnusedVertices()
Remove unused vertices and rebuild mesh connectivity.
Definition: mesh.cpp:12872
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
Definition: optparser.cpp:151
void PrintUsage(std::ostream &out) const
Print the usage message.
Definition: optparser.cpp:462
void PrintOptions(std::ostream &out) const
Print the options.
Definition: optparser.cpp:331
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 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition: optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Definition: optparser.hpp:159
Type
Ordering methods:
Definition: fespace.hpp:34
Base class for vector Coefficients that optionally depend on time and space.
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,...
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Vector coefficient defined by a vector GridFunction.
Vector data type.
Definition: vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:832
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:227
real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12, int max_its=10)
Bregman projection of ρ = sigmoid(ψ) onto the subspace ∫_Ω ρ dx = θ vol(Ω) as follows:
Definition: ex37.cpp:72
int visport
char vishost[]
int main()
real_t b
Definition: lissajous.cpp:42
real_t a
Definition: lissajous.cpp:41
real_t f(const Vector &p)
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:648
Geometry Geometries
Definition: fe.cpp:49
void add(const Vector &v1, const Vector &v2, Vector &v)
Definition: vector.cpp:316
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:472
float real_t
Definition: config.hpp:43
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53
RefCoord s[3]
void GetHexPermutation(Array< int > const &h1, Array< int > const &h2, std::vector< int > &perm)
Definition: reflector.cpp:166
real_t GetElementEdgeMin(Mesh const &mesh, const int elem)
Definition: reflector.cpp:645
void ReflectPoint(Vector &p, Vector const &origin, Vector const &normal)
Definition: reflector.cpp:37
bool GetMeshElementOrder(Mesh const &mesh, Vector const &origin, Vector const &normal, std::vector< int > &elOrder)
Definition: reflector.cpp:704
Mesh * ReflectHighOrderMesh(Mesh &mesh, Vector origin, Vector normal)
Definition: reflector.cpp:765
void FindElementsTouchingPlane(Mesh const &mesh, Vector const &origin, Vector const &normal, std::vector< int > &el)
Definition: reflector.cpp:672