24 MFEM_VERIFY(order > 0,
"Invalid input");
30 MFEM_VERIFY(order > 0,
"Invalid input");
39 GenerateLSVector(Tr,
LvlSet);
45 LevelSet2D ls(pe,lsvec);
46 auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
50 for (
size_t i=0; i<q.nodes.size(); i++)
53 ip.
Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
58 LevelSet3D ls(pe,lsvec);
59 auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
64 for (
size_t i=0; i<q.nodes.size(); i++)
67 ip.
Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
77 GenerateLSVector(Tr,
LvlSet);
83 LevelSet2D ls(pe,lsvec);
84 auto q = Algoim::quadGen<2>(ls,Algoim::BoundingBox<real_t,2>(0.0,1.0),
88 for (
size_t i=0; i<q.nodes.size(); i++)
91 ip.
Set2w(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].w);
96 LevelSet3D ls(pe,lsvec);
97 auto q = Algoim::quadGen<3>(ls,Algoim::BoundingBox<real_t,3>(0.0,1.0),
100 result.
SetSize(q.nodes.size());
102 for (
size_t i=0; i<q.nodes.size(); i++)
105 ip.
Set(q.nodes[i].x(0),q.nodes[i].x(1),q.nodes[i].x(2),q.nodes[i].w);
115 GenerateLSVector(Tr,
LvlSet);
149 if (currentLvlSet==
lvlset)
167 pe=
new H1Pos_QuadrilateralElement(
lsOrder);
168 le=
new H1_QuadrilateralElement(
lsOrder);
172 pe=
new H1Pos_HexahedronElement(
lsOrder);
173 le=
new H1_HexahedronElement(
lsOrder);
177 MFEM_ABORT(
"Currently MFEM + Algoim supports only quads and hexes.");
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++)
191 const IntegrationPoint &ip = ir.IntPoint(i);
193 lsfun(i)=
lvlset->Eval(Tr,ip);
200#ifdef MFEM_USE_LAPACK
205 Init(order, levelset, lsO);
284 Mat.SetCol(ip, shape);
313 for (
int face = 0; face < me->
GetNFaces(); face++)
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];
333 Mesh local_mesh(2,4,1,0,3);
351 for (
int ip = 0; ip < FaceRule.
GetNPoints(); ip++)
367 for (
int ip = 0; ip < FaceRule.
GetNPoints(); ip++)
402 ip2.
x = (ip1.
x + ip0.
x) / 2.;
433 if (LvlSet->
Eval(Tr, ip0) * LvlSet->
Eval(Tr, ip1) < 0.)
437 while (LvlSet->
Eval(Tr, ip2) > 1e-12
438 || LvlSet->
Eval(Tr, ip2) < -1e-12)
440 if (LvlSet->
Eval(Tr, ip0) * LvlSet->
Eval(Tr, ip2) < 0.)
449 ip2.
x = (ip1.
x + ip0.
x) / 2.;
454 else if (LvlSet->
Eval(Tr, ip0) > 0. && LvlSet->
Eval(Tr, ip1) <= 1e-12)
459 else if (LvlSet->
Eval(Tr, ip1) > 0. && LvlSet->
Eval(Tr, ip0) <= 1e-12)
540 bool element_int =
false;
541 bool interior =
true;
552 for (
int edge = 0; edge < me->
GetNEdges(); edge++)
554 enum class Layout {inside, intersected, outside};
567 edgelength(edge) = edgevec.
Norml2();
583 layout = Layout::inside;
588 layout = Layout::intersected;
593 layout = Layout::intersected;
601 layout = Layout::outside;
605 if (layout == Layout::intersected)
636 PointA.
SetRow(edge, pointA);
637 PointB.
SetRow(edge, pointB);
639 if ((layout == Layout::inside || layout == Layout::intersected))
650 for (
int edge = 0; edge < me->
GetNEdges(); edge++)
652 if (edge_int[edge] && !interior)
656 PointA.
GetRow(edge, point0);
657 PointB.
GetRow(edge, point1);
666 if (edge == 0 || edge == 2)
670 if (edge == 1 || edge == 3)
674 if (edge == 0 || edge == 3)
679 for (
int ip = 0; ip < ir2->
GetNPoints(); ip++)
697 for (
int dof = 0; dof <
nBasis; dof++)
701 * dist.
Norml2() / edgelength(edge);
708 if (element_int && !interior)
729 for (
int dof = 0; dof < fe->
GetDof(); dof++)
731 dshape.
GetRow(dof, gradi);
732 gradi *= LevelSet(dofs[dof]);
735 normal *= (-1. / normal.
Norml2());
740 for (
int dof = 0; dof <
nBasis; dof++)
744 Mat(dof, ip) = (grad * normal);
755 for (
int i = 0; i <
nBasis; i++)
769 intp.
weight = ElemWeights(ip);
790 bool element_int =
false;
791 bool interior =
true;
802 for (
int edge = 0; edge < me->
GetNEdges(); edge++)
804 enum class Layout {inside, intersected, outside};
817 edgelength(edge) = edgevec.
Norml2();
833 layout = Layout::inside;
838 layout = Layout::intersected;
843 layout = Layout::intersected;
851 layout = Layout::outside;
854 if (layout == Layout::intersected)
886 PointA.
SetRow(edge, pointA);
887 PointB.
SetRow(edge, pointB);
889 if ((layout == Layout::inside || layout == Layout::intersected))
900 for (
int edge = 0; edge < me->
GetNEdges(); edge++)
902 if (edge_int[edge] && !interior)
906 PointA.
GetRow(edge, point0);
907 PointB.
GetRow(edge, point1);
915 if (edge == 0 || edge == 2)
919 if (edge == 1 || edge == 3)
923 if (edge == 0 || edge == 3)
928 for (
int ip = 0; ip < ir2->
GetNPoints(); ip++)
949 * dist.
Norml2() / edgelength(edge);
957 if (element_int && !interior)
972 for (
int ip = 0; ip < sir->
GetNPoints(); ip++)
978 for (
int dof = 0; dof < fe->
GetDof(); dof++)
980 dshape.
GetRow(dof, gradi);
981 gradi *= LevelSet(dofs[dof]);
984 normal *= (-1. / normal.
Norml2());
1015 intp.
weight = ElemWeights(ip);
1052 bool element_int =
false;
1054 bool interior =
true;
1059 for (
int face = 0; face < me->
GetNFaces(); face++)
1110 if (face == 0 || face == 5)
1114 if (face == 1 || face == 3)
1118 if (face == 2 || face == 4)
1122 if (face == 0 || face == 1 || face == 4)
1136 for (
int dof = 0; dof <
nBasis; dof++)
1140 RHS(dof) -= (grad * normal) *
FaceWeights(ip, faces[face]);
1146 if (element_int && !interior)
1168 for (
int dof = 0; dof < fe->
GetDof(); dof++)
1170 dshape.
GetRow(dof, gradi);
1171 gradi *= LevelSet(dofs[dof]);
1174 normal *= (-1. / normal.
Norml2());
1179 for (
int dof = 0; dof <
nBasis; dof++)
1182 shapes.
GetRow(dof, grad);
1183 Mat(dof, ip) = (grad * normal);
1194 for (
int i = 0; i <
nBasis; i++)
1208 intp.
weight = ElemWeights(ip);
1234 bool element_int =
false;
1236 bool interior =
true;
1241 for (
int face = 0; face < me->
GetNFaces(); face++)
1293 if (face == 0 || face == 5)
1297 if (face == 1 || face == 3)
1301 if (face == 2 || face == 4)
1305 if (face == 0 || face == 1 || face == 4)
1323 RHS(dof) += (adiv * normal) *
FaceWeights(ip, faces[face]);
1330 if (element_int && !interior)
1346 for (
int ip = 0; ip < sir->
GetNPoints(); ip++)
1352 for (
int dof = 0; dof < fe->
GetDof(); dof++)
1354 dshape.
GetRow(dof, gradi);
1355 gradi *= LevelSet(dofs[dof]);
1358 normal *= (-1. / normal.
Norml2());
1366 shapes.
GetRow(dof, adiv);
1388 intp.
weight = ElemWeights(ip);
1413 X(0) = -1. + 2. * ip.
x;
1414 X(1) = -1. + 2. * ip.
y;
1416 for (
int c = 0; c <=
Order; c++)
1420 a(1) = pow(X(0), (
real_t)(c));
1424 b(0) = pow(X(1), (
real_t)(c));
1431 int count = 2 *
Order+ 2;
1432 for (
int c = 1; c <=
Order; c++)
1434 const int* factorial = poly.
Binom(c);
1435 for (
int expo = c; expo > 0; expo--)
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));
1467 for (
int i = 0; i <
nBasis; i++)
1468 for (
int j = 0; j < 2; j++)
1470 shapeMFN(i, j,
p) = shapeN(i, j);
1475 for (
int count = 1; count <
nBasis; count++)
1477 mGSStep(shape, shapeMFN, count);
1485 X(0) = -1. + 2. * ip.
x;
1486 X(1) = -1. + 2. * ip.
y;
1487 X(2) = -1. + 2. * ip.
z;
1497 for (
int count = step; count < shape.
Height(); count++)
1502 for (
int ip = 0; ip < ir_->
GetNPoints(); ip++)
1507 shapeMFN(ip).GetRow(count,
u);
1508 shapeMFN(ip).GetRow(step - 1, v);
1514 real_t coeff = num / den;
1518 shape.
GetRow(step - 1, s);
1524 for (
int ip = 0; ip < ir_->
GetNPoints(); ip++)
1526 shapeMFN(ip).GetRow(step - 1, s);
1527 shapeMFN(ip).GetRow(count, t);
1530 shapeMFN(ip).SetRow(count, t);
1540 X(0) = -1. + 2. * ip.
x;
1541 X(1) = -1. + 2. * ip.
y;
1544 for (
int c = 0; c <=
Order; c++)
1546 for (
int expo = 0; expo <= c; expo++)
1548 shape(count) = pow(X(0), (
real_t)(expo))
1549 * pow(X(1), (
real_t)(c - expo));
1561 X(0) = -1. + 2. * ip.
x;
1562 X(1) = -1. + 2. * ip.
y;
1565 for (
int c = 0; c <=
Order; c++)
1567 for (
int expo = 0; expo <= c; expo++)
1569 shape(count, 0) = .25 * pow(X(0), (
real_t)(expo + 1))
1570 * pow(X(1), (
real_t)(c - expo))
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);
1585 X(0) = -1. + 2. * ip.
x;
1586 X(1) = -1. + 2. * ip.
y;
1587 X(2) = -1. + 2. * ip.
z;
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++)
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));
1607 X(0) = -1. + 2. * ip.
x;
1608 X(1) = -1. + 2. * ip.
y;
1609 X(2) = -1. + 2. * ip.
z;
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++)
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));;
1715 else if (sir == NULL)
1727 else { SIR = *sir; }
1767 bool computeweights =
false;
1768 for (
int ip = 0; ip < sir.
GetNPoints(); ip++)
1772 computeweights =
true;
1796 for (
int ip = 0; ip < sir.
GetNPoints(); ip++)
1804 for (
int dof = 0; dof < fe->
GetDof(); dof++)
1806 dshape.
GetRow(dof, gradi);
1807 gradi *= LevelSet(dofs[dof]);
1811 normal *= (-1. / normal.
Norml2());
1813 weights(ip) = normphys / normref;
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.
int Size() const
Return the logical size of the array.
void DeleteAll()
Delete the whole array.
int Append(const T &el)
Append element 'el' to array, resize if necessary.
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.
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.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void SetRow(int r, const real_t *row)
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
void GetRow(int r, Vector &row) const
Rank 3 tensor (array of matrices)
Abstract data type element.
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
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
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...
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Abstract class for all finite elements.
int GetDim() const
Returns the reference space dimension for the finite element.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
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.
Class for grid function - Vector with associated FE space.
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.
Class for integration point with weight.
void Set(const real_t *p, const int dim)
void Set2w(const real_t x1, const real_t x2, const real_t w)
Class for an integration rule - an Array of IntegrationPoint.
int GetOrder() const
Returns the order of the integration rule.
int GetNPoints() const
Returns the number of the points in the integration rule.
void SetOrder(const int order)
Sets the order of the integration rule. This is only for keeping order information,...
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Container class for integration rules.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
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.
const Element * GetElement(int i) const
Return pointer to the i'th element object.
int GetNFaces() const
Return the number of faces in a 3D mesh.
int AddVertex(real_t x, real_t y=0.0, real_t z=0.0)
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...
void FinalizeQuadMesh(int generate_edges=0, int refine=0, bool fix_orientation=true)
Finalize the construction of a quadrilateral Mesh.
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
int SpaceDimension() const
Dimension of the physical space containing the mesh.
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
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().
Class for computing 1D special polynomials and their associated basis functions.
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "pchoose k" for k=0,...
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
real_t Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
void SetSize(int s)
Resize the vector to size s.
real_t lvlset(const Vector &X)
Level-set function defining the implicit interface.
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)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void subtract(const Vector &x, const Vector &y, Vector &z)
double bisect(ElementTransformation &Tr, Coefficient *LvlSet)
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
real_t p(const Vector &x, real_t t)