16#include "moonolith_build_quadrature.hpp"
18using namespace mfem::internal;
25template <
class Polytope>
class CutGeneric :
public Cut
28 using Point =
typename Polytope::Point;
29 static const int Dim = Point::n_dims;
30 using Quadrature_t = moonolith::Quadrature<double, Dim>;
31 using BuildQuadrature_t = moonolith::BuildQuadrature<Polytope>;
34 const int from_elem_idx,
39 void Describe()
const override;
42 inline Quadrature_t &GetQRule() {
return q_rule_; }
43 inline int GetOrder()
const {
return order_; }
45 inline void SetOrder(
const int order) { order_ = order; }
47 virtual void MakePolytope(
Mesh &mesh,
const int elem_idx,
48 Polytope &polygon) = 0;
50 bool IsValidPhysicalPoint(
const Vector &p_mfem)
const
54 for(
int d = 0; d < Dim; ++d) {
58 moonolith::HPolytope<double, Dim> poly;
61 if(!poly.contains(
p, 1e-8))
return false;
65 return poly.contains(
p, 1e-8);
78 assert(
p.x <= 1 + 1e-8);
79 assert(
p.y <= 1 + 1e-8);
80 assert(
p.z <= 1 + 1e-8);
84 ok = ok && (
p.x ==
p.x);
85 ok = ok && (
p.y ==
p.y);
86 ok = ok && (
p.z ==
p.z);
88 ok = ok && (
p.x >= -1e-8);
89 ok = ok && (
p.y >= -1e-8);
90 ok = ok && (
p.z >= -1e-8);
92 ok = ok && (
p.x <= 1 + 1e-8);
93 ok = ok && (
p.y <= 1 + 1e-8);
94 ok = ok && (
p.z <= 1 + 1e-8);
98 static void ConvertQRule(
const IntegrationRule &ir, Quadrature_t &result)
100 const int size = ir.
Size();
104 for (
int k = 0; k < size; ++k)
107 result.points[k][0] = qp.x;
108 result.points[k][1] = qp.y;
109 if constexpr (Dim == 3) {
110 result.points[k][12] = qp.z;
113 result.weights[k] = qp.weight;
117 inline const Polytope & from()
const {
return from_; }
118 inline const Polytope & to()
const {
return to_; }
122 Quadrature_t q_rule_;
123 Quadrature_t physical_quadrature_;
124 BuildQuadrature_t builder_;
126 double intersection_measure_{0};
129class Cut2D :
public CutGeneric<moonolith::Polygon<double, 2>>
132 using Polygon_t = moonolith::Polygon<double, 2>;
133 using Quadrature_t = moonolith::Quadrature<double, 2>;
134 using BuildQuadrature_t = moonolith::BuildQuadrature<Polygon_t>;
136 void SetIntegrationOrder(
const int order)
override;
140 void MakePolytope(
Mesh &mesh,
const int elem_idx,
141 Polygon_t &polygon)
override;
147class Cut3D :
public CutGeneric<moonolith::Polyhedron<double>>
150 using Polyhedron_t = moonolith::Polyhedron<double>;
151 using Quadrature_t = moonolith::Quadrature<double, 3>;
152 using BuildQuadrature_t = moonolith::BuildQuadrature<Polyhedron_t>;
154 void SetIntegrationOrder(
const int order)
override;
158 void MakePolytope(
Mesh &mesh,
const int elem_idx,
159 Polyhedron_t &polyhedron)
override;
168 const Vector &physical_p,
const double &w,
175 assert(ref_p.
x >= -1e-8);
176 assert(ref_p.
y >= -1e-8);
177 assert(ref_p.
z >= -1e-8);
179 assert(ref_p.
x <= 1 + 1e-8);
180 assert(ref_p.
y <= 1 + 1e-8);
181 assert(ref_p.
z <= 1 + 1e-8);
183#ifdef MFEM_DEBUG_MOONOLITH
188 physical_test_p -= physical_p;
190 assert(physical_test_p.
Norml2() < 1e-10);
206 assert(
dim == 2 ||
dim == 3);
210template <
class Polytope>
212 const int from_elem_idx,
214 const int to_elem_idx,
219 MakePolytope(*from_space.
GetMesh(), from_elem_idx, from_);
220 MakePolytope(*to_space.
GetMesh(), to_elem_idx, to_);
222 if (!builder_.apply(q_rule_, from_, to_, physical_quadrature_))
228#ifdef MFEM_DEBUG_MOONOLITH
231 static int counter = 0;
233 moonolith::MatlabScripter script;
235 script.plot(from_,
"\'g-\'");
236 script.plot(from_,
"\'g.\'");
237 script.plot(to_,
"\'r-\'");
238 script.plot(to_,
"\'r.\'");
239 script.plot(physical_quadrature_.points,
"\'*b\'");
241 std::cout <<
"measure(" << counter <<
"): " << moonolith::measure(physical_quadrature_) <<
"\n";
243 script.save(
"out_" + std::to_string(counter++) +
".m");
251 const int n_qp = physical_quadrature_.n_points();
261 double from_measure = moonolith::measure(from_);
262 double to_measure = moonolith::measure(to_);
265 for (
int qp = 0; qp < n_qp; ++qp)
268 for (
int d = 0; d < Dim; ++d)
270 p(d) = physical_quadrature_.points[qp][d];
273 double w = physical_quadrature_.weights[qp];
274 intersection_measure_ += w;
276 assert(IsValidPhysicalPoint(
p));
279 from_quadrature[qp]);
285 assert(IsValidQPoint(from_quadrature[qp]));
286 assert(IsValidQPoint(to_quadrature[qp]));
292template <
class Polytope>
void CutGeneric<Polytope>::Describe()
const
294 mfem::out <<
"Cut measure " << intersection_measure_ <<
'\n';
297template class CutGeneric<::moonolith::Polygon<double, 2>>;
298template class CutGeneric<::moonolith::Polyhedron<double>>;
300void Cut2D::MakePolytope(
Mesh &mesh,
const int elem_idx, Polygon_t &polygon)
305 const int n_points = buffer_pts.Width();
306 polygon.resize(n_points);
308 for (
int k = 0; k < n_points; ++k)
310 for (
int d = 0; d < 2; ++d)
312 polygon.points[k][d] = buffer_pts(d, k);
316 assert(polygon.check_convexity());
317 assert(::moonolith::measure(polygon) > 0.0);
324 const int size = ir.
Size();
325 auto &q_rule = this->GetQRule();
330 for (
int k = 0; k < size; ++k)
333 q_rule.points[k][0] = qp.x;
334 q_rule.points[k][1] = qp.y;
335 q_rule.weights[k] = qp.weight;
343void Cut2D::SetIntegrationOrder(
const int order)
345 if (this->GetOrder() != order)
349 SetQuadratureRule(ir);
353void Cut3D::SetIntegrationOrder(
const int order)
355 if (this->GetOrder() != order)
359 SetQuadratureRule(ir);
363void Cut3D::MakePolytope(
Mesh &mesh,
const int elem_idx,
364 Polyhedron_t &polyhedron)
372 const int e_type = e.
GetType();
378 const int n_faces = buffer_faces.Size();
380 polyhedron.el_ptr.resize(n_faces + 1);
381 polyhedron.points.resize(buffer_vertices.Size());
382 polyhedron.el_index.resize(MaxVertsXFace(e_type) * n_faces);
383 polyhedron.el_ptr[0] = 0;
385 for (
int i = 0; i < buffer_vertices.Size(); ++i)
387 for (
int j = 0; j <
dim; ++j)
389 polyhedron.points[i][j] = buffer_pts(j, i);
394 for (
int i = 0; i < buffer_faces.Size(); ++i)
397 const int eptr = polyhedron.el_ptr[i];
399 for (
int j = 0; j < f2v.
Size(); ++j)
401 const int v_offset = buffer_vertices.Find(f2v[j]);
402 polyhedron.el_index[eptr + j] = v_offset;
405 polyhedron.el_ptr[i + 1] = polyhedron.el_ptr[i] + f2v.
Size();
408 polyhedron.fix_ordering();
414 const int size = ir.
Size();
416 auto &q_rule = this->GetQRule();
422 for (
int k = 0; k < size; ++k)
425 q_rule.points[k][0] = qp.x;
426 q_rule.points[k][1] = qp.y;
427 q_rule.points[k][2] = qp.z;
428 q_rule.weights[k] = qp.weight;
440 return std::make_shared<Cut2D>();
444 return std::make_shared<Cut3D>();
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
All subclasses of Cut will implement intersection routines and quadrature point generation within the...
Data type dense matrix using column-major storage.
Abstract data type element.
virtual Type GetType() const =0
Returns element's type.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th element in t...
Mesh * GetMesh() const
Returns the mesh.
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetOrder() const
Returns the order of the integration rule.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
const Element * GetElement(int i) const
Return pointer to the i'th element object.
int Dimension() const
Dimension of the reference space used within the elements.
void GetElementFaces(int i, Array< int > &faces, Array< int > &ori) const
Return the indices and the orientations of all faces of element i.
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
void GetPointMatrix(int i, DenseMatrix &pointmat) const
real_t Norml2() const
Returns the l2 norm of the vector.
int Size() const
Returns the size of the vector.
std::shared_ptr< Cut > NewCut(const int dim)
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void TransformToReference(ElementTransformation &Trans, int type, const Vector &physical_p, const double &w, IntegrationPoint &ref_p)
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)