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