MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
eltrans.cpp
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 #include "../mesh/mesh_headers.hpp"
13 #include "fem.hpp"
14 #include <cmath>
15 
16 namespace mfem
17 {
18 
20  : IntPoint(static_cast<IntegrationPoint *>(NULL)),
21  EvalState(0),
22  Attribute(-1),
23  ElementNo(-1)
24 { }
25 
27 {
28  MFEM_ASSERT((EvalState & WEIGHT_MASK) == 0, "");
29  Jacobian();
31  return (Wght = (dFdx.Width() == 0) ? 1.0 : dFdx.Weight());
32 }
33 
35 {
36  MFEM_ASSERT((EvalState & ADJUGATE_MASK) == 0, "");
37  Jacobian();
39  if (dFdx.Width() > 0) { CalcAdjugate(dFdx, adjJ); }
41  return adjJ;
42 }
43 
45 {
46  // TODO: compute as invJ = / adjJ/Weight, if J is square,
47  // \ adjJ/Weight^2, otherwise.
48  MFEM_ASSERT((EvalState & INVERSE_MASK) == 0, "");
49  Jacobian();
51  if (dFdx.Width() > 0) { CalcInverse(dFdx, invJ); }
53  return invJ;
54 }
55 
56 
58  const Vector& pt, const IntegrationRule &ir)
59 {
60  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
61  MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point");
62 
63  DenseMatrix physPts;
64  T->Transform(ir, physPts);
65 
66  // Initialize distance and index of closest point
67  int minIndex = -1;
68  double minDist = std::numeric_limits<double>::max();
69 
70  // Check all integration points in ir
71  const int npts = ir.GetNPoints();
72  for (int i = 0; i < npts; ++i)
73  {
74  double dist = pt.DistanceTo(physPts.GetColumn(i));
75  if (dist < minDist)
76  {
77  minDist = dist;
78  minIndex = i;
79  }
80  }
81  return minIndex;
82 }
83 
85  const Vector& pt, const IntegrationRule &ir)
86 {
87  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
88  MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point");
89 
90  // Initialize distance and index of closest point
91  int minIndex = -1;
92  double minDist = std::numeric_limits<double>::max();
93 
94  // Check all integration points in ir using the local metric at each point
95  // induced by the transformation.
96  Vector dp(T->GetSpaceDim()), dr(T->GetDimension());
97  const int npts = ir.GetNPoints();
98  for (int i = 0; i < npts; ++i)
99  {
100  const IntegrationPoint &ip = ir.IntPoint(i);
101  T->Transform(ip, dp);
102  dp -= pt;
103  T->SetIntPoint(&ip);
104  T->InverseJacobian().Mult(dp, dr);
105  double dist = dr.Norml2();
106  // double dist = dr.Normlinf();
107  if (dist < minDist)
108  {
109  minDist = dist;
110  minIndex = i;
111  }
112  }
113  return minIndex;
114 }
115 
117 {
118  std::ostream &out = mfem::out;
119 
120  // separator:
121  switch (mode%3)
122  {
123  case 0: out << ", "; break;
124  case 1: out << "Newton: "; break;
125  case 2: out << " "; break;
126  // "Newton: iter = xx, "
127  }
128  switch ((mode/3)%4)
129  {
130  case 0: out << "iter = " << std::setw(2) << int(val); break;
131  case 1: out << "delta_ref = " << std::setw(11) << val; break;
132  case 2: out << " err_phys = " << std::setw(11) << val; break;
133  case 3: break;
134  }
135  // ending:
136  switch ((mode/12)%4)
137  {
138  case 0: break;
139  case 1: out << '\n'; break;
140  case 2: out << " (converged)\n"; break;
141  case 3: out << " (actual)\n"; break;
142  }
143 }
144 
146  const Vector &pt,
147  const char *suffix)
148 {
149  std::ostream &out = mfem::out;
150 
151  out << prefix << " = (";
152  for (int j = 0; j < pt.Size(); j++)
153  {
154  out << (j > 0 ? ", " : "") << pt(j);
155  }
156  out << ')' << suffix;
157 }
158 
160  IntegrationPoint &ip)
161 {
162  MFEM_ASSERT(pt.Size() == T->GetSpaceDim(), "invalid point");
163 
164  const double phys_tol = phys_rtol*pt.Normlinf();
165 
166  const int geom = T->GetGeometryType();
167  const int dim = T->GetDimension();
168  const int sdim = T->GetSpaceDim();
169  IntegrationPoint xip, prev_xip;
170  double xd[3], yd[3], dxd[3], dx_norm = -1.0, err_phys, real_dx_norm;
171  Vector x(xd, dim), y(yd, sdim), dx(dxd, dim);
172  bool hit_bdr = false, prev_hit_bdr;
173 
174  // Use ip0 as initial guess:
175  xip = *ip0;
176  xip.Get(xd, dim); // xip -> x
177  if (print_level >= 3)
178  {
179  NewtonPrint(1, 0.); // iter 0
180  NewtonPrintPoint(", ref_pt", x, "\n");
181  }
182 
183  for (int it = 0; true; )
184  {
185  // Remarks:
186  // If f(x) := 1/2 |pt-F(x)|^2, then grad(f)(x) = -J^t(x) [pt-F(x)].
187  // Linearize F(y) at y=x: F(y) ~ L[x](y) := F(x) + J(x) [y-x].
188  // Newton iteration for F(y)=b is given by L[x_old](x_new) = b, i.e.
189  // F(x_old) + J(x_old) [x_new-x_old] = b.
190  //
191  // To minimize: 1/2 |F(y)-b|^2, subject to: l(y) >= 0, we may consider the
192  // iteration: minimize: |L[x_old](x_new)-b|^2, subject to l(x_new) >= 0,
193  // i.e. minimize: |F(x_old) + J(x_old) [x_new-x_old] - b|^2.
194 
195  // This method uses:
196  // Newton iteration: x := x + J(x)^{-1} [pt-F(x)]
197  // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)]
198 
199  // Compute the physical coordinates of the current point:
200  T->Transform(xip, y);
201  if (print_level >= 3)
202  {
203  NewtonPrint(11, 0.); // continuation line
204  NewtonPrintPoint("approx_pt", y, ", ");
205  NewtonPrintPoint("exact_pt", pt, "\n");
206  }
207  subtract(pt, y, y); // y = pt-y
208 
209  // Check for convergence in physical coordinates:
210  err_phys = y.Normlinf();
211  if (err_phys < phys_tol)
212  {
213  if (print_level >= 1)
214  {
215  NewtonPrint(1, (double)it);
216  NewtonPrint(3, dx_norm);
217  NewtonPrint(30, err_phys);
218  }
219  ip = xip;
220  if (solver_type != Newton) { return Inside; }
221  return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside;
222  }
223  if (print_level >= 1)
224  {
225  if (it == 0 || print_level >= 2)
226  {
227  NewtonPrint(1, (double)it);
228  NewtonPrint(3, dx_norm);
229  NewtonPrint(18, err_phys);
230  }
231  }
232 
233  if (hit_bdr)
234  {
235  xip.Get(xd, dim); // xip -> x
236  if (prev_hit_bdr || it == max_iter || print_level >= 2)
237  {
238  prev_xip.Get(dxd, dim); // prev_xip -> dx
239  subtract(x, dx, dx); // dx = xip - prev_xip
240  real_dx_norm = dx.Normlinf();
241  if (print_level >= 2)
242  {
243  NewtonPrint(41, real_dx_norm);
244  }
245  if (prev_hit_bdr && real_dx_norm < ref_tol)
246  {
247  if (print_level >= 0)
248  {
249  if (print_level <= 1)
250  {
251  NewtonPrint(1, (double)it);
252  NewtonPrint(3, dx_norm);
253  NewtonPrint(18, err_phys);
254  NewtonPrint(41, real_dx_norm);
255  }
256  mfem::out << "Newton: *** stuck on boundary!\n";
257  }
258  return Outside;
259  }
260  }
261  }
262 
263  if (it == max_iter) { break; }
264 
265  // Perform a Newton step:
266  T->SetIntPoint(&xip);
267  T->InverseJacobian().Mult(y, dx);
268  x += dx;
269  it++;
270  if (solver_type != Newton)
271  {
272  prev_xip = xip;
273  prev_hit_bdr = hit_bdr;
274  }
275  xip.Set(xd, dim); // x -> xip
276 
277  // Perform projection based on solver_type:
278  switch (solver_type)
279  {
280  case Newton: break;
282  hit_bdr = !Geometry::ProjectPoint(geom, prev_xip, xip); break;
284  hit_bdr = !Geometry::ProjectPoint(geom, xip); break;
285  default: MFEM_ABORT("invalid solver type");
286  }
287  if (print_level >= 3)
288  {
289  NewtonPrint(1, double(it));
290  xip.Get(xd, dim); // xip -> x
291  NewtonPrintPoint(", ref_pt", x, "\n");
292  }
293 
294  // Check for convergence in reference coordinates:
295  dx_norm = dx.Normlinf();
296  if (dx_norm < ref_tol)
297  {
298  if (print_level >= 1)
299  {
300  NewtonPrint(1, (double)it);
301  NewtonPrint(27, dx_norm);
302  }
303  ip = xip;
304  if (solver_type != Newton) { return Inside; }
305  return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside;
306  }
307  }
308  if (print_level >= 0)
309  {
310  if (print_level <= 1)
311  {
312  NewtonPrint(1, (double)max_iter);
313  NewtonPrint(3, dx_norm);
314  NewtonPrint(18, err_phys);
315  if (hit_bdr) { NewtonPrint(41, real_dx_norm); }
316  }
317  mfem::out << "Newton: *** iteration did not converge!\n";
318  }
319  ip = xip;
320  return Unknown;
321 }
322 
324  IntegrationPoint &ip)
325 {
326  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
327 
328  // Select initial guess ...
329  switch (init_guess_type)
330  {
331  case Center:
333  break;
334 
335  case ClosestPhysNode:
336  case ClosestRefNode:
337  {
338  const int order = std::max(T->Order()+rel_qpts_order, 0);
339  if (order == 0)
340  {
342  }
343  else
344  {
345  const int old_type = GlobGeometryRefiner.GetType();
347  RefinedGeometry &RefG =
349  int closest_idx = (init_guess_type == ClosestPhysNode) ?
350  FindClosestPhysPoint(pt, RefG.RefPts) :
351  FindClosestRefPoint(pt, RefG.RefPts);
352  ip0 = &RefG.RefPts.IntPoint(closest_idx);
353  GlobGeometryRefiner.SetType(old_type);
354  }
355  break;
356  }
357 
358  case GivenPoint:
359  break;
360 
361  default:
362  MFEM_ABORT("invalid initial guess type");
363  }
364 
365  // Call the solver ...
366  return NewtonSolve(pt, ip);
367 }
368 
369 
371 {
372  switch (GeomType)
373  {
374  case Geometry::POINT : FElem = &PointFE; break;
375  case Geometry::SEGMENT : FElem = &SegmentFE; break;
376  case Geometry::TRIANGLE : FElem = &TriangleFE; break;
377  case Geometry::SQUARE : FElem = &QuadrilateralFE; break;
378  case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
379  case Geometry::CUBE : FElem = &HexahedronFE; break;
380  default:
381  MFEM_ABORT("unknown Geometry::Type!");
382  }
383  int dim = FElem->GetDim();
384  int dof = FElem->GetDof();
385  const IntegrationRule &nodes = FElem->GetNodes();
386  PointMat.SetSize(dim, dof);
387  for (int j = 0; j < dof; j++)
388  {
389  nodes.IntPoint(j).Get(&PointMat(0,j), dim);
390  }
391  geom = GeomType;
392  space_dim = dim;
393 }
394 
395 const DenseMatrix &IsoparametricTransformation::EvalJacobian()
396 {
397  MFEM_ASSERT(space_dim == PointMat.Height(),
398  "the IsoparametricTransformation has not been finalized;"
399  " call FinilizeTransformation() after setup");
400  MFEM_ASSERT((EvalState & JACOBIAN_MASK) == 0, "");
401 
402  dshape.SetSize(FElem->GetDof(), FElem->GetDim());
403  dFdx.SetSize(PointMat.Height(), dshape.Width());
404  if (dshape.Width() > 0)
405  {
406  FElem->CalcDShape(*IntPoint, dshape);
407  Mult(PointMat, dshape, dFdx);
408  }
410 
411  return dFdx;
412 }
413 
415 {
416  switch (FElem->Space())
417  {
418  case FunctionSpace::Pk:
419  return (FElem->GetOrder()-1);
420  case FunctionSpace::Qk:
421  return (FElem->GetOrder());
422  default:
423  mfem_error("IsoparametricTransformation::OrderJ()");
424  }
425  return 0;
426 }
427 
429 {
430  switch (FElem->Space())
431  {
432  case FunctionSpace::Pk:
433  return (FElem->GetOrder() - 1) * FElem->GetDim();
434  case FunctionSpace::Qk:
435  return (FElem->GetOrder() * FElem->GetDim() - 1);
436  default:
437  mfem_error("IsoparametricTransformation::OrderW()");
438  }
439  return 0;
440 }
441 
443 {
444  if (FElem->Space() == fe->Space())
445  {
446  int k = FElem->GetOrder();
447  int d = FElem->GetDim();
448  int l = fe->GetOrder();
449  switch (fe->Space())
450  {
451  case FunctionSpace::Pk:
452  return ((k-1)*(d-1)+(l-1));
453  case FunctionSpace::Qk:
454  return (k*(d-1)+(l-1));
455  }
456  }
457  mfem_error("IsoparametricTransformation::OrderGrad(...)");
458  return 0;
459 }
460 
462  Vector &trans)
463 {
464  shape.SetSize(FElem->GetDof());
465  trans.SetSize(PointMat.Height());
466 
467  FElem -> CalcShape(ip, shape);
468  PointMat.Mult(shape, trans);
469 }
470 
472  DenseMatrix &tr)
473 {
474  int dof, n, dim, i, j, k;
475 
476  dim = PointMat.Height();
477  dof = FElem->GetDof();
478  n = ir.GetNPoints();
479 
480  shape.SetSize(dof);
481  tr.SetSize(dim, n);
482 
483  for (j = 0; j < n; j++)
484  {
485  FElem -> CalcShape (ir.IntPoint(j), shape);
486  for (i = 0; i < dim; i++)
487  {
488  tr(i, j) = 0.0;
489  for (k = 0; k < dof; k++)
490  {
491  tr(i, j) += PointMat(i, k) * shape(k);
492  }
493  }
494  }
495 }
496 
498  DenseMatrix &result)
499 {
500  MFEM_ASSERT(matrix.Height() == GetDimension(), "invalid input");
501  result.SetSize(PointMat.Height(), matrix.Width());
502 
503  IntegrationPoint ip;
504  Vector col;
505 
506  for (int j = 0; j < matrix.Width(); j++)
507  {
508  ip.Set(matrix.GetColumn(j), matrix.Height());
509 
510  result.GetColumnReference(j, col);
511  Transform(ip, col);
512  }
513 }
514 
516  IntegrationPoint &ip2)
517 {
518  double vec[3];
519  Vector v (vec, Transf.GetPointMat().Height());
520 
521  Transf.Transform (ip1, v);
522  ip2.Set(vec, v.Size());
523 }
524 
526  IntegrationRule &ir2)
527 {
528  int i, n;
529 
530  n = ir1.GetNPoints();
531  for (i = 0; i < n; i++)
532  {
533  Transform (ir1.IntPoint(i), ir2.IntPoint(i));
534  }
535 }
536 
537 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
Abstract class for Finite Elements.
Definition: fe.hpp:47
void Set(const double *p, const int dim)
Definition: intrules.hpp:32
int NewtonSolve(const Vector &pt, IntegrationPoint &ip)
Definition: eltrans.cpp:159
void Get(double *p, const int dim) const
Definition: intrules.hpp:45
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:116
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
const IntegrationPoint * IntPoint
Definition: eltrans.hpp:26
ElementTransformation * T
Definition: eltrans.hpp:154
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:320
const Geometry::Type geom
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
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:468
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:42
int GetDimension() const
Return the dimension of the reference element.
Definition: eltrans.hpp:88
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:52
RefinedGeometry * Refine(int Geom, int Times, int ETimes=1)
Definition: geom.cpp:802
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:3048
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
int Size() const
Returns the size of the vector.
Definition: vector.hpp:113
Polynomials of order k.
Definition: fe.hpp:34
const DenseMatrix & EvalAdjugateJ()
Definition: eltrans.cpp:34
int Space() const
Returns the type of space on each element.
Definition: fe.hpp:128
const DenseMatrix & InverseJacobian()
Definition: eltrans.hpp:75
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:707
virtual void Transform(const IntegrationPoint &, Vector &)
Definition: eltrans.cpp:461
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:61
void NewtonPrintPoint(const char *prefix, const Vector &pt, const char *suffix)
Definition: eltrans.cpp:145
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
double Weight() const
Definition: densemat.cpp:492
PointFiniteElement PointFE
Definition: point.cpp:30
Geometry Geometries
Definition: geom.cpp:759
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
int dim
Definition: ex3.cpp:47
The point is inside the element.
Definition: eltrans.hpp:147
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:36
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
Definition: geom.hpp:233
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
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:52
int GetSpaceDim() const
Get the dimension of the target (physical) space.
Definition: eltrans.hpp:93
const IntegrationRule & GetNodes() const
Definition: fe.hpp:167
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1188
IntegrationRule RefPts
Definition: geom.hpp:210
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:2312
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:408
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:3156
const DenseMatrix & EvalInverseJ()
Definition: eltrans.cpp:44
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:445
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:122
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:390
void mfem_error(const char *msg)
Definition: error.cpp:107
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:515
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:243
Linear3DFiniteElement TetrahedronFE
Linear2DFiniteElement TriangleFE
Definition: triangle.cpp:198
Class for integration point with weight.
Definition: intrules.hpp:25
The algorithm failed to determine where the point is.
Definition: eltrans.hpp:149
int GetType() const
Get the Quadrature1D type of points used for subdivision.
Definition: geom.hpp:235
int FindClosestPhysPoint(const Vector &pt, const IntegrationRule &ir)
Find the IntegrationPoint mapped closest to pt.
Definition: eltrans.cpp:57
IsoparametricTransformation Transf
Definition: eltrans.hpp:338
Vector data type.
Definition: vector.hpp:41
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:143
virtual void Transform(const IntegrationPoint &, Vector &)=0
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:296
const IntegrationPoint * ip0
Definition: eltrans.hpp:157
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:86
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:64
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
BiLinear2DFiniteElement QuadrilateralFE
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
Tensor products of polynomials of order k.
Definition: fe.hpp:35