MFEM v4.8.0
Finite element discretization library
Loading...
Searching...
No Matches
intrules_cut.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2025, 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// Implementation of Surface and Cutcell IntegrationRule(s) classes
13
14#include "fem.hpp"
15#include <cmath>
16
17using namespace std;
18
19namespace mfem
20{
21
23{
24 MFEM_VERIFY(order > 0, "Invalid input");
25 Order = order;
26}
27
29{
30 MFEM_VERIFY(order > 0, "Invalid input");
31 lsOrder = order;
32}
33
34#ifdef MFEM_USE_ALGOIM
36 &Tr,
37 IntegrationRule &result)
38{
39 GenerateLSVector(Tr,LvlSet);
40
41 const int dim=pe->GetDim();
42 int np1d=CutIntegrationRules::Order/2+1;
43 if (dim==2)
44 {
45 LevelSet2D ls(pe,lsvec);
46 auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
47 2, -1, np1d);
48 result.SetSize(q.nodes.size());
50 for (size_t i=0; i<q.nodes.size(); i++)
51 {
52 IntegrationPoint& ip=result.IntPoint(i);
53 ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
54 }
55 }
56 else
57 {
58 LevelSet3D ls(pe,lsvec);
59 auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
60 3, -1, np1d);
61
62 result.SetSize(q.nodes.size());
64 for (size_t i=0; i<q.nodes.size(); i++)
65 {
66 IntegrationPoint& ip=result.IntPoint(i);
67 ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
68 }
69 }
70
71}
72
74 IntegrationRule &result,
75 const IntegrationRule *sir)
76{
77 GenerateLSVector(Tr,LvlSet);
78
79 const int dim=pe->GetDim();
80 int np1d=CutIntegrationRules::Order/2+1;
81 if (dim==2)
82 {
83 LevelSet2D ls(pe,lsvec);
84 auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
85 -1, -1, np1d);
86 result.SetSize(q.nodes.size());
88 for (size_t i=0; i<q.nodes.size(); i++)
89 {
90 IntegrationPoint& ip=result.IntPoint(i);
91 ip.Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
92 }
93 }
94 else
95 {
96 LevelSet3D ls(pe,lsvec);
97 auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
98 -1, -1, np1d);
99
100 result.SetSize(q.nodes.size());
102 for (size_t i=0; i<q.nodes.size(); i++)
103 {
104 IntegrationPoint& ip=result.IntPoint(i);
105 ip.Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
106 }
107 }
108
109}
110
112 const IntegrationRule &sir,
113 Vector &weights)
114{
115 GenerateLSVector(Tr,LvlSet);
116
117 DenseMatrix bmat; // gradients of the shape functions in isoparametric space
118 DenseMatrix pmat; // gradients of the shape functions in physical space
119 Vector inormal; // normal to the level set in isoparametric space
120 Vector tnormal; // normal to the level set in physical space
121 bmat.SetSize(pe->GetDof(),pe->GetDim());
122 pmat.SetSize(pe->GetDof(),pe->GetDim());
123 inormal.SetSize(pe->GetDim());
124 tnormal.SetSize(pe->GetDim());
125
126 weights.SetSize(sir.GetNPoints());
127
128 for (int j = 0; j < sir.GetNPoints(); j++)
129 {
130 const IntegrationPoint &ip = sir.IntPoint(j);
131 Tr.SetIntPoint(&ip);
132 pe->CalcDShape(ip,bmat);
133 Mult(bmat, Tr.InverseJacobian(), pmat);
134 // compute the normal to the LS in isoparametric space
135 bmat.MultTranspose(lsvec,inormal);
136 // compute the normal to the LS in physical space
137 pmat.MultTranspose(lsvec,tnormal);
138 weights[j]= tnormal.Norml2() / inormal.Norml2();
139 }
140
141}
142
143void AlgoimIntegrationRules::GenerateLSVector(ElementTransformation &Tr,
145{
146 //check if the coefficient is already projected
147 if (currentElementNo==Tr.ElementNo)
148 {
149 if (currentLvlSet==lvlset)
150 {
151 if (currentGeometry==Tr.GetGeometryType())
152 {
153 return;
154 }
155 }
156 }
157
158 currentElementNo=Tr.ElementNo;
159
160 if (currentGeometry!=Tr.GetGeometryType())
161 {
162 delete le;
163 delete pe;
164 currentGeometry=Tr.GetGeometryType();
166 {
167 pe=new H1Pos_QuadrilateralElement(lsOrder);
168 le=new H1_QuadrilateralElement(lsOrder);
169 }
171 {
172 pe=new H1Pos_HexahedronElement(lsOrder);
173 le=new H1_HexahedronElement(lsOrder);
174 }
175 else
176 {
177 MFEM_ABORT("Currently MFEM + Algoim supports only quads and hexes.");
178 }
179
180 T.SetSize(pe->GetDof());
181 pe->Project(*le,Tr,T);
182 //The transformation matrix depends only on the geometry for change of basis
183 }
184
185 currentLvlSet=lvlset;
186 const IntegrationRule &ir=le->GetNodes();
187 lsvec.SetSize(ir.GetNPoints());
188 lsfun.SetSize(ir.GetNPoints());
189 for (int i=0; i<ir.GetNPoints(); i++)
190 {
191 const IntegrationPoint &ip = ir.IntPoint(i);
192 Tr.SetIntPoint(&ip);
193 lsfun(i)=lvlset->Eval(Tr,ip);
194 }
195 T.Mult(lsfun,lsvec);
196}
197
198#endif
199
200#ifdef MFEM_USE_LAPACK
201
203 int lsO, ElementTransformation& Tr)
204{
205 Init(order, levelset, lsO);
206 dim = Tr.GetDimension();
207 if (Tr.GetDimension() == 1)
208 {
209 nBasis = -1;
211 ir = irs.Get(Geometry::SEGMENT, 0);
212 }
213 else
214 {
215 if (Tr.GetDimension() == 2)
216 {
217 nBasis = 2 * (Order + 1) + static_cast<int>(Order * (Order + 1) / 2);
218 }
219 else if (Tr.GetDimension() == 3)
220 {
221 if (Order== 0) { nBasis = 3; }
222 else if (Order== 1) { nBasis = 11; }
223 else if (Order== 2) { nBasis = 26; }
224 else if (Order== 3) { nBasis = 50; }
225 else if (Order== 4) { nBasis = 85; }
226 else if (Order== 5) { nBasis = 133; }
227 else if (Order== 6) { nBasis = 196; }
228 else if (Order>= 7) { nBasis = 276; Order = 7; }
229 }
230
231 // compute the quadrature points
232 int qorder = 0;
234 ir = irs.Get(Tr.GetGeometryType(), qorder);
235 for (; ir.GetNPoints() <= nBasis; qorder++)
236 {
237 ir = irs.Get(Tr.GetGeometryType(), qorder);
238 }
239 }
240}
241
243 int lsO, ElementTransformation& Tr)
244{
245 order++;
246 InitSurface(order, levelset, lsO, Tr);
247 Order--;
248
249 nBasisVolume = 0;
250 if (Tr.GetDimension() == 1)
251 {
252 nBasisVolume = -1;
255 }
256 else
257 {
258 if (Tr.GetDimension() == 2)
259 {
260 nBasisVolume = (int)((Order + 1) * (Order + 2) / 2);
261 }
262 else if (Tr.GetDimension() == 3)
263 {
264 for (int p = 0; p <= Order; p++)
265 {
266 nBasisVolume +=(int)((p + 1) * (p + 2) / 2);
267 }
268 }
269
270 // assemble the matrix
272 for (int ip = 0; ip < ir.GetNPoints(); ip++)
273 {
274 Vector shape;
275 if (Tr.GetDimension() == 2)
276 {
277 Basis2D(ir.IntPoint(ip), shape);
278 }
279 else if (Tr.GetDimension() == 3)
280 {
281 Basis3D(ir.IntPoint(ip), shape);
282 }
283
284 Mat.SetCol(ip, shape);
285 }
286
287 // compute the SVD for the matrix
288 VolumeSVD = new DenseMatrixSVD(Mat, 'A', 'A');
289 VolumeSVD->Eval(Mat);
290 }
291
292}
293
295{
296 int elem = Tr.ElementNo;
297 const Mesh *mesh = Tr.mesh;
298
299 if (FaceIP.Size() == 0)
300 {
302 FaceWeightsComp = 0.;
303 }
304
305 const Element* me = mesh->GetElement(elem);
307 mesh->GetElementTransformation(elem, &Trafo);
308
309 Array<int> faces;
310 Array<int> cor;
311 mesh->GetElementFaces(elem, faces, cor);
312
313 for (int face = 0; face < me->GetNFaces(); face++)
314 {
315 if (FaceWeightsComp(faces[face]) == 0.)
316 {
317 FaceWeightsComp(faces[face]) = 1.;
318 Array<int> verts;
319 mesh->GetFaceVertices(faces[face], verts);
320 Vector pointA(mesh->SpaceDimension());
321 Vector pointB(mesh->SpaceDimension());
322 Vector pointC(mesh->SpaceDimension());
323 Vector pointD(mesh->SpaceDimension());
324 for (int d = 0; d < mesh->SpaceDimension(); d++)
325 {
326 pointA(d) = (mesh->GetVertex(verts[0]))[d];
327 pointB(d) = (mesh->GetVertex(verts[1]))[d];
328 pointC(d) = (mesh->GetVertex(verts[2]))[d];
329 pointD(d) = (mesh->GetVertex(verts[3]))[d];
330 }
331
332 // TODO - don't we lose the curvature with this local mesh setup?
333 Mesh local_mesh(2,4,1,0,3);
334 local_mesh.AddVertex(pointA);
335 local_mesh.AddVertex(pointB);
336 local_mesh.AddVertex(pointC);
337 local_mesh.AddVertex(pointD);
338 local_mesh.AddQuad(0,1,2,3);
339 local_mesh.FinalizeQuadMesh(1);
341 local_mesh.GetElementTransformation(0, &faceTrafo);
342
343 // The 3D face integrals are computed as 2D volumetric integrals.
344 // The 2D face integrals are computed as 1D volumetric integrals.
346 IntegrationRule FaceRule;
347 FaceRules.GetVolumeIntegrationRule(faceTrafo, FaceRule);
348 if (FaceIP.Size() != FaceRule.Size())
349 {
350 FaceIP.SetSize(FaceRule.Size());
351 for (int ip = 0; ip < FaceRule.GetNPoints(); ip++)
352 {
353 FaceIP[ip].index = ip;
354 IntegrationPoint &intp = FaceIP[ip];
355 intp.x = FaceRule.IntPoint(ip).x;
356 intp.y = FaceRule.IntPoint(ip).y;
357 intp.weight = 0.;
358 }
359
360 FaceWeights.SetSize(FaceRule.GetNPoints(), mesh->GetNFaces());
361 FaceWeights = 0.;
362
363 FaceWeightsComp = 0.;
364 FaceWeightsComp(faces[face]) = 1.;
365 }
366
367 for (int ip = 0; ip < FaceRule.GetNPoints(); ip++)
368 {
369 FaceWeights(ip, faces[face]) = FaceRule.IntPoint(ip).weight;
370 }
371 }
372 }
373
374 mesh->GetElementTransformation(elem, &Trafo);
375}
376
378{
379 IntegrationPoint& intp = ir.IntPoint(0);
380
382 ip0.x = 0.;
384 ip1.x = 1.;
385 Tr.SetIntPoint(&ip0);
386 if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip1) < 0.)
387 {
389 ip2.x = .5;
390 while (LvlSet->Eval(Tr, ip2) > tol_1
391 || LvlSet->Eval(Tr, ip2) < -tol_1)
392 {
393 if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip2) < 0.)
394 {
395 ip1.x = ip2.x;
396 }
397 else
398 {
399 ip0.x = ip2.x;
400 }
401
402 ip2.x = (ip1.x + ip0.x) / 2.;
403 }
404 intp.x = ip2.x;
405 intp.weight = 1. / Tr.Weight();
406 }
407 else if (LvlSet->Eval(Tr, ip0) > 0. && LvlSet->Eval(Tr, ip1) <= tol_1)
408 {
409 intp.x = 1.;
410 intp.weight = 1. / Tr.Weight();
411 }
412 else if (LvlSet->Eval(Tr, ip1) > 0. && LvlSet->Eval(Tr, ip0) <= tol_1)
413 {
414 intp.x = 0.;
415 intp.weight = 1. / Tr.Weight();
416 }
417 else
418 {
419 intp.x = .5;
420 intp.weight = 0.;
421 }
422}
423
425{
426 IntegrationPoint intp;
427
429 ip0.x = 0.;
431 ip1.x = 1.;
432 Tr.SetIntPoint(&ip0);
433 if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip1) < 0.)
434 {
436 ip2.x = .5;
437 while (LvlSet->Eval(Tr, ip2) > 1e-12
438 || LvlSet->Eval(Tr, ip2) < -1e-12)
439 {
440 if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip2) < 0.)
441 {
442 ip1.x = ip2.x;
443 }
444 else
445 {
446 ip0.x = ip2.x;
447 }
448
449 ip2.x = (ip1.x + ip0.x) / 2.;
450 }
451 intp.x = ip2.x;
452 intp.weight = 1. / Tr.Weight();
453 }
454 else if (LvlSet->Eval(Tr, ip0) > 0. && LvlSet->Eval(Tr, ip1) <= 1e-12)
455 {
456 intp.x = 1.;
457 intp.weight = 1. / Tr.Weight();
458 }
459 else if (LvlSet->Eval(Tr, ip1) > 0. && LvlSet->Eval(Tr, ip0) <= 1e-12)
460 {
461 intp.x = 0.;
462 intp.weight = 1. / Tr.Weight();
463 }
464 else
465 {
466 intp.x = .5;
467 intp.weight = 0.;
468 }
469
470 return intp.x;
471}
472
474{
477
479 ip0.x = 0.;
481 ip1.x = 1.;
482 Tr.SetIntPoint(&ip0);
483 if (LvlSet->Eval(Tr, ip0) * LvlSet->Eval(Tr, ip1) < 0.)
484 {
485 Vector tempX(ir.GetNPoints());
486 real_t length;
487 if (LvlSet->Eval(Tr, ip0) > 0.)
488 {
489 length = bisect(Tr, LvlSet);
490 for (int ip = 0; ip < ir.GetNPoints(); ip++)
491 {
492 IntegrationPoint &intp = ir.IntPoint(ip);
493 intp.x = ir2.IntPoint(ip).x * length;
494 intp.weight = ir2.IntPoint(ip).weight * length;
495 }
496 }
497 else
498 {
499 length = 1. - bisect(Tr, LvlSet);
500 for (int ip = 0; ip < ir.GetNPoints(); ip++)
501 {
502 IntegrationPoint &intp = ir.IntPoint(ip);
503 intp.x = bisect(Tr, LvlSet) + ir2.IntPoint(ip).x * length;
504 intp.weight = ir2.IntPoint(ip).weight * length;
505 }
506 }
507 }
508 else if (LvlSet->Eval(Tr, ip0) <= -tol_1
509 || LvlSet->Eval(Tr, ip1) <= -tol_1)
510 {
511 for (int ip = 0; ip < ir.GetNPoints(); ip++)
512 {
513 IntegrationPoint &intp = ir.IntPoint(ip);
514 intp.x = ir2.IntPoint(ip).x;
515 intp.weight = 0.;
516 }
517 }
518 else
519 {
520 ir = ir2;
521 }
522}
523
525{
526 int elem = Tr.ElementNo;
527 const Mesh* mesh = Tr.mesh;
528
529 const Element* me = mesh->GetElement(elem);
531 mesh->GetElementTransformation(elem, &Trafo);
532
534 Mat = 0.;
535 Vector RHS(nBasis);
536 RHS = 0.;
537 Vector ElemWeights(ir.GetNPoints());
538 ElemWeights = 0.;
539
540 bool element_int = false;
541 bool interior = true;
542 Array<bool> edge_int;
543
544 DenseMatrix PointA(me->GetNEdges(), Trafo.GetSpaceDim());
545 DenseMatrix PointB(me->GetNEdges(), Trafo.GetSpaceDim());
546 Vector edgelength(me->GetNEdges());
547
548 Array<int> verts;
549 mesh->GetElementVertices(elem, verts);
550
551 // find the edges that are intersected by the surface and inside the area
552 for (int edge = 0; edge < me->GetNEdges(); edge++)
553 {
554 enum class Layout {inside, intersected, outside};
555 Layout layout;
556
557 const int* vert = me->GetEdgeVertices(edge);
558 Vector pointA(Trafo.GetSpaceDim());
559 Vector pointB(Trafo.GetSpaceDim());
560 for (int d = 0; d < Trafo.GetSpaceDim(); d++)
561 {
562 pointA(d) = (Trafo.mesh->GetVertex(verts[vert[0]]))[d];
563 pointB(d) = (Trafo.mesh->GetVertex(verts[vert[1]]))[d];
564 }
565 Vector edgevec(Trafo.GetSpaceDim());
566 subtract(pointA, pointB, edgevec);
567 edgelength(edge) = edgevec.Norml2();
568
570 Trafo.TransformBack(pointA, ipA);
572 Trafo.TransformBack(pointB, ipB);
573
574 if (LvlSet->Eval(Trafo, ipA) < -tol_1
575 || LvlSet->Eval(Trafo, ipB) < -tol_1)
576 {
577 interior = false;
578 }
579
580 if (LvlSet->Eval(Trafo, ipA) > -tol_1
581 && LvlSet->Eval(Trafo, ipB) > -tol_1)
582 {
583 layout = Layout::inside;
584 }
585 else if (LvlSet->Eval(Trafo, ipA) > tol_2
586 && LvlSet->Eval(Trafo, ipB) <= 0.)
587 {
588 layout = Layout::intersected;
589 }
590 else if (LvlSet->Eval(Trafo, ipA) <= 0.
591 && LvlSet->Eval(Trafo, ipB) > tol_2)
592 {
593 layout = Layout::intersected;
594 Vector temp(pointA.Size());
595 temp = pointA;
596 pointA = pointB;
597 pointB = temp;
598 }
599 else
600 {
601 layout = Layout::outside;
602 }
603
604 // Store the end points of the (1D) intersected edge.
605 if (layout == Layout::intersected)
606 {
607 Vector pointC(pointA.Size());
608 Vector mid(pointA.Size());
609 pointC = pointA;
610 mid = pointC;
611 mid += pointB;
612 mid /= 2.;
613
615 Trafo.TransformBack(mid, ip);
616
617 while (LvlSet->Eval(Trafo, ip) > tol_1
618 || LvlSet->Eval(Trafo, ip) < -tol_1)
619 {
620 if (LvlSet->Eval(Trafo, ip) > tol_1)
621 {
622 pointC = mid;
623 }
624 else
625 {
626 pointB = mid;
627 }
628
629 mid = pointC;
630 mid += pointB;
631 mid /= 2.;
632 Trafo.TransformBack(mid, ip);
633 }
634 pointB = mid;
635 }
636 PointA.SetRow(edge, pointA);
637 PointB.SetRow(edge, pointB);
638
639 if ((layout == Layout::inside || layout == Layout::intersected))
640 {
641 edge_int.Append(true);
642 }
643 else
644 {
645 edge_int.Append(false);
646 }
647 }
648
649 // Integrate over the 1D edges.
650 for (int edge = 0; edge < me->GetNEdges(); edge++)
651 {
652 if (edge_int[edge] && !interior)
653 {
654 Vector point0(Trafo.GetSpaceDim());
655 Vector point1(Trafo.GetSpaceDim());
656 PointA.GetRow(edge, point0);
657 PointB.GetRow(edge, point1);
658
659 element_int = true;
660
662 2*Order+1);
663
664 Vector normal(Trafo.GetDimension());
665 normal = 0.;
666 if (edge == 0 || edge == 2)
667 {
668 normal(1) = 1.;
669 }
670 if (edge == 1 || edge == 3)
671 {
672 normal(0) = 1.;
673 }
674 if (edge == 0 || edge == 3)
675 {
676 normal *= -1.;
677 }
678
679 for (int ip = 0; ip < ir2->GetNPoints(); ip++)
680 {
681 Vector dist(Trafo.GetSpaceDim());
682 dist = point1;
683 dist -= point0;
684
685 Vector point(Trafo.GetSpaceDim());
686 point = dist;
687 point *= ir2->IntPoint(ip).x;
688 point += point0;
689
690 IntegrationPoint intpoint;
691 Trafo.TransformBack(point, intpoint);
692 Trafo.SetIntPoint(&intpoint);
693 DenseMatrix shapes;
694 OrthoBasis2D(intpoint, shapes);
695 Vector grad(Trafo.GetDimension());
696
697 for (int dof = 0; dof < nBasis; dof++)
698 {
699 shapes.GetRow(dof, grad);
700 RHS(dof) -= (grad * normal) * ir2->IntPoint(ip).weight
701 * dist.Norml2() / edgelength(edge);
702 }
703 }
704 }
705 }
706
707 // do integration over the area for integral over interface
708 if (element_int && !interior)
709 {
710 H1_FECollection fec(lsOrder, 2);
711 FiniteElementSpace fes(const_cast<Mesh*>(Tr.mesh), &fec);
712 GridFunction LevelSet(&fes);
713 LevelSet.ProjectCoefficient(*LvlSet);
714 mesh->GetElementTransformation(elem, &Trafo);
715
716 const FiniteElement* fe = fes.GetFE(elem);
717 Vector normal(Trafo.GetDimension());
718 Vector gradi(Trafo.GetDimension());
719 DenseMatrix dshape(fe->GetDof(), Trafo.GetDimension());
720 Array<int> dofs;
721 fes.GetElementDofs(elem, dofs);
722
723 for (int ip = 0; ip < ir.GetNPoints(); ip++)
724 {
725 Trafo.SetIntPoint(&(ir.IntPoint(ip)));
726
727 normal = 0.;
728 fe->CalcDShape(ir.IntPoint(ip), dshape);
729 for (int dof = 0; dof < fe->GetDof(); dof++)
730 {
731 dshape.GetRow(dof, gradi);
732 gradi *= LevelSet(dofs[dof]);
733 normal += gradi;
734 }
735 normal *= (-1. / normal.Norml2());
736
737 DenseMatrix shapes;
738 OrthoBasis2D(ir.IntPoint(ip), shapes);
739
740 for (int dof = 0; dof < nBasis; dof++)
741 {
742 Vector grad(Trafo.GetSpaceDim());
743 shapes.GetRow(dof, grad);
744 Mat(dof, ip) = (grad * normal);
745 }
746 }
747
748 // solve the underdetermined linear system
749 Vector temp(nBasis);
750 Vector temp2(ir.GetNPoints());
751 DenseMatrixSVD SVD(Mat, 'A', 'A');
752 SVD.Eval(Mat);
753 SVD.LeftSingularvectors().MultTranspose(RHS, temp);
754 temp2 = 0.;
755 for (int i = 0; i < nBasis; i++)
756 {
757 if (SVD.Singularvalue(i) > tol_1)
758 {
759 temp2(i) = temp(i) / SVD.Singularvalue(i);
760 }
761 }
762 SVD.RightSingularvectors().MultTranspose(temp2, ElemWeights);
763 }
764
765 // save the weights
766 for (int ip = 0; ip < ir.GetNPoints(); ip++)
767 {
768 IntegrationPoint& intp = ir.IntPoint(ip);
769 intp.weight = ElemWeights(ip);
770 }
771
772 mesh->GetElementTransformation(elem, &Trafo);
773}
774
776 const IntegrationRule* sir)
777{
778 int elem = Tr.ElementNo;
779 const Mesh* mesh = Tr.mesh;
780
781 const Element* me = mesh->GetElement(elem);
783 mesh->GetElementTransformation(elem, &Trafo);
784
785 Vector RHS(nBasisVolume);
786 RHS = 0.;
787 Vector ElemWeights(ir.GetNPoints());
788 ElemWeights = 0.;
789
790 bool element_int = false;
791 bool interior = true;
792 Array<bool> edge_int;
793
794 DenseMatrix PointA(me->GetNEdges(), Trafo.GetSpaceDim());
795 DenseMatrix PointB(me->GetNEdges(), Trafo.GetSpaceDim());
796 Vector edgelength(me->GetNEdges());
797
798 Array<int> verts;
799 mesh->GetElementVertices(elem, verts);
800
801 // find the edges that are intersected by he surface and inside the area
802 for (int edge = 0; edge < me->GetNEdges(); edge++)
803 {
804 enum class Layout {inside, intersected, outside};
805 Layout layout;
806
807 const int* vert = me->GetEdgeVertices(edge);
808 Vector pointA(Trafo.GetSpaceDim());
809 Vector pointB(Trafo.GetSpaceDim());
810 for (int d = 0; d < Trafo.GetSpaceDim(); d++)
811 {
812 pointA(d) = (Trafo.mesh->GetVertex(verts[vert[0]]))[d];
813 pointB(d) = (Trafo.mesh->GetVertex(verts[vert[1]]))[d];
814 }
815 Vector edgevec(Trafo.GetSpaceDim());
816 subtract(pointA, pointB, edgevec);
817 edgelength(edge) = edgevec.Norml2();
818
820 Trafo.TransformBack(pointA, ipA);
822 Trafo.TransformBack(pointB, ipB);
823
824 if (LvlSet->Eval(Trafo, ipA) < -tol_1
825 || LvlSet->Eval(Trafo, ipB) < -tol_1)
826 {
827 interior = false;
828 }
829
830 if (LvlSet->Eval(Trafo, ipA) > -tol_1
831 && LvlSet->Eval(Trafo, ipB) > -tol_1)
832 {
833 layout = Layout::inside;
834 }
835 else if (LvlSet->Eval(Trafo, ipA) > tol_2
836 && LvlSet->Eval(Trafo, ipB) <= 0.)
837 {
838 layout = Layout::intersected;
839 }
840 else if (LvlSet->Eval(Trafo, ipA) <= 0.
841 && LvlSet->Eval(Trafo, ipB) > tol_2)
842 {
843 layout = Layout::intersected;
844 Vector temp(pointA.Size());
845 temp = pointA;
846 pointA = pointB;
847 pointB = temp;
848 }
849 else
850 {
851 layout = Layout::outside;
852 }
853
854 if (layout == Layout::intersected)
855 {
856 Vector pointC(pointA.Size());
857 Vector mid(pointA.Size());
858 pointC = pointA;
859 mid = pointC;
860 mid += pointB;
861 mid /= 2.;
862
864 Trafo.TransformBack(mid, ip);
865
866 while (LvlSet->Eval(Trafo, ip) > tol_1
867 || LvlSet->Eval(Trafo, ip) < -tol_1)
868 {
869 if (LvlSet->Eval(Trafo, ip) > tol_1)
870 {
871 pointC = mid;
872 }
873 else
874 {
875 pointB = mid;
876 }
877
878 mid = pointC;
879 mid += pointB;
880 mid /= 2.;
881 Trafo.TransformBack(mid, ip);
882 }
883 pointB = mid;
884 }
885
886 PointA.SetRow(edge, pointA);
887 PointB.SetRow(edge, pointB);
888
889 if ((layout == Layout::inside || layout == Layout::intersected))
890 {
891 edge_int.Append(true);
892 }
893 else
894 {
895 edge_int.Append(false);
896 }
897 }
898
899 // do the integration over the edges
900 for (int edge = 0; edge < me->GetNEdges(); edge++)
901 {
902 if (edge_int[edge] && !interior)
903 {
904 Vector point0(Trafo.GetSpaceDim());
905 Vector point1(Trafo.GetSpaceDim());
906 PointA.GetRow(edge, point0);
907 PointB.GetRow(edge, point1);
908
909 element_int = true;
910
912 2*Order+1);
913 Vector normal(Trafo.GetDimension());
914 normal = 0.;
915 if (edge == 0 || edge == 2)
916 {
917 normal(1) = 1.;
918 }
919 if (edge == 1 || edge == 3)
920 {
921 normal(0) = 1.;
922 }
923 if (edge == 0 || edge == 3)
924 {
925 normal *= -1.;
926 }
927
928 for (int ip = 0; ip < ir2->GetNPoints(); ip++)
929 {
930 Vector dist(Trafo.GetSpaceDim());
931 dist = point1;
932 dist -= point0;
933
934 Vector point(Trafo.GetSpaceDim());
935 point = dist;
936 point *= ir2->IntPoint(ip).x;
937 point += point0;
938
939 IntegrationPoint intpoint;
940 Trafo.TransformBack(point, intpoint);
941 DenseMatrix shapes;
942 BasisAD2D(intpoint, shapes);
943 Vector adiv(Trafo.GetDimension());
944
945 for (int dof = 0; dof < nBasisVolume; dof++)
946 {
947 shapes.GetRow(dof, adiv);
948 RHS(dof) += (adiv * normal) * ir2->IntPoint(ip).weight
949 * dist.Norml2() / edgelength(edge);
950 }
951 }
952 }
953 }
954
955 // Integrate over the interface using the already computed surface rule, and
956 // solve the linear system for the weights.
957 if (element_int && !interior)
958 {
959 H1_FECollection fec(lsOrder, 2);
960 FiniteElementSpace fes(const_cast<Mesh*>(Tr.mesh), &fec);
961 GridFunction LevelSet(&fes);
962 LevelSet.ProjectCoefficient(*LvlSet);
963 mesh->GetElementTransformation(elem, &Trafo);
964
965 const FiniteElement* fe = fes.GetFE(elem);
966 Vector normal(Trafo.GetDimension());
967 Vector gradi(Trafo.GetDimension());
968 DenseMatrix dshape(fe->GetDof(), Trafo.GetDimension());
969 Array<int> dofs;
970 fes.GetElementDofs(elem, dofs);
971
972 for (int ip = 0; ip < sir->GetNPoints(); ip++)
973 {
974 Trafo.SetIntPoint(&(sir->IntPoint(ip)));
975
976 normal = 0.;
977 fe->CalcDShape(sir->IntPoint(ip), dshape);
978 for (int dof = 0; dof < fe->GetDof(); dof++)
979 {
980 dshape.GetRow(dof, gradi);
981 gradi *= LevelSet(dofs[dof]);
982 normal += gradi;
983 }
984 normal *= (-1. / normal.Norml2());
985
986 DenseMatrix shapes;
987 BasisAD2D(sir->IntPoint(ip), shapes);
988
989 for (int dof = 0; dof < nBasisVolume; dof++)
990 {
991 Vector adiv(2);
992 shapes.GetRow(dof, adiv);
993 RHS(dof) += (adiv * normal) * sir->IntPoint(ip).weight;
994 }
995 }
996
997 // solve the underdetermined linear system
998 Vector temp(nBasisVolume);
999 Vector temp2(ir.GetNPoints());
1000 temp2 = 0.;
1002 for (int i = 0; i < nBasisVolume; i++)
1003 {
1004 if (VolumeSVD->Singularvalue(i) > tol_1)
1005 {
1006 temp2(i) = temp(i) / VolumeSVD->Singularvalue(i);
1007 }
1008 }
1009 VolumeSVD->RightSingularvectors().MultTranspose(temp2, ElemWeights);
1010 }
1011
1012 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1013 {
1014 IntegrationPoint& intp = ir.IntPoint(ip);
1015 intp.weight = ElemWeights(ip);
1016 }
1017
1018 if (interior)
1019 {
1020 int qorder = 0;
1022 IntegrationRule ir2 = irs.Get(Trafo.GetGeometryType(), qorder);
1023 for (; ir2.GetNPoints() < ir.GetNPoints(); qorder++)
1024 {
1025 ir2 = irs.Get(Trafo.GetGeometryType(), qorder);
1026 }
1027 ir = ir2;
1028 }
1029
1030 mesh->GetElementTransformation(elem, &Trafo);
1031}
1032
1034{
1036
1037 int elem = Tr.ElementNo;
1038 const Mesh* mesh = Tr.mesh;
1039
1040 const Element* me = mesh->GetElement(elem);
1042 mesh->GetElementTransformation(elem, &Trafo);
1043
1045 Mat = 0.;
1046 Vector RHS(nBasis);
1047 RHS = 0.;
1048 Vector ElemWeights(ir.GetNPoints());
1049 ElemWeights = 0.;
1050
1051 // Does the element have a positive vertex?
1052 bool element_int = false;
1053 // Are all element vertices positive?
1054 bool interior = true;
1055
1056 Array<int> verts;
1057 mesh->GetElementVertices(elem, verts);
1058
1059 for (int face = 0; face < me->GetNFaces(); face++)
1060 {
1061 const int* vert = me->GetFaceVertices(face);
1062 Vector pointA(Trafo.GetSpaceDim());
1063 Vector pointB(Trafo.GetSpaceDim());
1064 Vector pointC(Trafo.GetSpaceDim());
1065 Vector pointD(Trafo.GetSpaceDim());
1066 for (int d = 0; d < Trafo.GetSpaceDim(); d++)
1067 {
1068 pointA(d) = (Trafo.mesh->GetVertex(verts[vert[0]]))[d];
1069 pointB(d) = (Trafo.mesh->GetVertex(verts[vert[1]]))[d];
1070 pointC(d) = (Trafo.mesh->GetVertex(verts[vert[2]]))[d];
1071 pointD(d) = (Trafo.mesh->GetVertex(verts[vert[3]]))[d];
1072 }
1073
1074 IntegrationPoint ipA;
1075 Trafo.TransformBack(pointA, ipA);
1076 IntegrationPoint ipB;
1077 Trafo.TransformBack(pointB, ipB);
1078 IntegrationPoint ipC;
1079 Trafo.TransformBack(pointC, ipC);
1080 IntegrationPoint ipD;
1081 Trafo.TransformBack(pointD, ipD);
1082
1083 if (LvlSet->Eval(Trafo, ipA) < -tol_1
1084 || LvlSet->Eval(Trafo, ipB) < -tol_1
1085 || LvlSet->Eval(Trafo, ipC) < -tol_1
1086 || LvlSet->Eval(Trafo, ipD) < -tol_1)
1087 {
1088 interior = false;
1089 }
1090
1091 if (LvlSet->Eval(Trafo, ipA) > -tol_1
1092 || LvlSet->Eval(Trafo, ipB) > -tol_1
1093 || LvlSet->Eval(Trafo, ipC) > -tol_1
1094 || LvlSet->Eval(Trafo, ipD) > -tol_1)
1095 {
1096 element_int = true;
1097 }
1098
1099 Array<int> faces;
1100 Array<int> cor;
1101 mesh->GetElementFaces(elem, faces, cor);
1102
1105 Trafo.mesh->GetFaceElementTransformations(faces[face], FTrans, Tr1, Tr2);
1106 FTrans.SetIntPoint(&(FaceIP[0]));
1107
1108 Vector normal(Trafo.GetDimension());
1109 normal = 0.;
1110 if (face == 0 || face == 5)
1111 {
1112 normal(2) = 1.;
1113 }
1114 if (face == 1 || face == 3)
1115 {
1116 normal(1) = 1.;
1117 }
1118 if (face == 2 || face == 4)
1119 {
1120 normal(0) = 1.;
1121 }
1122 if (face == 0 || face == 1 || face == 4)
1123 {
1124 normal *= -1.;
1125 }
1126
1127 for (int ip = 0; ip < FaceIP.Size(); ip++)
1128 {
1129 DenseMatrix shape;
1130 Vector point(3);
1131 IntegrationPoint ipoint;
1132 FTrans.Transform(FaceIP[ip], point);
1133 Trafo.TransformBack(point, ipoint);
1134 OrthoBasis3D(ipoint, shape);
1135
1136 for (int dof = 0; dof < nBasis; dof++)
1137 {
1138 Vector grad(Trafo.GetSpaceDim());
1139 shape.GetRow(dof, grad);
1140 RHS(dof) -= (grad * normal) * FaceWeights(ip, faces[face]);
1141 }
1142 }
1143 }
1144
1145 // If the element is intersected, form the matrix and solve for the weights.
1146 if (element_int && !interior)
1147 {
1148 H1_FECollection fec(lsOrder, 3);
1149 FiniteElementSpace fes(const_cast<Mesh*>(Tr.mesh), &fec);
1150 GridFunction LevelSet(&fes);
1151 LevelSet.ProjectCoefficient(*LvlSet);
1152 mesh->GetElementTransformation(elem, &Trafo);
1153
1154 const FiniteElement* fe = fes.GetFE(elem);
1155 Vector normal(Trafo.GetDimension());
1156 Vector gradi(Trafo.GetDimension());
1157 DenseMatrix dshape(fe->GetDof(), Trafo.GetDimension());
1158 Array<int> dofs;
1159 fes.GetElementDofs(elem, dofs);
1160
1161 // Form the matrix.
1162 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1163 {
1164 Trafo.SetIntPoint(&(ir.IntPoint(ip)));
1165
1166 normal = 0.;
1167 fe->CalcDShape(ir.IntPoint(ip), dshape);
1168 for (int dof = 0; dof < fe->GetDof(); dof++)
1169 {
1170 dshape.GetRow(dof, gradi);
1171 gradi *= LevelSet(dofs[dof]);
1172 normal += gradi;
1173 }
1174 normal *= (-1. / normal.Norml2());
1175
1176 DenseMatrix shapes;
1177 OrthoBasis3D(ir.IntPoint(ip), shapes);
1178
1179 for (int dof = 0; dof < nBasis; dof++)
1180 {
1181 Vector grad(Trafo.GetSpaceDim());
1182 shapes.GetRow(dof, grad);
1183 Mat(dof, ip) = (grad * normal);
1184 }
1185 }
1186
1187 // solve the underdetermined linear system
1188 Vector temp(nBasis);
1189 Vector temp2(ir.GetNPoints());
1190 DenseMatrixSVD SVD(Mat, 'A', 'A');
1191 SVD.Eval(Mat);
1192 SVD.LeftSingularvectors().MultTranspose(RHS, temp);
1193 temp2 = 0.;
1194 for (int i = 0; i < nBasis; i++)
1195 {
1196 if (SVD.Singularvalue(i) > tol_1)
1197 {
1198 temp2(i) = temp(i) / SVD.Singularvalue(i);
1199 }
1200 }
1201 SVD.RightSingularvectors().MultTranspose(temp2, ElemWeights);
1202 }
1203
1204 // scale the weights
1205 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1206 {
1207 IntegrationPoint& intp = ir.IntPoint(ip);
1208 intp.weight = ElemWeights(ip);
1209 }
1210
1211 mesh->GetElementTransformation(elem, &Trafo);
1212}
1213
1215 const IntegrationRule* sir)
1216{
1217 Order++;
1219 Order--;
1220
1221 int elem = Tr.ElementNo;
1222 const Mesh* mesh = Tr.mesh;
1223
1224 const Element* me = mesh->GetElement(elem);
1226 mesh->GetElementTransformation(elem, &Trafo);
1227
1228 Vector RHS(nBasisVolume);
1229 RHS = 0.;
1230 Vector ElemWeights(ir.GetNPoints());
1231 ElemWeights = 0.;
1232
1233 // Does the element have a positive vertex?
1234 bool element_int = false;
1235 // Are all element vertices positive?
1236 bool interior = true;
1237
1238 Array<int> verts;
1239 mesh->GetElementVertices(elem, verts);
1240
1241 for (int face = 0; face < me->GetNFaces(); face++)
1242 {
1243 const int* vert = me->GetFaceVertices(face);
1244 Vector pointA(Trafo.GetSpaceDim());
1245 Vector pointB(Trafo.GetSpaceDim());
1246 Vector pointC(Trafo.GetSpaceDim());
1247 Vector pointD(Trafo.GetSpaceDim());
1248 for (int d = 0; d < Trafo.GetSpaceDim(); d++)
1249 {
1250 pointA(d) = (Trafo.mesh->GetVertex(verts[vert[0]]))[d];
1251 pointB(d) = (Trafo.mesh->GetVertex(verts[vert[1]]))[d];
1252 pointC(d) = (Trafo.mesh->GetVertex(verts[vert[2]]))[d];
1253 pointD(d) = (Trafo.mesh->GetVertex(verts[vert[3]]))[d];
1254 }
1255
1256 IntegrationPoint ipA;
1257 Trafo.TransformBack(pointA, ipA);
1258 IntegrationPoint ipB;
1259 Trafo.TransformBack(pointB, ipB);
1260 IntegrationPoint ipC;
1261 Trafo.TransformBack(pointC, ipC);
1262 IntegrationPoint ipD;
1263 Trafo.TransformBack(pointD, ipD);
1264
1265 if (LvlSet->Eval(Trafo, ipA) < -tol_1
1266 || LvlSet->Eval(Trafo, ipB) < -tol_1
1267 || LvlSet->Eval(Trafo, ipC) < -tol_1
1268 || LvlSet->Eval(Trafo, ipD) < -tol_1)
1269 {
1270 interior = false;
1271 }
1272
1273 if (LvlSet->Eval(Trafo, ipA) > -tol_1
1274 || LvlSet->Eval(Trafo, ipB) > -tol_1
1275 || LvlSet->Eval(Trafo, ipC) > -tol_1
1276 || LvlSet->Eval(Trafo, ipD) > -tol_1)
1277 {
1278 element_int = true;
1279 }
1280
1281 Array<int> faces;
1282 Array<int> cor;
1283 mesh->GetElementFaces(elem, faces, cor);
1284
1287 Trafo.mesh->GetFaceElementTransformations(faces[face], FTrans, Tr1, Tr2);
1288
1289 FTrans.SetIntPoint(&(FaceIP[0]));
1290
1291 Vector normal(Trafo.GetDimension());
1292 normal = 0.;
1293 if (face == 0 || face == 5)
1294 {
1295 normal(2) = 1.;
1296 }
1297 if (face == 1 || face == 3)
1298 {
1299 normal(1) = 1.;
1300 }
1301 if (face == 2 || face == 4)
1302 {
1303 normal(0) = 1.;
1304 }
1305 if (face == 0 || face == 1 || face == 4)
1306 {
1307 normal *= -1.;
1308 }
1309
1310 for (int ip = 0; ip < FaceIP.Size(); ip++)
1311 {
1312 DenseMatrix shape;
1313 Vector point(3);
1314 IntegrationPoint ipoint;
1315 FTrans.Transform(FaceIP[ip], point);
1316 Trafo.TransformBack(point, ipoint);
1317 BasisAD3D(ipoint, shape);
1318
1319 for (int dof = 0; dof < nBasisVolume; dof++)
1320 {
1321 Vector adiv(Trafo.GetSpaceDim());
1322 shape.GetRow(dof, adiv);
1323 RHS(dof) += (adiv * normal) * FaceWeights(ip, faces[face]);
1324 }
1325 }
1326 }
1327
1328 // If the element is intersected, integrate over the cut surface (using the
1329 // already computed rule) and solve the matrix for the weights.
1330 if (element_int && !interior)
1331 {
1332 H1_FECollection fec(lsOrder, 3);
1333 FiniteElementSpace fes(const_cast<Mesh*>(Tr.mesh), &fec);
1334 GridFunction LevelSet(&fes);
1335 LevelSet.ProjectCoefficient(*LvlSet);
1336 mesh->GetElementTransformation(elem, &Trafo);
1337
1338 const FiniteElement* fe = fes.GetFE(elem);
1339 Vector normal(Trafo.GetDimension());
1340 Vector gradi(Trafo.GetDimension());
1341 DenseMatrix dshape(fe->GetDof(), Trafo.GetDimension());
1342 Array<int> dofs;
1343 fes.GetElementDofs(elem, dofs);
1344
1345 // Integrate over the cut surface using the already computed rule.
1346 for (int ip = 0; ip < sir->GetNPoints(); ip++)
1347 {
1348 Trafo.SetIntPoint(&(sir->IntPoint(ip)));
1349
1350 normal = 0.;
1351 fe->CalcDShape(sir->IntPoint(ip), dshape);
1352 for (int dof = 0; dof < fe->GetDof(); dof++)
1353 {
1354 dshape.GetRow(dof, gradi);
1355 gradi *= LevelSet(dofs[dof]);
1356 normal += gradi;
1357 }
1358 normal *= (-1. / normal.Norml2());
1359
1360 DenseMatrix shapes;
1361 BasisAD3D(sir->IntPoint(ip), shapes);
1362
1363 for (int dof = 0; dof < nBasisVolume; dof++)
1364 {
1365 Vector adiv(Trafo.GetSpaceDim());
1366 shapes.GetRow(dof, adiv);
1367 RHS(dof) += (adiv * normal) * sir->IntPoint(ip).weight;
1368 }
1369 }
1370
1371 // solve the underdetermined linear system
1372 Vector temp(nBasisVolume);
1373 Vector temp2(ir.GetNPoints());
1375 temp2 = 0.;
1376 for (int i = 0; i < nBasisVolume; i++)
1377 if (VolumeSVD->Singularvalue(i) > tol_1)
1378 {
1379 temp2(i) = temp(i) / VolumeSVD->Singularvalue(i);
1380 }
1381 VolumeSVD->RightSingularvectors().MultTranspose(temp2, ElemWeights);
1382 }
1383
1384 // scale the weights
1385 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1386 {
1387 IntegrationPoint& intp = ir.IntPoint(ip);
1388 intp.weight = ElemWeights(ip);
1389 }
1390
1391 // Fully inside the subdomain -> standard integration.
1392 if (interior)
1393 {
1394 int qorder = 0;
1396 IntegrationRule ir2 = irs.Get(Trafo.GetGeometryType(), qorder);
1397 for (; ir2.GetNPoints() < ir.GetNPoints(); qorder++)
1398 {
1399 ir2 = irs.Get(Trafo.GetGeometryType(), qorder);
1400 }
1401 ir = ir2;
1402 }
1403
1404 mesh->GetElementTransformation(elem, &Trafo);
1405}
1406
1408 DenseMatrix& shape)
1409{
1410 shape.SetSize(nBasis, 2);
1411
1412 Vector X(2);
1413 X(0) = -1. + 2. * ip.x;
1414 X(1) = -1. + 2. * ip.y;
1415
1416 for (int c = 0; c <= Order; c++)
1417 {
1418 Vector a(2);
1419 a = 0.;
1420 a(1) = pow(X(0), (real_t)(c));
1421
1422 Vector b(2);
1423 b = 0.;
1424 b(0) = pow(X(1), (real_t)(c));
1425
1426 shape.SetRow(2 * c, a);
1427 shape.SetRow(2 * c + 1, b);
1428 }
1429
1430 Poly_1D poly;
1431 int count = 2 * Order+ 2;
1432 for (int c = 1; c <= Order; c++)
1433 {
1434 const int* factorial = poly.Binom(c);
1435 for (int expo = c; expo > 0; expo--)
1436 {
1437 Vector a(2);
1438 a(0) = (real_t)(factorial[expo]) * pow(X(0), (real_t)(expo))
1439 * pow(X(1), (real_t)(c - expo));
1440 a(1) = -1. * (real_t)(factorial[expo - 1])
1441 * pow(X(0), (real_t)(expo - 1))
1442 * pow(X(1), (real_t)(c - expo + 1));
1443
1444 shape.SetRow(count, a);
1445 count++;
1446 }
1447 }
1448}
1449
1451 DenseMatrix& shape)
1452{
1454
1455 shape.SetSize(nBasis, 2);
1456
1457 // evaluate basis in the point
1458 DenseMatrix preshape(nBasis, 2);
1459 DivFreeBasis2D(ip, shape);
1460
1461 // evaluate basis for quadrature points
1462 DenseTensor shapeMFN(nBasis, 2, ir_->GetNPoints());
1463 for (int p = 0; p < ir_->GetNPoints(); p++)
1464 {
1465 DenseMatrix shapeN(nBasis, 2);
1466 DivFreeBasis2D(ir_->IntPoint(p), shapeN);
1467 for (int i = 0; i < nBasis; i++)
1468 for (int j = 0; j < 2; j++)
1469 {
1470 shapeMFN(i, j, p) = shapeN(i, j);
1471 }
1472 }
1473
1474 // do modified Gram-Schmidt orthogonalization
1475 for (int count = 1; count < nBasis; count++)
1476 {
1477 mGSStep(shape, shapeMFN, count);
1478 }
1479}
1480
1482 DenseMatrix& shape)
1483{
1484 Vector X(3);
1485 X(0) = -1. + 2. * ip.x;
1486 X(1) = -1. + 2. * ip.y;
1487 X(2) = -1. + 2. * ip.z;
1488
1490}
1491
1493 int step)
1494{
1496
1497 for (int count = step; count < shape.Height(); count++)
1498 {
1499 real_t den = 0.;
1500 real_t num = 0.;
1501
1502 for (int ip = 0; ip < ir_->GetNPoints(); ip++)
1503 {
1504 Vector u(2);
1505 Vector v(2);
1506
1507 shapeMFN(ip).GetRow(count, u);
1508 shapeMFN(ip).GetRow(step - 1, v);
1509
1510 den += v * v * ir_->IntPoint(ip).weight;
1511 num += u * v * ir_->IntPoint(ip).weight;
1512 }
1513
1514 real_t coeff = num / den;
1515
1516 Vector s(2);
1517 Vector t(2);
1518 shape.GetRow(step - 1, s);
1519 shape.GetRow(count, t);
1520 s *= coeff;
1521 t += s;
1522 shape.SetRow(count, t);
1523
1524 for (int ip = 0; ip < ir_->GetNPoints(); ip++)
1525 {
1526 shapeMFN(ip).GetRow(step - 1, s);
1527 shapeMFN(ip).GetRow(count, t);
1528 s *= coeff;
1529 t += s;
1530 shapeMFN(ip).SetRow(count, t);
1531 }
1532 }
1533}
1534
1536{
1537 shape.SetSize(nBasisVolume);
1538
1539 Vector X(2);
1540 X(0) = -1. + 2. * ip.x;
1541 X(1) = -1. + 2. * ip.y;
1542
1543 int count = 0;
1544 for (int c = 0; c <= Order; c++)
1545 {
1546 for (int expo = 0; expo <= c; expo++)
1547 {
1548 shape(count) = pow(X(0), (real_t)(expo))
1549 * pow(X(1), (real_t)(c - expo));
1550 count++;
1551 }
1552 }
1553}
1554
1556 DenseMatrix& shape)
1557{
1558 shape.SetSize(nBasisVolume, 2);
1559
1560 Vector X(2);
1561 X(0) = -1. + 2. * ip.x;
1562 X(1) = -1. + 2. * ip.y;
1563
1564 int count = 0;
1565 for (int c = 0; c <= Order; c++)
1566 {
1567 for (int expo = 0; expo <= c; expo++)
1568 {
1569 shape(count, 0) = .25 * pow(X(0), (real_t)(expo + 1))
1570 * pow(X(1), (real_t)(c - expo))
1571 / (real_t)(expo + 1);
1572 shape(count, 1) = .25 * pow(X(0), (real_t)(expo))
1573 * pow(X(1), (real_t)(c - expo + 1))
1574 / (real_t)(c - expo + 1);
1575 count++;
1576 }
1577 }
1578}
1579
1581{
1582 shape.SetSize(nBasisVolume);
1583
1584 Vector X(3);
1585 X(0) = -1. + 2. * ip.x;
1586 X(1) = -1. + 2. * ip.y;
1587 X(2) = -1. + 2. * ip.z;
1588
1589 int count = 0;
1590 for (int c = 0; c <= Order; c++)
1591 for (int expo = 0; expo <= c; expo++)
1592 for (int expo2 = 0; expo2 <= c - expo; expo2++)
1593 {
1594 shape(count) = pow(X(0), (real_t)(expo))
1595 * pow(X(1), (real_t)(expo2))
1596 * pow(X(2), (real_t)(c - expo - expo2));
1597 count++;
1598 }
1599}
1600
1602 DenseMatrix& shape)
1603{
1604 shape.SetSize(nBasisVolume, 3);
1605
1606 Vector X(3);
1607 X(0) = -1. + 2. * ip.x;
1608 X(1) = -1. + 2. * ip.y;
1609 X(2) = -1. + 2. * ip.z;
1610
1611 int count = 0;
1612 for (int c = 0; c <= Order; c++)
1613 for (int expo = 0; expo <= c; expo++)
1614 for (int expo2 = 0; expo2 <= c - expo; expo2++)
1615 {
1616 shape(count, 0) = pow(X(0), (real_t)(expo + 1))
1617 * pow(X(1), (real_t)(expo2))
1618 * pow(X(2), (real_t)(c - expo - expo2))
1619 / (6. * (real_t)(expo + 1));
1620 shape(count, 1) = pow(X(0), (real_t)(expo))
1621 * pow(X(1), (real_t)(expo2 + 1))
1622 * pow(X(2), (real_t)(c - expo - expo2))
1623 / (6. * (real_t)(expo2 + 1));;
1624 shape(count, 2) = pow(X(0), (real_t)(expo))
1625 * pow(X(1), (real_t)(expo2))
1626 * pow(X(2), (real_t)(c - expo - expo2 + 1))
1627 / (6. * (real_t)(c - expo + expo2 + 1));;
1628 count++;
1629 }
1630}
1631
1633{
1634 dim = -1;
1635 nBasis = -1;
1636 nBasisVolume = -1;
1637 delete VolumeSVD;
1638 VolumeSVD = NULL;
1639 FaceIP.DeleteAll();
1640 FaceWeights = 0.;
1641 FaceWeightsComp = 0.;
1642}
1643
1645{
1646 if (order != Order) { Clear(); }
1647 Order = order;
1648}
1649
1651 IntegrationRule& result)
1652{
1653 if (nBasis == -1 || dim != Tr.GetDimension())
1654 {
1655 Clear();
1657 }
1658
1659 if (Tr.GetDimension() == 3)
1660 {
1661 FaceIP.DeleteAll();
1662 FaceWeights = 0.;
1663 FaceWeightsComp = 0.;
1664 }
1665
1666 if (Tr.GetDimension() == 1)
1667 {
1669 }
1670 else if (Tr.GetDimension() == 2)
1671 {
1673 }
1674 else if (Tr.GetDimension() == 3)
1675 {
1677 }
1678
1679 result.SetSize(ir.GetNPoints());
1680 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1681 {
1682 result.IntPoint(ip).index = ip;
1683 IntegrationPoint &intp = result.IntPoint(ip);
1684 intp.x = ir.IntPoint(ip).x;
1685 intp.y = ir.IntPoint(ip).y;
1686 intp.z = ir.IntPoint(ip).z;
1687 intp.weight = ir.IntPoint(ip).weight;
1688 }
1689}
1690
1692 IntegrationRule& result,
1693 const IntegrationRule* sir)
1694{
1695 if (nBasis == -1 || nBasisVolume == -1 || dim != Tr.GetDimension())
1696 {
1697 Clear();
1699 }
1700
1701 if (Tr.GetDimension() == 3)
1702 {
1703 FaceIP.DeleteAll();
1704 FaceWeights = 0.;
1705 FaceWeightsComp = 0.;
1706 }
1707
1708 IntegrationRule SIR;
1709
1710 if (Tr.GetDimension() == 1)
1711 {
1712 Clear();
1714 }
1715 else if (sir == NULL)
1716 {
1717 Order++;
1719 Order--;
1720 }
1721 else if (sir->GetOrder() - 1 != ir.GetOrder())
1722 {
1723 Order++;
1725 Order--;
1726 }
1727 else { SIR = *sir; }
1728
1729 if (Tr.GetDimension() == 1)
1730 {
1732 }
1733 else if (Tr.GetDimension() == 2)
1734 {
1735 ComputeVolumeWeights2D(Tr, &SIR);
1736 }
1737 else if (Tr.GetDimension() == 3)
1738 {
1739 ComputeVolumeWeights3D(Tr, &SIR);
1740 }
1741
1742 result.SetSize(ir.GetNPoints());
1743 for (int ip = 0; ip < ir.GetNPoints(); ip++)
1744 {
1745 result.IntPoint(ip).index = ip;
1746 IntegrationPoint &intp = result.IntPoint(ip);
1747 intp.x = ir.IntPoint(ip).x;
1748 intp.y = ir.IntPoint(ip).y;
1749 intp.z = ir.IntPoint(ip).z;
1750 intp.weight = ir.IntPoint(ip).weight;
1751 }
1752}
1753
1755 const IntegrationRule &sir,
1756 Vector &weights)
1757{
1758 if (nBasis == -1 || dim != Tr.GetDimension())
1759 {
1760 Clear();
1762 }
1763
1764 weights.SetSize(sir.GetNPoints());
1765 weights = 0.0;
1766
1767 bool computeweights = false;
1768 for (int ip = 0; ip < sir.GetNPoints(); ip++)
1769 {
1770 if (sir.IntPoint(ip).weight != 0.)
1771 {
1772 computeweights = true;
1773 }
1774 }
1775
1776 if (Tr.GetDimension() > 1 && computeweights)
1777 {
1778 int elem = Tr.ElementNo;
1779 const Mesh* mesh = Tr.mesh;
1781 FiniteElementSpace fes(const_cast<Mesh*>(Tr.mesh), &fec);
1782 GridFunction LevelSet(&fes);
1783 LevelSet.ProjectCoefficient(*LvlSet);
1784
1786 mesh->GetElementTransformation(elem, &Trafo);
1787
1788 const FiniteElement* fe = fes.GetFE(elem);
1789 Vector normal(Tr.GetDimension());
1790 Vector normal2(Tr.GetSpaceDim());
1791 Vector gradi(Tr.GetDimension());
1792 DenseMatrix dshape(fe->GetDof(), Tr.GetDimension());
1793 Array<int> dofs;
1794 fes.GetElementDofs(elem, dofs);
1795
1796 for (int ip = 0; ip < sir.GetNPoints(); ip++)
1797 {
1798 Trafo.SetIntPoint(&(sir.IntPoint(ip)));
1799 LevelSet.GetGradient(Trafo, normal2);
1800 real_t normphys = normal2.Norml2();
1801
1802 normal = 0.;
1803 fe->CalcDShape(sir.IntPoint(ip), dshape);
1804 for (int dof = 0; dof < fe->GetDof(); dof++)
1805 {
1806 dshape.GetRow(dof, gradi);
1807 gradi *= LevelSet(dofs[dof]);
1808 normal += gradi;
1809 }
1810 real_t normref = normal.Norml2();
1811 normal *= (-1. / normal.Norml2());
1812
1813 weights(ip) = normphys / normref;
1814 }
1815 }
1816}
1817
1818#endif // MFEM_USE_LAPACK
1819
1820}
virtual void GetVolumeIntegrationRule(ElementTransformation &Tr, IntegrationRule &result, const IntegrationRule *sir=nullptr) override
Construct a cut-volume IntegrationRule.
virtual void GetSurfaceIntegrationRule(ElementTransformation &Tr, IntegrationRule &result) override
Construct a cut-surface IntegrationRule.
virtual void GetSurfaceWeights(ElementTransformation &Tr, const IntegrationRule &sir, Vector &weights) override
Compute transformation quadrature weights for surface integration.
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:758
int Size() const
Return the logical size of the array.
Definition array.hpp:147
void DeleteAll()
Delete the whole array.
Definition array.hpp:925
int Append(const T &el)
Append element 'el' to array, resize if necessary.
Definition array.hpp:830
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
static constexpr real_t tol_1
virtual void SetLevelSetProjectionOrder(int order)
Coefficient * LvlSet
The zero level set of this Coefficient defines the cut surface.
int lsOrder
Space order for the LS projection.
int Order
Order of the IntegrationRule.
static constexpr real_t tol_2
virtual void SetOrder(int order)
Change the order of the constructed IntegrationRule.
Class for Singular Value Decomposition of a DenseMatrix.
Definition densemat.hpp:964
DenseMatrix & RightSingularvectors()
Return right singular vectors.
DenseMatrix & LeftSingularvectors()
Return left singular vectors.
real_t Singularvalue(int i)
Return specific singular value.
void Eval(DenseMatrix &M)
Evaluate the SVD.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:125
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:172
void SetRow(int r, const real_t *row)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:108
void GetRow(int r, Vector &row) const
Rank 3 tensor (array of matrices)
const Mesh * mesh
The Mesh object containing the element.
Definition eltrans.hpp:97
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition eltrans.hpp:175
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:158
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:144
int GetDimension() const
Return the topological dimension of the reference element.
Definition eltrans.hpp:178
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition eltrans.hpp:106
Abstract data type element.
Definition element.hpp:29
virtual MFEM_DEPRECATED int GetNFaces(int &nFaceVertices) const =0
virtual const int * GetFaceVertices(int fi) const =0
virtual int GetNEdges() const =0
virtual const int * GetEdgeVertices(int) const =0
A specialized ElementTransformation class representing a face and its two neighboring elements.
Definition eltrans.hpp:750
void SetIntPoint(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition eltrans.cpp:609
void Transform(const IntegrationPoint &, Vector &) override
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition eltrans.cpp:663
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:244
DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element 'elem'. The returned indices are offsets into an ldo...
Definition fespace.cpp:3513
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Definition fespace.cpp:3835
Abstract class for all finite elements.
Definition fe_base.hpp:244
int GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:321
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:400
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition fe_base.hpp:334
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
void GetGradient(ElementTransformation &tr, Vector &grad) const
Gradient of a scalar function at a quadrature point.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:275
Class for integration point with weight.
Definition intrules.hpp:35
void Set(const real_t *p, const int dim)
Definition intrules.hpp:46
void Set2w(const real_t x1, const real_t x2, const real_t w)
Definition intrules.hpp:84
Class for an integration rule - an Array of IntegrationPoint.
Definition intrules.hpp:100
int GetOrder() const
Returns the order of the integration rule.
Definition intrules.hpp:249
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition intrules.hpp:256
void SetOrder(const int order)
Sets the order of the integration rule. This is only for keeping order information,...
Definition intrules.hpp:253
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition intrules.hpp:259
Container class for integration rules.
Definition intrules.hpp:422
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
A standard isoparametric element transformation.
Definition eltrans.hpp:629
int TransformBack(const Vector &v, IntegrationPoint &ip, const real_t phys_rel_tol=tol_0) override
Transform a point pt from physical space to a point ip in reference space and optionally can set a so...
Definition eltrans.hpp:717
int GetSpaceDim() const override
Get the dimension of the target (physical) space.
Definition eltrans.hpp:709
Mesh data type.
Definition mesh.hpp:64
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:1004
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1508
int AddQuad(int v1, int v2, int v3, int v4, int attr=1)
Adds a quadrilateral to the mesh given by 4 vertices v1 through v4.
Definition mesh.cpp:1958
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1339
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition mesh.hpp:1291
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
Definition mesh.cpp:1873
void GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
Definition mesh.cpp:2404
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
Definition mesh.cpp:7514
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1219
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition mesh.hpp:1526
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1321
Class for subdomain IntegrationRules by means of moment-fitting.
void BasisAD2D(const IntegrationPoint &ip, DenseMatrix &shape)
Antiderivatives of the monomial basis on the element [-1,1]^2.
void mGSStep(DenseMatrix &shape, DenseTensor &shapeMFN, int step)
A step of the modified Gram-Schmidt algorithm.
void OrthoBasis3D(const IntegrationPoint &ip, DenseMatrix &shape)
A orthogonalized divergence free basis on the element [-1,1]^3.
Vector FaceWeightsComp
Indicates the already computed face IntegrationRules.
void ComputeVolumeWeights1D(ElementTransformation &Tr)
Compute the 1D quadrature weights.
Array< IntegrationPoint > FaceIP
Array of face integration points.
int nBasis
Number of divergence-free basis functions for surface integration.
void ComputeSurfaceWeights2D(ElementTransformation &Tr)
Compute 2D quadrature weights.
void Init(int order, Coefficient &levelset, int lsO)
Initialize the MomentFittingIntRules.
void Basis2D(const IntegrationPoint &ip, Vector &shape)
Monomial basis on the element [-1,1]^2.
void InitSurface(int order, Coefficient &levelset, int lsO, ElementTransformation &Tr)
Initialization for surface IntegrationRule.
void InitVolume(int order, Coefficient &levelset, int lsO, ElementTransformation &Tr)
Initialization for volume IntegrationRule.
DenseMatrix FaceWeights
Column-wise Matrix for the face quadrature weights.
void ComputeSurfaceWeights1D(ElementTransformation &Tr)
Compute 1D quadrature weights.
void SetOrder(int order) override
Change the order of the constructed IntegrationRule.
int dim
Space Dimension of the element.
void OrthoBasis2D(const IntegrationPoint &ip, DenseMatrix &shape)
A orthogonalized divergence free basis on the element [-1,1]^2.
void ComputeVolumeWeights2D(ElementTransformation &Tr, const IntegrationRule *sir)
Compute the 2D quadrature weights.
DenseMatrixSVD * VolumeSVD
SVD of the matrix for volumetric IntegrationRules.
void Basis3D(const IntegrationPoint &ip, Vector &shape)
Monomial basis on the element [-1,1]^3.
void ComputeFaceWeights(ElementTransformation &Tr)
Compute the IntegrationRules on the faces.
void GetSurfaceIntegrationRule(ElementTransformation &Tr, IntegrationRule &result) override
Construct a cut-surface IntegrationRule.
void GetSurfaceWeights(ElementTransformation &Tr, const IntegrationRule &sir, Vector &weights) override
Compute transformation quadrature weights for surface integration.
void ComputeSurfaceWeights3D(ElementTransformation &Tr)
Compute 2D quadrature weights.
void DivFreeBasis2D(const IntegrationPoint &ip, DenseMatrix &shape)
A divergence free basis on the element [-1,1]^2.
IntegrationRule ir
IntegrationRule representing the reused IntegrationPoints.
int nBasisVolume
Number of basis functions for volume integration.
void Clear()
Clear stored data of the MomentFittingIntRules.
void BasisAD3D(const IntegrationPoint &ip, DenseMatrix &shape)
Antiderivatives of the monomial basis on the element [-1,1]^3.
void GetVolumeIntegrationRule(ElementTransformation &Tr, IntegrationRule &result, const IntegrationRule *sir=nullptr) override
Construct a cut-volume IntegrationRule.
void ComputeVolumeWeights3D(ElementTransformation &Tr, const IntegrationRule *sir)
Compute the 3D quadrature weights.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
Class for computing 1D special polynomials and their associated basis functions.
Definition fe_base.hpp:988
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "pchoose k" for k=0,...
Definition fe_base.cpp:2044
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_pos.cpp:25
Vector data type.
Definition vector.hpp:82
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:931
int Size() const
Returns the size of the vector.
Definition vector.hpp:226
void SetSize(int s)
Resize the vector to size s.
Definition vector.hpp:558
int dim
Definition ex24.cpp:53
real_t lvlset(const Vector &X)
Level-set function defining the implicit interface.
Definition ex38.cpp:49
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
void GetDivFree3DBasis(const Vector &X, DenseMatrix &shape, int Order)
3 dimensional divergence free basis functions on [-1,1]^3
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:548
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition vector.cpp:547
float real_t
Definition config.hpp:43
double bisect(ElementTransformation &Tr, Coefficient *LvlSet)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:492
real_t p(const Vector &x, real_t t)