MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
eltrans.hpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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_ELEMENTTRANSFORM
13 #define MFEM_ELEMENTTRANSFORM
14 
15 #include "../config/config.hpp"
16 #include "../linalg/linalg.hpp"
17 #include "intrules.hpp"
18 #include "fe.hpp"
19 
20 namespace mfem
21 {
22 
24 {
25 protected:
29  double Wght;
30  int EvalState;
32  {
38  };
40  int space_dim;
41 
42  // Evaluate the Jacobian of the transformation at the IntPoint and store it
43  // in dFdx.
44  virtual const DenseMatrix &EvalJacobian() = 0;
45  virtual const DenseMatrix &EvalHessian() = 0;
46 
47  double EvalWeight();
48  const DenseMatrix &EvalAdjugateJ();
49  const DenseMatrix &EvalInverseJ();
50 
51 public:
53 
55 
56  void SetIntPoint(const IntegrationPoint *ip)
57  { IntPoint = ip; EvalState = 0; }
58  const IntegrationPoint &GetIntPoint() { return *IntPoint; }
59 
60  virtual void Transform(const IntegrationPoint &, Vector &) = 0;
61  virtual void Transform(const IntegrationRule &, DenseMatrix &) = 0;
62 
63  /// Transform columns of 'matrix', store result in 'result'.
64  virtual void Transform(const DenseMatrix &matrix, DenseMatrix &result) = 0;
65 
66  /** @brief Return the Jacobian matrix of the transformation at the currently
67  set IntegrationPoint, using the method SetIntPoint(). */
68  /** The dimensions of the Jacobian matrix are physical-space-dim by
69  reference-space-dim. The first column contains the x derivatives of the
70  transformation, the second -- the y derivatives, etc. */
72  { return (EvalState & JACOBIAN_MASK) ? dFdx : EvalJacobian(); }
73 
75  { return (EvalState & HESSIAN_MASK) ? d2Fdx2 : EvalHessian(); }
76 
77  double Weight() { return (EvalState & WEIGHT_MASK) ? Wght : EvalWeight(); }
78 
80  { return (EvalState & ADJUGATE_MASK) ? adjJ : EvalAdjugateJ(); }
81 
83  { return (EvalState & INVERSE_MASK) ? invJ : EvalInverseJ(); }
84 
85  virtual int Order() = 0;
86  virtual int OrderJ() = 0;
87  virtual int OrderW() = 0;
88  /// Order of adj(J)^t.grad(fi)
89  virtual int OrderGrad(const FiniteElement *fe) = 0;
90 
91  /// Return the Geometry::Type of the reference element.
92  Geometry::Type GetGeometryType() const { return geom; }
93 
94  /// Return the dimension of the reference element.
95  int GetDimension() const { return Geometry::Dimension[geom]; }
96 
97  /// Get the dimension of the target (physical) space.
98  /** We support 2D meshes embedded in 3D; in this case the function will
99  return "3". */
100  int GetSpaceDim() const { return space_dim; }
101 
102  /** @brief Transform a point @a pt from physical space to a point @a ip in
103  reference space. */
104  /** Attempt to find the IntegrationPoint that is transformed into the given
105  point in physical space. If the inversion fails a non-zero value is
106  returned. This method is not 100 percent reliable for non-linear
107  transformations. */
108  virtual int TransformBack(const Vector &pt, IntegrationPoint &ip) = 0;
109 
111 };
112 
113 
114 /// The inverse transformation of a given ElementTransformation.
116 {
117 public:
118  /// Algorithms for selecting an initial guess.
120  {
121  Center = 0, ///< Use the center of the reference element.
122  ClosestPhysNode = 1, /**<
123  Use the point returned by FindClosestPhysPoint() from a reference-space
124  grid of type and size controlled by SetInitGuessPointsType() and
125  SetInitGuessRelOrder(), respectively. */
126  ClosestRefNode = 2, /**<
127  Use the point returned by FindClosestRefPoint() from a reference-space
128  grid of type and size controlled by SetInitGuessPointsType() and
129  SetInitGuessRelOrder(), respectively. */
130  GivenPoint = 3 ///< Use a specific point, set with SetInitialGuess().
131  };
132 
133  /// Solution strategy.
135  {
136  Newton = 0, /**<
137  Use Newton's algorithm, without restricting the reference-space points
138  (iterates) to the reference element. */
140  Use Newton's algorithm, restricting the reference-space points to the
141  reference element by scaling back the Newton increments, i.e.
142  projecting new iterates, x_new, lying outside the element, to the
143  intersection of the line segment [x_old, x_new] with the boundary. */
145  Use Newton's algorithm, restricting the reference-space points to the
146  reference element by projecting new iterates, x_new, lying outside the
147  element, to the point on the boundary closest (in reference-space) to
148  x_new. */
149  };
150 
151  /// Values returned by Transform().
153  {
154  Inside = 0, ///< The point is inside the element
155  Outside = 1, ///< The point is _probably_ outside the element
156  Unknown = 2 ///< The algorithm failed to determine where the point is
157  };
158 
159 protected:
160  // Pointer to the forward transformation. Not owned.
162 
163  // Parameters of the inversion algorithms:
165  int init_guess_type; // algorithm to use
166  int qpts_type; // Quadrature1D type for the initial guess type
167  int rel_qpts_order; // num_1D_qpts = max(trans_order+rel_qpts_order,0)+1
168  int solver_type; // solution strategy to use
169  int max_iter; // max. number of Newton iterations
170  double ref_tol; // reference space tolerance
171  double phys_rtol; // physical space tolerance (relative)
172  double ip_tol; // tolerance for checking if a point is inside the ref. elem.
174 
175  void NewtonPrint(int mode, double val);
176  void NewtonPrintPoint(const char *prefix, const Vector &pt,
177  const char *suffix);
178  int NewtonSolve(const Vector &pt, IntegrationPoint &ip);
179 
180 public:
181  /// Construct the InverseElementTransformation with default parameters.
182  /** Some practical considerations regarding the choice of initial guess type
183  and solver type:
184  1. The combination of #Center and #NewtonSegmentProject provides the
185  fastest way to estimate if a point lies inside an element, assuming
186  that most queried elements are not very deformed, e.g. if most
187  elements are convex.
188  2. [Default] The combination of #Center and #NewtonElementProject
189  provides a somewhat slower alternative to 1 with the benefit of being
190  more reliable in the case when the query point is inside the element
191  but potentially slower in the case when the query point is outside the
192  element.
193  3. The combination of #ClosestPhysNode and #NewtonElementProject is
194  slower than 1 and 2 but much more reliable, especially in the case of
195  highly distorted elements which do not have very high aspect ratios.
196  4. The combination of #ClosestRefNode and #NewtonElementProject should
197  generally be the most reliable, coming at a bit higher computational
198  cost than 3 while performing better on distorted meshes with elements
199  having high aspect ratios.
200  @note None of these choices provide a guarantee that if a query point is
201  inside the element then it will be found. The only guarantee is that if
202  the Transform() method returns #Inside then the point lies inside the
203  element up to one of the specified physical- or reference-space
204  tolerances. */
206  : T(Trans),
207  ip0(NULL),
209  qpts_type(Quadrature1D::OpenHalfUniform),
210  rel_qpts_order(-1),
212  max_iter(16),
213  ref_tol(1e-15),
214  phys_rtol(1e-15),
215  ip_tol(1e-8),
216  print_level(-1)
217  { }
218 
220 
221  /// Set a new forward ElementTransformation, @a Trans.
223 
224  /** @brief Choose how the initial guesses for subsequent calls to Transform()
225  will be selected. */
227 
228  /** @brief Set the initial guess for subsequent calls to Transform(),
229  switching to the #GivenPoint #InitGuessType at the same time. */
230  void SetInitialGuess(const IntegrationPoint &init_ip)
231  { ip0 = &init_ip; SetInitialGuessType(GivenPoint); }
232 
233  /// Set the Quadrature1D type used for the `Closest*` initial guess types.
234  void SetInitGuessPointsType(int q_type) { qpts_type = q_type; }
235 
236  /// Set the relative order used for the `Closest*` initial guess types.
237  /** The number of points in each spatial direction is given by the formula
238  max(trans_order+order,0)+1, where trans_order is the order of the current
239  ElementTransformation. */
240  void SetInitGuessRelOrder(int order) { rel_qpts_order = order; }
241 
242  /** @brief Specify which algorithm to use for solving the transformation
243  equation, i.e. when calling the Transform() method. */
244  void SetSolverType(SolverType stype) { solver_type = stype; }
245 
246  /// Set the maximum number of iterations when solving for a reference point.
247  void SetMaxIter(int max_it) { max_iter = max_it; }
248 
249  /// Set the reference-space convergence tolerance.
250  void SetReferenceTol(double ref_sp_tol) { ref_tol = ref_sp_tol; }
251 
252  /// Set the relative physical-space convergence tolerance.
253  void SetPhysicalRelTol(double phys_rel_tol) { phys_rtol = phys_rel_tol; }
254 
255  /** @brief Set the tolerance used to determine if a point lies inside or
256  outside of the reference element. */
257  /** This tolerance is used only with the pure #Newton solver. */
258  void SetElementTol(double el_tol) { ip_tol = el_tol; }
259 
260  /// Set the desired print level, useful for debugging.
261  /** The valid options are: -1 - never print (default); 0 - print only errors;
262  1 - print the first and last last iterations; 2 - print every iteration;
263  and 3 - print every iteration including point coordinates. */
264  void SetPrintLevel(int pr_level) { print_level = pr_level; }
265 
266  /** @brief Find the IntegrationPoint mapped closest to @a pt. */
267  /** This function uses the given IntegrationRule, @a ir, maps its points to
268  physical space and finds the one that is closest to the point @a pt.
269 
270  @param pt The query point.
271  @param ir The IntegrationRule, i.e. the set of reference points to map
272  to physical space and check.
273  @return The index of the IntegrationPoint in @a ir whose mapped point is
274  closest to @a pt.
275  @see FindClosestRefPoint(). */
276  int FindClosestPhysPoint(const Vector& pt, const IntegrationRule &ir);
277 
278  /** @brief Find the IntegrationPoint mapped closest to @a pt, using a norm
279  that approximates the (unknown) distance in reference coordinates. */
280  /** @see FindClosestPhysPoint(). */
281  int FindClosestRefPoint(const Vector& pt, const IntegrationRule &ir);
282 
283  /** @brief Given a point, @a pt, in physical space, find its reference
284  coordinates, @a ip.
285 
286  @returns A value of type #TransformResult. */
287  virtual int Transform(const Vector &pt, IntegrationPoint &ip);
288 };
289 
290 
292 {
293 private:
294  DenseMatrix dshape,d2shape;
295  Vector shape;
296 
297  const FiniteElement *FElem;
298  DenseMatrix PointMat; // dim x dof
299 
300  // Evaluate the Jacobian of the transformation at the IntPoint and store it
301  // in dFdx.
302  virtual const DenseMatrix &EvalJacobian();
303  // Evaluate the Hessian of the transformation at the IntPoint and store it
304  // in d2Fdx2.
305  virtual const DenseMatrix &EvalHessian();
306 public:
307  void SetFE(const FiniteElement *FE) { FElem = FE; geom = FE->GetGeomType(); }
308  const FiniteElement* GetFE() const { return FElem; }
309 
310  /** @brief Read and write access to the underlying point matrix describing
311  the transformation. */
312  /** The dimensions of the matrix are space-dim x dof. The transformation is
313  defined as
314 
315  x=F(xh)=P.phi(xh),
316 
317  where xh (x hat) is the reference point, x is the corresponding physical
318  point, P is the point matrix, and phi(xh) is the column-vector of all
319  basis functions evaluated at xh. The columns of P represent the control
320  points in physical space defining the transformation. */
321  DenseMatrix &GetPointMat() { return PointMat; }
322  void FinalizeTransformation() { space_dim = PointMat.Height(); }
323 
325 
326  virtual void Transform(const IntegrationPoint &, Vector &);
327  virtual void Transform(const IntegrationRule &, DenseMatrix &);
328  virtual void Transform(const DenseMatrix &matrix, DenseMatrix &result);
329 
330  virtual int Order() { return FElem->GetOrder(); }
331  virtual int OrderJ();
332  virtual int OrderW();
333  virtual int OrderGrad(const FiniteElement *fe);
334 
335  virtual int TransformBack(const Vector & v, IntegrationPoint & ip)
336  {
337  InverseElementTransformation inv_tr(this);
338  return inv_tr.Transform(v, ip);
339  }
340 
342 };
343 
345 {
346 public:
348  void Transform (const IntegrationPoint &, IntegrationPoint &);
349  void Transform (const IntegrationRule &, IntegrationRule &);
350 };
351 
353 {
354 public:
358 };
359 
360 /* Elem1(Loc1(x)) = Face(x) = Elem2(Loc2(x))
361 
362 
363  Physical Space
364 
365  *--------* ^ *--------*
366  Elem1No / / \ / \ / \ \ Elem2No
367  / / \ / \ / \ \
368  / / n \ / \ / \ \
369  *--------* ==> * ( ) * *--------*
370  \ \ / \ / \ / /
371  \ \ / \ / \ / /
372  \ \ / \ / \ / /
373  *--------* v *--------*
374 
375  ^ ^
376  | ^ |
377  Elem1 | | | Elem2
378  | | Face |
379  |
380  *--------* *--------*
381  / /| / /|
382  1 *--------* | 1 *--------* 1 *--------* |
383  | | | Loc1 | | Loc2 | | |
384  | | * <----- | x | -----> | | *
385  | |/ | | | |/
386  *--------* *--------* *--------*
387  0 1 0 1 0 1
388 
389  Reference Space
390 */
391 
392 }
393 
394 #endif
Abstract class for Finite Elements.
Definition: fe.hpp:232
virtual ~ElementTransformation()
Definition: eltrans.hpp:110
void SetPrintLevel(int pr_level)
Set the desired print level, useful for debugging.
Definition: eltrans.hpp:264
int NewtonSolve(const Vector &pt, IntegrationPoint &ip)
Definition: eltrans.cpp:159
const DenseMatrix & AdjugateJacobian()
Definition: eltrans.hpp:79
ElementTransformation * Face
Definition: eltrans.hpp:356
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
DenseMatrix & GetPointMat()
Read and write access to the underlying point matrix describing the transformation.
Definition: eltrans.hpp:321
Use a specific point, set with SetInitialGuess().
Definition: eltrans.hpp:130
void SetPhysicalRelTol(double phys_rel_tol)
Set the relative physical-space convergence tolerance.
Definition: eltrans.hpp:253
void SetInitialGuess(const IntegrationPoint &init_ip)
Set the initial guess for subsequent calls to Transform(), switching to the GivenPoint InitGuessType ...
Definition: eltrans.hpp:230
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:92
const IntegrationPoint * IntPoint
Definition: eltrans.hpp:26
ElementTransformation * T
Definition: eltrans.hpp:161
Use the center of the reference element.
Definition: eltrans.hpp:121
virtual int Transform(const Vector &pt, IntegrationPoint &ip)
Given a point, pt, in physical space, find its reference coordinates, ip.
Definition: eltrans.cpp:323
virtual int TransformBack(const Vector &v, IntegrationPoint &ip)
Transform a point pt from physical space to a point ip in reference space.
Definition: eltrans.hpp:335
int GetDimension() const
Return the dimension of the reference element.
Definition: eltrans.hpp:95
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
void NewtonPrint(int mode, double val)
Definition: eltrans.cpp:116
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:318
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
void SetInitialGuessType(InitGuessType itype)
Choose how the initial guesses for subsequent calls to Transform() will be selected.
Definition: eltrans.hpp:226
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:357
const DenseMatrix & EvalAdjugateJ()
Definition: eltrans.cpp:34
const DenseMatrix & InverseJacobian()
Definition: eltrans.hpp:82
const FiniteElement * GetFE() const
Definition: eltrans.hpp:308
A class container for 1D quadrature type constants.
Definition: intrules.hpp:286
virtual void Transform(const IntegrationPoint &, Vector &)
Definition: eltrans.cpp:483
const IntegrationPoint & GetIntPoint()
Definition: eltrans.hpp:58
The inverse transformation of a given ElementTransformation.
Definition: eltrans.hpp:115
ElementTransformation * Elem2
Definition: eltrans.hpp:356
void SetSolverType(SolverType stype)
Specify which algorithm to use for solving the transformation equation, i.e. when calling the Transfo...
Definition: eltrans.hpp:244
void NewtonPrintPoint(const char *prefix, const Vector &pt, const char *suffix)
Definition: eltrans.cpp:145
void SetTransformation(ElementTransformation &Trans)
Set a new forward ElementTransformation, Trans.
Definition: eltrans.hpp:222
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:311
int FindClosestRefPoint(const Vector &pt, const IntegrationRule &ir)
Find the IntegrationPoint mapped closest to pt, using a norm that approximates the (unknown) distance...
Definition: eltrans.cpp:84
The point is inside the element.
Definition: eltrans.hpp:154
static const int Dimension[NumGeom]
Definition: geom.hpp:46
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:357
The point is probably outside the element.
Definition: eltrans.hpp:155
virtual int OrderGrad(const FiniteElement *fe)
Order of adj(J)^t.grad(fi)
Definition: eltrans.cpp:464
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:71
int GetSpaceDim() const
Get the dimension of the target (physical) space.
Definition: eltrans.hpp:100
void SetReferenceTol(double ref_sp_tol)
Set the reference-space convergence tolerance.
Definition: eltrans.hpp:250
void SetInitGuessPointsType(int q_type)
Set the Quadrature1D type used for the Closest* initial guess types.
Definition: eltrans.hpp:234
InverseElementTransformation(ElementTransformation *Trans=NULL)
Construct the InverseElementTransformation with default parameters.
Definition: eltrans.hpp:205
virtual const DenseMatrix & EvalHessian()=0
void SetElementTol(double el_tol)
Set the tolerance used to determine if a point lies inside or outside of the reference element...
Definition: eltrans.hpp:258
const DenseMatrix & EvalInverseJ()
Definition: eltrans.cpp:44
void SetIdentityTransformation(Geometry::Type GeomType)
Definition: eltrans.cpp:370
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:537
Class for integration point with weight.
Definition: intrules.hpp:25
The algorithm failed to determine where the point is.
Definition: eltrans.hpp:156
SolverType
Solution strategy.
Definition: eltrans.hpp:134
virtual int OrderGrad(const FiniteElement *fe)=0
Order of adj(J)^t.grad(fi)
TransformResult
Values returned by Transform().
Definition: eltrans.hpp:152
int FindClosestPhysPoint(const Vector &pt, const IntegrationRule &ir)
Find the IntegrationPoint mapped closest to pt.
Definition: eltrans.cpp:57
void SetFE(const FiniteElement *FE)
Definition: eltrans.hpp:307
IsoparametricTransformation Transf
Definition: eltrans.hpp:347
ElementTransformation * Elem1
Definition: eltrans.hpp:356
virtual const DenseMatrix & EvalJacobian()=0
InitGuessType
Algorithms for selecting an initial guess.
Definition: eltrans.hpp:119
Vector data type.
Definition: vector.hpp:48
virtual void Transform(const IntegrationPoint &, Vector &)=0
const IntegrationPoint * ip0
Definition: eltrans.hpp:164
const DenseMatrix & Hessian()
Definition: eltrans.hpp:74
void SetMaxIter(int max_it)
Set the maximum number of iterations when solving for a reference point.
Definition: eltrans.hpp:247
void SetInitGuessRelOrder(int order)
Set the relative order used for the Closest* initial guess types.
Definition: eltrans.hpp:240
virtual int TransformBack(const Vector &pt, IntegrationPoint &ip)=0
Transform a point pt from physical space to a point ip in reference space.