MFEM  v4.4.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
sbm_solver.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 #ifndef MFEM_SBM_SOLVER_HPP
13 #define MFEM_SBM_SOLVER_HPP
14 
15 #include "mfem.hpp"
16 #include "marking.hpp"
17 
18 namespace mfem
19 {
20 
21 /// ShiftedFunctionCoefficient, similar to FunctionCoefficient, but also takes
22 /// into account a displacement vector if specified.
24 {
25 protected:
26  std::function<double(const Vector &)> Function;
27  double constant = 0.0;
29 
30 public:
31  ShiftedFunctionCoefficient(std::function<double(const Vector &v)> F)
32  : Function(std::move(F)), constantcoefficient(false) { }
33  ShiftedFunctionCoefficient(double constant_)
34  : constant(constant_), constantcoefficient(true) { }
35 
36  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
37  {
38  if (constantcoefficient) { return constant; }
39 
40  Vector D(T.GetSpaceDim());
41  D = 0.;
42  return (this)->Eval(T, ip, D);
43  }
44 
45  /// Evaluate the coefficient at @a ip + @a D.
46  double Eval(ElementTransformation &T,
47  const IntegrationPoint &ip,
48  const Vector &D);
49 };
50 
52 {
53 protected:
54  std::function<void(const Vector &, Vector &)> Function;
55 
56 public:
58  std::function<void(const Vector &, Vector &)> F)
59  : VectorCoefficient(dim), Function(std::move(F)) { }
60 
62  virtual void Eval(Vector &V, ElementTransformation &T,
63  const IntegrationPoint &ip)
64  {
65  Vector D(vdim);
66  D = 0.;
67  return (this)->Eval(V, T, ip, D);
68  }
69 
70  /// Evaluate the coefficient at @a ip + @a D.
71  void Eval(Vector &V,
73  const IntegrationPoint &ip,
74  const Vector &D);
75 };
76 
77 /// BilinearFormIntegrator for the high-order extension of shifted boundary
78 /// method.
79 /// A(u, w) = -<nabla u.n, w>
80 /// -<u + nabla u.d + h.o.t, nabla w.n>
81 /// +<alpha h^{-1} (u + nabla u.d + h.o.t), w + nabla w.d + h.o.t>
82 /// where h.o.t include higher-order derivatives (nabla^k u) due to Taylor
83 /// expansion. Since this interior face integrator is applied to the surrogate
84 /// boundary (see marking.hpp for notes on how the surrogate faces are
85 /// determined and elements are marked), this integrator adds contribution to
86 /// only the element that is adjacent to that face (Trans.Elem1 or Trans.Elem2)
87 /// and is part of the surrogate domain.
89 {
90 protected:
91  double alpha;
92  VectorCoefficient *vD; // Distance function coefficient
93  Array<int> *elem_marker; // marker indicating whether element is inside,
94  //cut, or outside the domain.
95  bool include_cut_cell; // include element cut by true boundary
96  int nterms; // Number of terms in addition to the gradient
97  // term from Taylor expansion that should be included. (0 by default).
98  int NEproc; // Number of elements on the current MPI rank
100  Array<int> cut_marker; // Array with marker values for cut-cell
101  // corresponding to the level set that BilinearForm applies to.
102 
103  // these are not thread-safe!
106 
107 
108 public:
110  const double a,
111  VectorCoefficient &vD_,
112  Array<int> &elem_marker_,
113  Array<int> &cut_marker_,
114  bool include_cut_cell_ = false,
115  int nterms_ = 0)
116  : alpha(a), vD(&vD_),
117  elem_marker(&elem_marker_),
118  include_cut_cell(include_cut_cell_),
119  nterms(nterms_),
120  NEproc(pmesh->GetNE()),
122  cut_marker(cut_marker_) { }
123 
125  virtual void AssembleFaceMatrix(const FiniteElement &el1,
126  const FiniteElement &el2,
128  DenseMatrix &elmat);
129 
131 };
132 
133 /// LinearFormIntegrator for the high-order extension of shifted boundary
134 /// method.
135 /// (u, w) = -<u_D, nabla w.n >
136 /// +<alpha h^{-1} u_D, w + nabla w.d + h.o.t>
137 /// where h.o.t include higher-order derivatives (nabla^k u) due to Taylor
138 /// expansion. Since this interior face integrator is applied to the surrogate
139 /// boundary (see marking.hpp for notes on how the surrogate faces are
140 /// determined and elements are marked), this integrator adds contribution to
141 /// only the element that is adjacent to that face (Trans.Elem1 or Trans.Elem2)
142 /// and is part of the surrogate domain.
143 /// Note that u_D is evaluated at the true boundary using the distance function
144 /// and ShiftedFunctionCoefficient, i.e. u_D(x_true) = u_D(x_surrogate + D),
145 /// where x_surrogate is the location of the integration point on the surrogate
146 /// boundary and D is the distance vector from the surrogate boundary to the
147 /// true boundary.
149 {
150 protected:
152  double alpha; // Nitsche parameter
153  VectorCoefficient *vD; // Distance function coefficient
154  Array<int> *elem_marker; // marker indicating whether element is inside,
155  //cut, or outside the domain.
156  bool include_cut_cell; // include element cut by true boundary
157  int nterms; // Number of terms in addition to the gradient
158  // term from Taylor expansion that should be included. (0 by default).
159  int NEproc; // Number of elements on the current MPI rank
161  int ls_cut_marker; // Flag used for the cut-cell corresponding to the
162  // level set.
163 
164  // these are not thread-safe!
167 
168 public:
171  const double alpha_,
172  VectorCoefficient &vD_,
173  Array<int> &elem_marker_,
174  bool include_cut_cell_ = false,
175  int nterms_ = 0,
176  int ls_cut_marker_ = ShiftedFaceMarker::SBElementType::CUT)
177  : uD(&u), alpha(alpha_), vD(&vD_),
178  elem_marker(&elem_marker_),
179  include_cut_cell(include_cut_cell_),
180  nterms(nterms_),
181  NEproc(pmesh->GetNE()),
183  ls_cut_marker(ls_cut_marker_) { }
184 
185  virtual void AssembleRHSElementVect(const FiniteElement &el,
187  Vector &elvect);
188  virtual void AssembleRHSElementVect(const FiniteElement &el,
190  Vector &elvect);
191  virtual void AssembleRHSElementVect(const FiniteElement &el1,
192  const FiniteElement &el2,
194  Vector &elvect);
195 };
196 
197 
198 /// BilinearFormIntegrator for Neumann boundaries using the shifted boundary
199 /// method.
200 /// A(u,w) = <[nabla u + nabla(nabla u).d + h.o.t.].nhat(n.nhat),w>-<grad u.n,w>
201 /// where h.o.t are the high-order terms due to Taylor expansion for nabla u,
202 /// nhat is the normal vector at the true boundary, n is the normal vector at
203 /// the surrogate boundary. Since this interior face integrator is applied to
204 /// the surrogate boundary (see marking.hpp for notes on how the surrogate faces
205 /// are determined and elements are marked), this integrator adds contribution
206 /// to only the element that is adjacent to that face (Trans.Elem1 or
207 /// Trans.Elem2) and is part of the surrogate domain.
209 {
210 protected:
211  ShiftedVectorFunctionCoefficient *vN; // Normal function coefficient
212  VectorCoefficient *vD; // Distance function coefficient
213  Array<int> *elem_marker; // Marker indicating whether element is inside,
214  //cut, or outside the domain.
216  int nterms; // Number of terms in addition to the gradient
217  // term from Taylor expansion that should be included. (0 by default).
218  int NEproc; // Number of elements on the current MPI rank
221 
222 
223  // these are not thread-safe!
226 
227 
228 public:
230  VectorCoefficient &vD_,
232  Array<int> &elem_marker_,
233  Array<int> &cut_marker_,
234  bool include_cut_cell_ = false,
235  int nterms_ = 1)
236  : vN(&vN_), vD(&vD_),
237  elem_marker(&elem_marker_),
238  include_cut_cell(include_cut_cell_),
239  nterms(nterms_),
240  NEproc(pmesh->GetNE()),
242  cut_marker(cut_marker_) { }
243 
245  virtual void AssembleFaceMatrix(const FiniteElement &el1,
246  const FiniteElement &el2,
248  DenseMatrix &elmat);
249 
250  bool GetTrimFlag() const { return include_cut_cell; }
251 
253 };
254 
255 /// LinearFormIntegrator for Neumann boundaries using the shifted boundary
256 /// method.
257 /// (u, w) = <nhat.n t_n, w>
258 /// where nhat is the normal vector at the true boundary, n is the normal vector
259 /// at the surrogate boundary, and t_n is the traction boundary condition.
260 /// Since this interior face integrator is applied to the surrogate boundary
261 /// (see marking.hpp for notes on how the surrogate faces are determined and
262 /// elements are marked), this integrator adds contribution to only the element
263 /// that is adjacent to that face (Trans.Elem1 or Trans.Elem2) and is part of
264 /// the surrogate domain.
265 /// Note that t_N is evaluated at the true boundary using the distance function
266 /// and ShiftedFunctionCoefficient, i.e. t_N(x_true) = t_N(x_surrogate + D),
267 /// where x_surrogate is the location of the integration point on the surrogate
268 /// boundary and D is the distance vector from the surrogate boundary to the
269 /// true boundary.
271 {
272 protected:
273  ShiftedVectorFunctionCoefficient *vN; // Normal function coefficient
274  ShiftedFunctionCoefficient *uN; // Neumann condition on true boundary
275  VectorCoefficient *vD; // Distance function coefficient
276  Array<int> *elem_marker; // Marker indicating whether element is inside,
277  int nterms; // Number of terms in addition to the gradient
278  // term from Taylor expansion that should be included. (0 by default).
280  int NEproc; // Number of elements on the current MPI rank
283 
284  // these are not thread-safe!
286 
287 public:
290  VectorCoefficient &vD_,
292  Array<int> &elem_marker_,
293  int nterms_ = 0,
294  bool include_cut_cell_ = false,
295  int ls_cut_marker_ = ShiftedFaceMarker::SBElementType::CUT)
296  : vN(&vN_), uN(&u), vD(&vD_),
297  elem_marker(&elem_marker_),
298  nterms(nterms_),
299  include_cut_cell(include_cut_cell_),
300  NEproc(pmesh->GetNE()),
302  ls_cut_marker(ls_cut_marker_) { }
303 
304  virtual void AssembleRHSElementVect(const FiniteElement &el,
306  Vector &elvect);
307  virtual void AssembleRHSElementVect(const FiniteElement &el,
309  Vector &elvect);
310  virtual void AssembleRHSElementVect(const FiniteElement &el1,
311  const FiniteElement &el2,
313  Vector &elvect);
314  bool GetTrimFlag() const { return include_cut_cell; }
315 };
316 
317 } // namespace mfem
318 
319 #endif
Abstract class for all finite elements.
Definition: fe_base.hpp:235
Base class for vector Coefficients that optionally depend on time and space.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
SBM2NeumannIntegrator(const ParMesh *pmesh, VectorCoefficient &vD_, ShiftedVectorFunctionCoefficient &vN_, Array< int > &elem_marker_, Array< int > &cut_marker_, bool include_cut_cell_=false, int nterms_=1)
Definition: sbm_solver.hpp:229
SBM2NeumannLFIntegrator(const ParMesh *pmesh, ShiftedFunctionCoefficient &u, VectorCoefficient &vD_, ShiftedVectorFunctionCoefficient &vN_, Array< int > &elem_marker_, int nterms_=0, bool include_cut_cell_=false, int ls_cut_marker_=ShiftedFaceMarker::SBElementType::CUT)
Definition: sbm_solver.hpp:288
ShiftedVectorFunctionCoefficient * vN
Definition: sbm_solver.hpp:273
VectorCoefficient * vD
Definition: sbm_solver.hpp:92
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
SBM2DirichletLFIntegrator(const ParMesh *pmesh, ShiftedFunctionCoefficient &u, const double alpha_, VectorCoefficient &vD_, Array< int > &elem_marker_, bool include_cut_cell_=false, int nterms_=0, int ls_cut_marker_=ShiftedFaceMarker::SBElementType::CUT)
Definition: sbm_solver.hpp:169
VectorCoefficient * vD
Definition: sbm_solver.hpp:275
Abstract base class LinearFormIntegrator.
Definition: lininteg.hpp:22
SBM2DirichletIntegrator(const ParMesh *pmesh, const double a, VectorCoefficient &vD_, Array< int > &elem_marker_, Array< int > &cut_marker_, bool include_cut_cell_=false, int nterms_=0)
Definition: sbm_solver.hpp:109
std::function< void(const Vector &, Vector &)> Function
Definition: sbm_solver.hpp:54
ShiftedFunctionCoefficient(double constant_)
Definition: sbm_solver.hpp:33
ShiftedFunctionCoefficient * uN
Definition: sbm_solver.hpp:274
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: sbm_solver.cpp:338
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: sbm_solver.cpp:646
ShiftedVectorFunctionCoefficient * vN
Definition: sbm_solver.hpp:211
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
Definition: sbm_solver.hpp:62
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
ShiftedVectorFunctionCoefficient(int dim, std::function< void(const Vector &, Vector &)> F)
Definition: sbm_solver.hpp:57
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
ShiftedFunctionCoefficient(std::function< double(const Vector &v)> F)
Definition: sbm_solver.hpp:31
VectorCoefficient * vD
Definition: sbm_solver.hpp:212
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: sbm_solver.hpp:36
std::function< double(const Vector &)> Function
Definition: sbm_solver.hpp:26
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
double a
Definition: lissajous.cpp:41
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:139
Class for integration point with weight.
Definition: intrules.hpp:25
int dim
Definition: ex24.cpp:53
virtual void AssembleFaceMatrix(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Trans, DenseMatrix &elmat)
Definition: sbm_solver.cpp:48
Vector data type.
Definition: vector.hpp:60
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
virtual void AssembleRHSElementVect(const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
Definition: sbm_solver.cpp:932
Class for parallel meshes.
Definition: pmesh.hpp:32
ShiftedFunctionCoefficient * uD
Definition: sbm_solver.hpp:151