MFEM  v4.5.2
Finite element discretization library
marking.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 "marking.hpp"
13 
14 namespace mfem
15 {
16 
18  Array<int> &elem_marker)
19 {
20  elem_marker.SetSize(pmesh.GetNE() + pmesh.GetNSharedFaces());
21  if (!initial_marking_done) { elem_marker = SBElementType::INSIDE; }
22  else { level_set_index += 1; }
23 
25 
26  // This tolerance is relevant for points that are exactly on the zero LS.
27  const double eps = 1e-10;
28  auto outside_of_domain = [&](double value)
29  {
30  if (include_cut_cell)
31  {
32  // Points on the zero LS are considered outside the domain.
33  return (value - eps < 0.0);
34  }
35  else
36  {
37  // Points on the zero LS are considered inside the domain.
38  return (value + eps < 0.0);
39  }
40  };
41 
42  Vector vals;
43  // Check elements on the current MPI rank
44  for (int i = 0; i < pmesh.GetNE(); i++)
45  {
46  const IntegrationRule &ir = pfes_sltn->GetFE(i)->GetNodes();
47  ls_func.GetValues(i, ir, vals);
48 
49  int count = 0;
50  for (int j = 0; j < ir.GetNPoints(); j++)
51  {
52  if (outside_of_domain(vals(j))) { count++; }
53  }
54 
55  if (count == ir.GetNPoints()) // completely outside
56  {
57  elem_marker[i] = SBElementType::OUTSIDE;
58  }
59  else if (count > 0) // partially outside
60  {
61  MFEM_VERIFY(elem_marker[i] <= SBElementType::OUTSIDE,
62  " One element cut by multiple level-sets.");
63  elem_marker[i] = SBElementType::CUT + level_set_index;
64  }
65  }
66 
67  // Check neighbors on the adjacent MPI rank
68  for (int i = pmesh.GetNE(); i < pmesh.GetNE()+pmesh.GetNSharedFaces(); i++)
69  {
70  int shared_fnum = i-pmesh.GetNE();
73  int Elem2NbrNo = tr->Elem2No - pmesh.GetNE();
74 
75  ElementTransformation *eltr =
77  const IntegrationRule &ir =
78  IntRulesLo.Get(pmesh.GetElementBaseGeometry(0), 4*eltr->OrderJ());
79 
80  const int nip = ir.GetNPoints();
81  vals.SetSize(nip);
82  int count = 0;
83  for (int j = 0; j < nip; j++)
84  {
85  const IntegrationPoint &ip = ir.IntPoint(j);
86  vals(j) = ls_func.GetValue(tr->Elem2No, ip);
87  if (outside_of_domain(vals(j))) { count++; }
88  }
89 
90  if (count == ir.GetNPoints()) // completely outside
91  {
92  MFEM_VERIFY(elem_marker[i] != SBElementType::OUTSIDE,
93  "An element cannot be excluded by more than 1 level-set.");
94  elem_marker[i] = SBElementType::OUTSIDE;
95  }
96  else if (count > 0) // partially outside
97  {
98  MFEM_VERIFY(elem_marker[i] <= SBElementType::OUTSIDE,
99  "An element cannot be cut by multiple level-sets.");
100  elem_marker[i] = SBElementType::CUT + level_set_index;
101  }
102  }
103  initial_marking_done = true;
104 }
105 
107  Array<int> &sface_dof_list) const
108 {
109  if (func_dof_marking)
110  {
111  ListShiftedFaceDofs2(elem_marker, sface_dof_list);
112  return;
113  }
114 
115  sface_dof_list.DeleteAll();
116  Array<int> dofs; // work array
117 
118  // First we check interior faces of the mesh (excluding interior faces that
119  // are on the processor boundaries)
120  for (int f = 0; f < pmesh.GetNumFaces(); f++)
121  {
123  if (tr != NULL)
124  {
125  int te1 = elem_marker[tr->Elem1No], te2 = elem_marker[tr->Elem2No];
126  if (!include_cut_cell &&
128  {
129  pfes_sltn->GetFaceDofs(f, dofs);
130  sface_dof_list.Append(dofs);
131  }
132  if (!include_cut_cell &&
134  {
135  pfes_sltn->GetFaceDofs(f, dofs);
136  sface_dof_list.Append(dofs);
137  }
138  if (include_cut_cell &&
139  te1 >= SBElementType::CUT && te2 == SBElementType::OUTSIDE)
140  {
141  pfes_sltn->GetFaceDofs(f, dofs);
142  sface_dof_list.Append(dofs);
143  }
144  if (include_cut_cell &&
145  te1 == SBElementType::OUTSIDE && te2 >= SBElementType::CUT)
146  {
147  pfes_sltn->GetFaceDofs(f, dofs);
148  sface_dof_list.Append(dofs);
149  }
150  }
151  }
152 
153  // Add boundary faces that we want to model as SBM faces.
154  if (include_cut_cell)
155  {
156  for (int i = 0; i < pmesh.GetNBE(); i++)
157  {
159  if (tr != NULL)
160  {
161  if (elem_marker[tr->Elem1No] >= SBElementType::CUT)
162  {
164  sface_dof_list.Append(dofs);
165  }
166  }
167  }
168  }
169 
170  // Now we add interior faces that are on processor boundaries.
171  for (int i = 0; i < pmesh.GetNSharedFaces(); i++)
172  {
174  if (tr != NULL)
175  {
176  int ne1 = tr->Elem1No;
177  int te1 = elem_marker[ne1];
178  int te2 = elem_marker[i+pmesh.GetNE()];
179  const int faceno = pmesh.GetSharedFace(i);
180  // Add if the element on this MPI rank is completely inside the domain
181  // and the element on other MPI rank is not.
182  if (!include_cut_cell &&
184  {
185  pfes_sltn->GetFaceDofs(faceno, dofs);
186  sface_dof_list.Append(dofs);
187  }
188  if (include_cut_cell &&
189  te2 == SBElementType::OUTSIDE && te1 >= SBElementType::CUT)
190  {
191  pfes_sltn->GetFaceDofs(faceno, dofs);
192  sface_dof_list.Append(dofs);
193  }
194  }
195  }
196 }
197 
198 // Determine the list of true (i.e. conforming) essential boundary dofs. To do
199 // this, we first make a list of all dofs that are on the real boundary of the
200 // mesh, then add all the dofs of the elements that are completely outside or
201 // intersect shifted boundary. Then we remove the dofs from SBM faces
203  const Array<int> &sface_dof_list,
204  Array<int> &ess_tdof_list,
205  Array<int> &ess_shift_bdr) const
206 {
207  // Change the attribute of SBM faces that are on the boundary.
208  Array<int> ess_bdr(pmesh.bdr_attributes.Max());
209  int pmesh_bdr_attr_max = 0;
210  if (ess_bdr.Size())
211  {
212  pmesh_bdr_attr_max = pmesh.bdr_attributes.Max();
213  ess_bdr = 1;
214  }
215  bool sbm_at_true_boundary = false;
216  if (include_cut_cell)
217  {
218  for (int i = 0; i < pmesh.GetNBE(); i++)
219  {
221  if (tr != NULL)
222  {
223  if (elem_marker[tr->Elem1No] >= SBElementType::CUT)
224  {
225  pmesh.SetBdrAttribute(i, pmesh_bdr_attr_max+1);
226  sbm_at_true_boundary = true;
227  }
228  }
229  }
230  }
231  bool sbm_at_true_boundary_global = false;
232  MPI_Allreduce(&sbm_at_true_boundary, &sbm_at_true_boundary_global, 1,
233  MPI_C_BOOL, MPI_LOR, MPI_COMM_WORLD);
234  if (sbm_at_true_boundary_global)
235  {
236  ess_bdr.Append(0);
238  }
239 
240  // Make a list of dofs on all boundaries
241  ess_shift_bdr.SetSize(ess_bdr.Size());
242  if (pmesh.bdr_attributes.Size())
243  {
244  for (int i = 0; i < ess_bdr.Size(); i++)
245  {
246  ess_shift_bdr[i] = 1 - ess_bdr[i];
247  }
248  }
249  Array<int> ess_vdofs_bdr;
250  pfes_sltn->GetEssentialVDofs(ess_bdr, ess_vdofs_bdr);
251 
252  // Get all dofs associated with elements outside the domain or intersected by
253  // the boundary.
254  Array<int> ess_vdofs(ess_vdofs_bdr.Size()), dofs;
255  ess_vdofs = 0;
256  for (int e = 0; e < pmesh.GetNE(); e++)
257  {
258  if (!include_cut_cell &&
259  (elem_marker[e] == SBElementType::OUTSIDE ||
260  elem_marker[e] >= SBElementType::CUT))
261  {
262  pfes_sltn->GetElementVDofs(e, dofs);
263  for (int i = 0; i < dofs.Size(); i++)
264  {
265  ess_vdofs[dofs[i]] = -1;
266  }
267  }
268  if (include_cut_cell &&
269  elem_marker[e] == SBElementType::OUTSIDE)
270  {
271  pfes_sltn->GetElementVDofs(e, dofs);
272  for (int i = 0; i < dofs.Size(); i++)
273  {
274  ess_vdofs[dofs[i]] = -1;
275  }
276  }
277  }
278 
279  // Combine the lists to mark essential dofs.
280  for (int i = 0; i < ess_vdofs.Size(); i++)
281  {
282  if (ess_vdofs_bdr[i] == -1) { ess_vdofs[i] = -1; }
283  }
284 
285  // Unmark dofs that are on SBM faces (but not on Dirichlet boundaries)
286  for (int i = 0; i < sface_dof_list.Size(); i++)
287  {
288  if (ess_vdofs_bdr[sface_dof_list[i]] != -1)
289  {
290  ess_vdofs[sface_dof_list[i]] = 0;
291  }
292  }
293 
294  // Synchronize
295  for (int i = 0; i < ess_vdofs.Size() ; i++) { ess_vdofs[i] += 1; }
296  pfes_sltn->Synchronize(ess_vdofs);
297  for (int i = 0; i < ess_vdofs.Size() ; i++) { ess_vdofs[i] -= 1; }
298 
299  // Convert to tdofs
300  Array<int> ess_tdofs;
301  pfes_sltn->GetRestrictionMatrix()->BooleanMult(ess_vdofs, ess_tdofs);
302  pfes_sltn->MarkerToList(ess_tdofs, ess_tdof_list);
303 }
304 
306  Array<int> &sface_dof_list) const
307 {
308  sface_dof_list.DeleteAll();
309 
310  L2_FECollection mat_coll(0, pmesh.Dimension());
311  ParFiniteElementSpace mat_fes(&pmesh, &mat_coll);
312  ParGridFunction mat(&mat_fes);
313  ParGridFunction marker_gf(pfes_sltn);
314  for (int i = 0; i < pmesh.GetNE(); i++)
315  {
316  // 0 is inside, 1 is outside.
317  mat(i) = 0.0;
318  if (elem_marker[i] == SBElementType::OUTSIDE)
319  { mat(i) = 1.0; }
320  if (elem_marker[i] >= SBElementType::CUT && include_cut_cell == false)
321  { mat(i) = 1.0; }
322  }
323 
324  GridFunctionCoefficient coeff_mat(&mat);
325  marker_gf.ProjectDiscCoefficient(coeff_mat, GridFunction::ARITHMETIC);
326 
327  for (int j = 0; j < marker_gf.Size(); j++)
328  {
329  if (marker_gf(j) > 0.1 && marker_gf(j) < 0.9)
330  {
331  sface_dof_list.Append(j);
332  }
333  }
334 
335  // Add boundary faces that we want to model as SBM faces.
336  if (include_cut_cell)
337  {
338  Array<int> dofs;
339  for (int i = 0; i < pmesh.GetNBE(); i++)
340  {
342  if (tr != NULL)
343  {
344  if (elem_marker[tr->Elem1No] >= SBElementType::CUT)
345  {
347  sface_dof_list.Append(dofs);
348  }
349  }
350  }
351  }
352 }
353 
354 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
IntegrationRules IntRulesLo(0, Quadrature1D::GaussLobatto)
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1108
int Dimension() const
Definition: mesh.hpp:1047
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5389
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
const bool include_cut_cell
Definition: marking.hpp:29
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition: array.cpp:68
Container class for integration rules.
Definition: intrules.hpp:311
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:546
Abstract parallel finite element space.
Definition: pfespace.hpp:28
FaceElementTransformations * GetInteriorFaceTransformations(int FaceNo)
See GetFaceElementTransformations().
Definition: mesh.hpp:1327
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:389
virtual int OrderJ() const =0
Return the order of the elements of the Jacobian of the transformation.
int GetNBE() const
Returns number of boundary elements.
Definition: mesh.hpp:939
void DeleteAll()
Delete the whole array.
Definition: array.hpp:851
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:598
void ListShiftedFaceDofs2(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
Definition: marking.cpp:305
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:1021
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:524
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3080
void Synchronize(Array< int > &ldof_marker) const
Given an integer array on the local degrees of freedom, perform a bitwise OR between the shared dofs...
Definition: pfespace.cpp:1009
ParFiniteElementSpace * pfes_sltn
Definition: marking.hpp:26
void ListEssentialTDofs(const Array< int > &elem_marker, const Array< int > &sface_dof_list, Array< int > &ess_tdof_list, Array< int > &ess_shift_bdr) const
Definition: marking.cpp:202
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1095
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3215
int GetBdrFace(int BdrElemNo) const
Return the local face index for the given boundary face.
Definition: mesh.cpp:1113
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
virtual int GetFaceDofs(int i, Array< int > &dofs, int variant=0) const
Definition: pfespace.cpp:519
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:275
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3196
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:535
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
void SetBdrAttribute(int i, int attr)
Set the attribute of boundary element i.
Definition: mesh.hpp:1553
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:264
void SetAttributes() override
Definition: pmesh.cpp:1580
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
void GetFaceNbrElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: pmesh.cpp:1976
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:936
Class for integration point with weight.
Definition: intrules.hpp:25
void ListShiftedFaceDofs(const Array< int > &elem_marker, Array< int > &sface_dof_list) const
Definition: marking.cpp:106
const bool func_dof_marking
Definition: marking.hpp:34
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
Vector data type.
Definition: vector.hpp:60
Class for parallel grid function.
Definition: pgridfunc.hpp:32
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:388
void MarkElements(const ParGridFunction &ls_func, Array< int > &elem_marker)
Definition: marking.cpp:17
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition: fe_coll.hpp:320