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