MFEM  v4.5.2
Finite element discretization library
cut.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 
18 using namespace mfem::internal;
19 
20 // #define MFEM_DEBUG_MOONOLITH
21 
22 namespace mfem
23 {
24 
25 template <class Polytope> class CutGeneric : public Cut
26 {
27 public:
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 
41 protected:
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 
120 private:
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 
129 class Cut2D : public CutGeneric<moonolith::Polygon<double, 2>>
130 {
131 public:
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 
138 protected:
139  void SetQuadratureRule(const IntegrationRule &ir) override;
140  void MakePolytope(Mesh &mesh, const int elem_idx,
141  Polygon_t &polygon) override;
142 
143 private:
144  DenseMatrix buffer_pts;
145 };
146 
147 class Cut3D : public CutGeneric<moonolith::Polyhedron<double>>
148 {
149 public:
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 
156 protected:
157  void SetQuadratureRule(const IntegrationRule &ir) override;
158  void MakePolytope(Mesh &mesh, const int elem_idx,
159  Polyhedron_t &polyhedron) override;
160 
161 private:
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 
210 template <class Polytope>
211 bool 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 
292 template <class Polytope> void CutGeneric<Polytope>::Describe() const
293 {
294  mfem::out << "Cut measure " << intersection_measure_ << '\n';
295 }
296 
297 template class CutGeneric<::moonolith::Polygon<double, 2>>;
298 template class CutGeneric<::moonolith::Polyhedron<double>>;
299 
300 void 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 
322 void 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 
343 void Cut2D::SetIntegrationOrder(const int order)
344 {
345  if (this->GetOrder() != order)
346  {
347  const IntegrationRule &ir = IntRules.Get(Geometry::TRIANGLE, order);
348  assert(ir.GetOrder() >= order);
349  SetQuadratureRule(ir);
350  }
351 }
352 
353 void Cut3D::SetIntegrationOrder(const int order)
354 {
355  if (this->GetOrder() != order)
356  {
357  const IntegrationRule &ir = IntRules.Get(Geometry::TETRAHEDRON, order);
358  assert(ir.GetOrder() >= order);
359  SetQuadratureRule(ir);
360  }
361 }
362 
363 void 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 
412 void 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 
436 std::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
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
STL namespace.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
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
std::shared_ptr< Cut > NewCut(const int dim)
Definition: cut.cpp:436
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:804
Vector data type.
Definition: vector.hpp:60
void TransformToReference(ElementTransformation &Trans, int type, const Vector &physical_p, const double &w, IntegrationPoint &ref_p)
Definition: cut.cpp:167