MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
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#include "cut.hpp"
13
14#include "transferutils.hpp"
15
16#include "moonolith_build_quadrature.hpp"
17
18using namespace mfem::internal;
19
20// #define MFEM_DEBUG_MOONOLITH
21
22namespace mfem
23{
24
25template <class Polytope> class CutGeneric : public Cut
26{
27public:
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>;
32
33 bool BuildQuadrature(const FiniteElementSpace &from_space,
34 const int from_elem_idx,
35 const FiniteElementSpace &to_space,
36 const int to_elem_idx, IntegrationRule &from_quadrature,
37 IntegrationRule &to_quadrature) override;
38
39 void Describe() const override;
40
41protected:
42 inline Quadrature_t &GetQRule() { return q_rule_; }
43 inline int GetOrder() const { return order_; }
44
45 inline void SetOrder(const int order) { order_ = order; }
46
47 virtual void MakePolytope(Mesh &mesh, const int elem_idx,
48 Polytope &polygon) = 0;
49
50 bool IsValidPhysicalPoint(const Vector &p_mfem) const
51 {
52 Point p;
53
54 for(int d = 0; d < Dim; ++d) {
55 p[d] = p_mfem[d];
56 }
57
58 moonolith::HPolytope<double, Dim> poly;
59 poly.make(from());
60
61 if(!poly.contains(p, 1e-8)) return false;
62
63 poly.make(to());
64
65 return poly.contains(p, 1e-8);
66 }
67
68 bool IsValidQPoint(const IntegrationPoint &p) const
69 {
70 assert(p.x == p.x);
71 assert(p.y == p.y);
72 assert(p.z == p.z);
73
74 assert(p.x >= -1e-8);
75 assert(p.y >= -1e-8);
76 assert(p.z >= -1e-8);
77
78 assert(p.x <= 1 + 1e-8);
79 assert(p.y <= 1 + 1e-8);
80 assert(p.z <= 1 + 1e-8);
81
82 bool ok = true;
83
84 ok = ok && (p.x == p.x);
85 ok = ok && (p.y == p.y);
86 ok = ok && (p.z == p.z);
87
88 ok = ok && (p.x >= -1e-8);
89 ok = ok && (p.y >= -1e-8);
90 ok = ok && (p.z >= -1e-8);
91
92 ok = ok && (p.x <= 1 + 1e-8);
93 ok = ok && (p.y <= 1 + 1e-8);
94 ok = ok && (p.z <= 1 + 1e-8);
95 return ok;
96 }
97
98 static void ConvertQRule(const IntegrationRule &ir, Quadrature_t &result)
99 {
100 const int size = ir.Size();
101
102 result.resize(size);
103
104 for (int k = 0; k < size; ++k)
105 {
106 auto &qp = ir[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;
111 }
112
113 result.weights[k] = qp.weight;
114 }
115 }
116
117 inline const Polytope & from() const { return from_; }
118 inline const Polytope & to() const { return to_; }
119
120private:
121 Polytope from_, to_;
122 Quadrature_t q_rule_;
123 Quadrature_t physical_quadrature_;
124 BuildQuadrature_t builder_;
125 int order_{-1};
126 double intersection_measure_{0};
127};
128
129class Cut2D : public CutGeneric<moonolith::Polygon<double, 2>>
130{
131public:
132 using Polygon_t = moonolith::Polygon<double, 2>;
133 using Quadrature_t = moonolith::Quadrature<double, 2>;
134 using BuildQuadrature_t = moonolith::BuildQuadrature<Polygon_t>;
135
136 void SetIntegrationOrder(const int order) override;
137
138protected:
139 void SetQuadratureRule(const IntegrationRule &ir) override;
140 void MakePolytope(Mesh &mesh, const int elem_idx,
141 Polygon_t &polygon) override;
142
143private:
144 DenseMatrix buffer_pts;
145};
146
147class Cut3D : public CutGeneric<moonolith::Polyhedron<double>>
148{
149public:
150 using Polyhedron_t = moonolith::Polyhedron<double>;
151 using Quadrature_t = moonolith::Quadrature<double, 3>;
152 using BuildQuadrature_t = moonolith::BuildQuadrature<Polyhedron_t>;
153
154 void SetIntegrationOrder(const int order) override;
155
156protected:
157 void SetQuadratureRule(const IntegrationRule &ir) override;
158 void MakePolytope(Mesh &mesh, const int elem_idx,
159 Polyhedron_t &polyhedron) override;
160
161private:
162 DenseMatrix buffer_pts;
163 Array<int> buffer_vertices;
164 Array<int> buffer_faces, buffer_cor;
165};
166
168 const Vector &physical_p, const double &w,
169 IntegrationPoint &ref_p)
170{
171
172 int dim = physical_p.Size();
173 Trans.TransformBack(physical_p, ref_p);
174
175 assert(ref_p.x >= -1e-8);
176 assert(ref_p.y >= -1e-8);
177 assert(ref_p.z >= -1e-8);
178
179 assert(ref_p.x <= 1 + 1e-8);
180 assert(ref_p.y <= 1 + 1e-8);
181 assert(ref_p.z <= 1 + 1e-8);
182
183#ifdef MFEM_DEBUG_MOONOLITH
184 {
185 Vector physical_test_p;
186 Trans.Transform(ref_p, physical_test_p);
187
188 physical_test_p -= physical_p;
189
190 assert(physical_test_p.Norml2() < 1e-10);
191 }
192#endif
193
194 ref_p.weight = w;
195
196 if (type == Geometry::TRIANGLE && dim == 2)
197 {
198 ref_p.weight *= 0.5;
199 }
200 else if (type == Geometry::TETRAHEDRON && dim == 3)
201 {
202 ref_p.weight *= 1. / 6;
203 }
204 else
205 {
206 assert(dim == 2 || dim == 3);
207 }
208}
209
210template <class Polytope>
211bool CutGeneric<Polytope>::BuildQuadrature(const FiniteElementSpace &from_space,
212 const int from_elem_idx,
213 const FiniteElementSpace &to_space,
214 const int to_elem_idx,
215 IntegrationRule &from_quadrature,
216 IntegrationRule &to_quadrature)
217{
218
219 MakePolytope(*from_space.GetMesh(), from_elem_idx, from_);
220 MakePolytope(*to_space.GetMesh(), to_elem_idx, to_);
221
222 if (!builder_.apply(q_rule_, from_, to_, physical_quadrature_))
223 {
224 return false;
225 }
226
227
228#ifdef MFEM_DEBUG_MOONOLITH
229
230 {
231 static int counter = 0;
232
233 moonolith::MatlabScripter script;
234 script.hold_on();
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\'");
240
241 std::cout << "measure(" << counter << "): " << moonolith::measure(physical_quadrature_) << "\n";
242
243 script.save("out_" + std::to_string(counter++) + ".m");
244 }
245
246#endif // MFEM_DEBUG_MOONOLITH
247
248 int from_type = from_space.GetFE(from_elem_idx)->GetGeomType();
249 int to_type = to_space.GetFE(to_elem_idx)->GetGeomType();
250
251 const int n_qp = physical_quadrature_.n_points();
252
253 from_quadrature.SetSize(n_qp);
254 to_quadrature.SetSize(n_qp);
255
256 ElementTransformation &from_trans =
257 *from_space.GetElementTransformation(from_elem_idx);
258 ElementTransformation &to_trans =
259 *to_space.GetElementTransformation(to_elem_idx);
260
261 double from_measure = moonolith::measure(from_);
262 double to_measure = moonolith::measure(to_);
263
264 Vector p(Dim);
265 for (int qp = 0; qp < n_qp; ++qp)
266 {
267
268 for (int d = 0; d < Dim; ++d)
269 {
270 p(d) = physical_quadrature_.points[qp][d];
271 }
272
273 double w = physical_quadrature_.weights[qp];
274 intersection_measure_ += w;
275
276 assert(IsValidPhysicalPoint(p));
277
278 TransformToReference(from_trans, from_type, p, w / from_measure,
279 from_quadrature[qp]);
280
281 TransformToReference(to_trans, to_type, p, w / to_measure,
282 to_quadrature[qp]);
283
284
285 assert(IsValidQPoint(from_quadrature[qp]));
286 assert(IsValidQPoint(to_quadrature[qp]));
287 }
288
289 return true;
290}
291
292template <class Polytope> void CutGeneric<Polytope>::Describe() const
293{
294 mfem::out << "Cut measure " << intersection_measure_ << '\n';
295}
296
297template class CutGeneric<::moonolith::Polygon<double, 2>>;
298template class CutGeneric<::moonolith::Polyhedron<double>>;
299
300void Cut2D::MakePolytope(Mesh &mesh, const int elem_idx, Polygon_t &polygon)
301{
302
303 mesh.GetPointMatrix(elem_idx, buffer_pts);
304
305 const int n_points = buffer_pts.Width();
306 polygon.resize(n_points);
307
308 for (int k = 0; k < n_points; ++k)
309 {
310 for (int d = 0; d < 2; ++d)
311 {
312 polygon.points[k][d] = buffer_pts(d, k);
313 }
314 }
315
316 assert(polygon.check_convexity());
317 assert(::moonolith::measure(polygon) > 0.0);
318}
319
320
321
322void Cut2D::SetQuadratureRule(const IntegrationRule &ir)
323{
324 const int size = ir.Size();
325 auto &q_rule = this->GetQRule();
326
327 q_rule.resize(size);
328
329 double rule_w = 0.0;
330 for (int k = 0; k < size; ++k)
331 {
332 auto &qp = ir[k];
333 q_rule.points[k][0] = qp.x;
334 q_rule.points[k][1] = qp.y;
335 q_rule.weights[k] = qp.weight;
336 rule_w += qp.weight;
337 }
338
339 this->SetOrder(ir.GetOrder());
340 q_rule.normalize();
341}
342
343void Cut2D::SetIntegrationOrder(const int order)
344{
345 if (this->GetOrder() != order)
346 {
348 assert(ir.GetOrder() >= order);
349 SetQuadratureRule(ir);
350 }
351}
352
353void Cut3D::SetIntegrationOrder(const int order)
354{
355 if (this->GetOrder() != order)
356 {
358 assert(ir.GetOrder() >= order);
359 SetQuadratureRule(ir);
360 }
361}
362
363void Cut3D::MakePolytope(Mesh &mesh, const int elem_idx,
364 Polyhedron_t &polyhedron)
365{
366 using namespace std;
367 const int dim = mesh.Dimension();
368
369 assert(mesh.GetElement(elem_idx));
370
371 const Element &e = *mesh.GetElement(elem_idx);
372 const int e_type = e.GetType();
373
374 mesh.GetElementFaces(elem_idx, buffer_faces, buffer_cor);
375 mesh.GetPointMatrix(elem_idx, buffer_pts);
376 mesh.GetElementVertices(elem_idx, buffer_vertices);
377
378 const int n_faces = buffer_faces.Size();
379 polyhedron.clear();
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;
384
385 for (int i = 0; i < buffer_vertices.Size(); ++i)
386 {
387 for (int j = 0; j < dim; ++j)
388 {
389 polyhedron.points[i][j] = buffer_pts(j, i);
390 }
391 }
392
393 Array<int> f2v;
394 for (int i = 0; i < buffer_faces.Size(); ++i)
395 {
396 mesh.GetFaceVertices(buffer_faces[i], f2v);
397 const int eptr = polyhedron.el_ptr[i];
398
399 for (int j = 0; j < f2v.Size(); ++j)
400 {
401 const int v_offset = buffer_vertices.Find(f2v[j]);
402 polyhedron.el_index[eptr + j] = v_offset;
403 }
404
405 polyhedron.el_ptr[i + 1] = polyhedron.el_ptr[i] + f2v.Size();
406 }
407
408 polyhedron.fix_ordering();
409}
410
411
412void Cut3D::SetQuadratureRule(const IntegrationRule &ir)
413{
414 const int size = ir.Size();
415
416 auto &q_rule = this->GetQRule();
417
418 q_rule.resize(size);
419
420 double rule_w = 0.0;
421
422 for (int k = 0; k < size; ++k)
423 {
424 auto &qp = ir[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;
429 rule_w += qp.weight;
430 }
431
432 this->SetOrder(ir.GetOrder());
433 q_rule.normalize();
434}
435
436std::shared_ptr<Cut> NewCut(const int dim)
437{
438 if (dim == 2)
439 {
440 return std::make_shared<Cut2D>();
441 }
442 else if (dim == 3)
443 {
444 return std::make_shared<Cut3D>();
445 }
446 else
447 {
448 assert(false);
449 return nullptr;
450 }
451}
452
453} // namespace mfem
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition array.hpp:697
int Size() const
Return the logical size of the array.
Definition array.hpp:144
int Find(const T &el) const
Return the first index where 'el' is found; return -1 if not found.
Definition array.hpp:828
All subclasses of Cut will implement intersection routines and quadrature point generation within the...
Definition cut.hpp:29
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
virtual int TransformBack(const Vector &pt, IntegrationPoint &ip, const real_t phys_tol=1e-15)=0
Transform a point pt from physical space to a point ip in reference space and optionally can set a so...
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Abstract data type element.
Definition element.hpp:29
virtual Type GetType() const =0
Returns element's type.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition fespace.hpp:773
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
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:326
Class for integration point with weight.
Definition intrules.hpp:35
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
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Mesh data type.
Definition mesh.hpp:56
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1438
const Element * GetElement(int i) const
Return pointer to the i'th element object.
Definition mesh.hpp:1283
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
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
void GetFaceVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of face i.
Definition mesh.hpp:1456
void GetPointMatrix(int i, DenseMatrix &pointmat) const
Definition mesh.cpp:7326
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Data type point element.
Definition point.hpp:23
Vector data type.
Definition vector.hpp:80
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
int dim
Definition ex24.cpp:53
std::shared_ptr< Cut > NewCut(const int dim)
Definition cut.cpp:436
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
void TransformToReference(ElementTransformation &Trans, int type, const Vector &physical_p, const double &w, IntegrationPoint &ref_p)
Definition cut.cpp:167
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition intrules.hpp:486
real_t p(const Vector &x, real_t t)