MFEM  v4.5.1
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-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 #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  geom(Geometry::INVALID),
23  Attribute(-1),
24  ElementNo(-1),
25  mesh(nullptr)
26 { }
27 
29 {
30  MFEM_ASSERT((EvalState & WEIGHT_MASK) == 0, "");
31  Jacobian();
33  return (Wght = (dFdx.Width() == 0) ? 1.0 : dFdx.Weight());
34 }
35 
37 {
38  MFEM_ASSERT((EvalState & ADJUGATE_MASK) == 0, "");
39  Jacobian();
41  if (dFdx.Width() > 0) { CalcAdjugate(dFdx, adjJ); }
43  return adjJ;
44 }
45 
47 {
48  MFEM_ASSERT((EvalState & TRANS_ADJUGATE_MASK) == 0, "");
49  Jacobian();
52  else { AdjugateJacobian(); adjJT.Transpose(adjJ); }
54  return adjJT;
55 }
56 
58 {
59  // TODO: compute as invJ = / adjJ/Weight, if J is square,
60  // \ adjJ/Weight^2, otherwise.
61  MFEM_ASSERT((EvalState & INVERSE_MASK) == 0, "");
62  Jacobian();
64  if (dFdx.Width() > 0) { CalcInverse(dFdx, invJ); }
66  return invJ;
67 }
68 
69 
71  const Vector& pt, const IntegrationRule &ir)
72 {
73  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
74  MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point");
75 
76  DenseMatrix physPts;
77  T->Transform(ir, physPts);
78 
79  // Initialize distance and index of closest point
80  int minIndex = -1;
81  double minDist = std::numeric_limits<double>::max();
82 
83  // Check all integration points in ir
84  const int npts = ir.GetNPoints();
85  for (int i = 0; i < npts; ++i)
86  {
87  double dist = pt.DistanceTo(physPts.GetColumn(i));
88  if (dist < minDist)
89  {
90  minDist = dist;
91  minIndex = i;
92  }
93  }
94  return minIndex;
95 }
96 
98  const Vector& pt, const IntegrationRule &ir)
99 {
100  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
101  MFEM_VERIFY(pt.Size() == T->GetSpaceDim(), "invalid point");
102 
103  // Initialize distance and index of closest point
104  int minIndex = -1;
105  double minDist = std::numeric_limits<double>::max();
106 
107  // Check all integration points in ir using the local metric at each point
108  // induced by the transformation.
109  Vector dp(T->GetSpaceDim()), dr(T->GetDimension());
110  const int npts = ir.GetNPoints();
111  for (int i = 0; i < npts; ++i)
112  {
113  const IntegrationPoint &ip = ir.IntPoint(i);
114  T->Transform(ip, dp);
115  dp -= pt;
116  T->SetIntPoint(&ip);
117  T->InverseJacobian().Mult(dp, dr);
118  double dist = dr.Norml2();
119  // double dist = dr.Normlinf();
120  if (dist < minDist)
121  {
122  minDist = dist;
123  minIndex = i;
124  }
125  }
126  return minIndex;
127 }
128 
130 {
131  std::ostream &os = mfem::out;
132 
133  // separator:
134  switch (mode%3)
135  {
136  case 0: os << ", "; break;
137  case 1: os << "Newton: "; break;
138  case 2: os << " "; break;
139  // "Newton: iter = xx, "
140  }
141  switch ((mode/3)%4)
142  {
143  case 0: os << "iter = " << std::setw(2) << int(val); break;
144  case 1: os << "delta_ref = " << std::setw(11) << val; break;
145  case 2: os << " err_phys = " << std::setw(11) << val; break;
146  case 3: break;
147  }
148  // ending:
149  switch ((mode/12)%4)
150  {
151  case 0: break;
152  case 1: os << '\n'; break;
153  case 2: os << " (converged)\n"; break;
154  case 3: os << " (actual)\n"; break;
155  }
156 }
157 
159  const Vector &pt,
160  const char *suffix)
161 {
162  std::ostream &os = mfem::out;
163 
164  os << prefix << " = (";
165  for (int j = 0; j < pt.Size(); j++)
166  {
167  os << (j > 0 ? ", " : "") << pt(j);
168  }
169  os << ')' << suffix;
170 }
171 
173  IntegrationPoint &ip)
174 {
175  MFEM_ASSERT(pt.Size() == T->GetSpaceDim(), "invalid point");
176 
177  const double phys_tol = phys_rtol*pt.Normlinf();
178 
179  const int geom = T->GetGeometryType();
180  const int dim = T->GetDimension();
181  const int sdim = T->GetSpaceDim();
182  IntegrationPoint xip, prev_xip;
183  double xd[3], yd[3], dxd[3], dx_norm = -1.0, err_phys, real_dx_norm = -1.0;
184  Vector x(xd, dim), y(yd, sdim), dx(dxd, dim);
185  bool hit_bdr = false, prev_hit_bdr = false;
186 
187  // Use ip0 as initial guess:
188  xip = *ip0;
189  xip.Get(xd, dim); // xip -> x
190  if (print_level >= 3)
191  {
192  NewtonPrint(1, 0.); // iter 0
193  NewtonPrintPoint(", ref_pt", x, "\n");
194  }
195 
196  for (int it = 0; true; )
197  {
198  // Remarks:
199  // If f(x) := 1/2 |pt-F(x)|^2, then grad(f)(x) = -J^t(x) [pt-F(x)].
200  // Linearize F(y) at y=x: F(y) ~ L[x](y) := F(x) + J(x) [y-x].
201  // Newton iteration for F(y)=b is given by L[x_old](x_new) = b, i.e.
202  // F(x_old) + J(x_old) [x_new-x_old] = b.
203  //
204  // To minimize: 1/2 |F(y)-b|^2, subject to: l(y) >= 0, we may consider the
205  // iteration: minimize: |L[x_old](x_new)-b|^2, subject to l(x_new) >= 0,
206  // i.e. minimize: |F(x_old) + J(x_old) [x_new-x_old] - b|^2.
207 
208  // This method uses:
209  // Newton iteration: x := x + J(x)^{-1} [pt-F(x)]
210  // or when dim != sdim: x := x + [J^t.J]^{-1}.J^t [pt-F(x)]
211 
212  // Compute the physical coordinates of the current point:
213  T->Transform(xip, y);
214  if (print_level >= 3)
215  {
216  NewtonPrint(11, 0.); // continuation line
217  NewtonPrintPoint("approx_pt", y, ", ");
218  NewtonPrintPoint("exact_pt", pt, "\n");
219  }
220  subtract(pt, y, y); // y = pt-y
221 
222  // Check for convergence in physical coordinates:
223  err_phys = y.Normlinf();
224  if (err_phys < phys_tol)
225  {
226  if (print_level >= 1)
227  {
228  NewtonPrint(1, (double)it);
229  NewtonPrint(3, dx_norm);
230  NewtonPrint(30, err_phys);
231  }
232  ip = xip;
233  if (solver_type != Newton) { return Inside; }
234  return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside;
235  }
236  if (print_level >= 1)
237  {
238  if (it == 0 || print_level >= 2)
239  {
240  NewtonPrint(1, (double)it);
241  NewtonPrint(3, dx_norm);
242  NewtonPrint(18, err_phys);
243  }
244  }
245 
246  if (hit_bdr)
247  {
248  xip.Get(xd, dim); // xip -> x
249  if (prev_hit_bdr || it == max_iter || print_level >= 2)
250  {
251  prev_xip.Get(dxd, dim); // prev_xip -> dx
252  subtract(x, dx, dx); // dx = xip - prev_xip
253  real_dx_norm = dx.Normlinf();
254  if (print_level >= 2)
255  {
256  NewtonPrint(41, real_dx_norm);
257  }
258  if (prev_hit_bdr && real_dx_norm < ref_tol)
259  {
260  if (print_level >= 0)
261  {
262  if (print_level <= 1)
263  {
264  NewtonPrint(1, (double)it);
265  NewtonPrint(3, dx_norm);
266  NewtonPrint(18, err_phys);
267  NewtonPrint(41, real_dx_norm);
268  }
269  mfem::out << "Newton: *** stuck on boundary!\n";
270  }
271  return Outside;
272  }
273  }
274  }
275 
276  if (it == max_iter) { break; }
277 
278  // Perform a Newton step:
279  T->SetIntPoint(&xip);
280  T->InverseJacobian().Mult(y, dx);
281  x += dx;
282  it++;
283  if (solver_type != Newton)
284  {
285  prev_xip = xip;
286  prev_hit_bdr = hit_bdr;
287  }
288  xip.Set(xd, dim); // x -> xip
289 
290  // Perform projection based on solver_type:
291  switch (solver_type)
292  {
293  case Newton: break;
295  hit_bdr = !Geometry::ProjectPoint(geom, prev_xip, xip); break;
297  hit_bdr = !Geometry::ProjectPoint(geom, xip); break;
298  default: MFEM_ABORT("invalid solver type");
299  }
300  if (print_level >= 3)
301  {
302  NewtonPrint(1, double(it));
303  xip.Get(xd, dim); // xip -> x
304  NewtonPrintPoint(", ref_pt", x, "\n");
305  }
306 
307  // Check for convergence in reference coordinates:
308  dx_norm = dx.Normlinf();
309  if (dx_norm < ref_tol)
310  {
311  if (print_level >= 1)
312  {
313  NewtonPrint(1, (double)it);
314  NewtonPrint(27, dx_norm);
315  }
316  ip = xip;
317  if (solver_type != Newton) { return Inside; }
318  return Geometry::CheckPoint(geom, ip, ip_tol) ? Inside : Outside;
319  }
320  }
321  if (print_level >= 0)
322  {
323  if (print_level <= 1)
324  {
325  NewtonPrint(1, (double)max_iter);
326  NewtonPrint(3, dx_norm);
327  NewtonPrint(18, err_phys);
328  if (hit_bdr) { NewtonPrint(41, real_dx_norm); }
329  }
330  mfem::out << "Newton: *** iteration did not converge!\n";
331  }
332  ip = xip;
333  return Unknown;
334 }
335 
337  IntegrationPoint &ip)
338 {
339  MFEM_VERIFY(T != NULL, "invalid ElementTransformation");
340 
341  // Select initial guess ...
342  switch (init_guess_type)
343  {
344  case Center:
346  break;
347 
348  case ClosestPhysNode:
349  case ClosestRefNode:
350  {
351  const int order = std::max(T->Order()+rel_qpts_order, 0);
352  if (order == 0)
353  {
355  }
356  else
357  {
358  const int old_type = GlobGeometryRefiner.GetType();
360  RefinedGeometry &RefG =
362  int closest_idx = (init_guess_type == ClosestPhysNode) ?
363  FindClosestPhysPoint(pt, RefG.RefPts) :
364  FindClosestRefPoint(pt, RefG.RefPts);
365  ip0 = &RefG.RefPts.IntPoint(closest_idx);
366  GlobGeometryRefiner.SetType(old_type);
367  }
368  break;
369  }
370 
371  case GivenPoint:
372  break;
373 
374  default:
375  MFEM_ABORT("invalid initial guess type");
376  }
377 
378  // Call the solver ...
379  return NewtonSolve(pt, ip);
380 }
381 
382 
384  Geometry::Type GeomType)
385 {
386  switch (GeomType)
387  {
388  case Geometry::POINT : FElem = &PointFE; break;
389  case Geometry::SEGMENT : FElem = &SegmentFE; break;
390  case Geometry::TRIANGLE : FElem = &TriangleFE; break;
391  case Geometry::SQUARE : FElem = &QuadrilateralFE; break;
392  case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
393  case Geometry::CUBE : FElem = &HexahedronFE; break;
394  case Geometry::PRISM : FElem = &WedgeFE; break;
395  case Geometry::PYRAMID : FElem = &PyramidFE; break;
396  default:
397  MFEM_ABORT("unknown Geometry::Type!");
398  }
399  int dim = FElem->GetDim();
400  int dof = FElem->GetDof();
401  const IntegrationRule &nodes = FElem->GetNodes();
402  PointMat.SetSize(dim, dof);
403  for (int j = 0; j < dof; j++)
404  {
405  nodes.IntPoint(j).Get(&PointMat(0,j), dim);
406  }
407  geom = GeomType;
408 }
409 
410 const DenseMatrix &IsoparametricTransformation::EvalJacobian()
411 {
412  MFEM_ASSERT((EvalState & JACOBIAN_MASK) == 0, "");
413 
414  dshape.SetSize(FElem->GetDof(), FElem->GetDim());
415  dFdx.SetSize(PointMat.Height(), dshape.Width());
416  if (dshape.Width() > 0)
417  {
418  FElem->CalcDShape(*IntPoint, dshape);
419  Mult(PointMat, dshape, dFdx);
420  }
422 
423  return dFdx;
424 }
425 
426 const DenseMatrix &IsoparametricTransformation::EvalHessian()
427 {
428  MFEM_ASSERT((EvalState & HESSIAN_MASK) == 0, "");
429 
430  int Dim = FElem->GetDim();
431  d2shape.SetSize(FElem->GetDof(), (Dim*(Dim+1))/2);
432  d2Fdx2.SetSize(PointMat.Height(), d2shape.Width());
433  if (d2shape.Width() > 0)
434  {
435  FElem->CalcHessian(*IntPoint, d2shape);
436  Mult(PointMat, d2shape, d2Fdx2);
437  }
439 
440  return d2Fdx2;
441 }
442 
444 {
445  switch (FElem->Space())
446  {
447  case FunctionSpace::Pk:
448  return (FElem->GetOrder()-1);
449  case FunctionSpace::Qk:
450  return (FElem->GetOrder());
451  default:
452  MFEM_ABORT("unsupported finite element");
453  }
454  return 0;
455 }
456 
458 {
459  switch (FElem->Space())
460  {
461  case FunctionSpace::Pk:
462  return (FElem->GetOrder() - 1) * FElem->GetDim();
463  case FunctionSpace::Qk:
464  return (FElem->GetOrder() * FElem->GetDim() - 1);
465  default:
466  MFEM_ABORT("unsupported finite element");
467  }
468  return 0;
469 }
470 
472 {
473  if (FElem->Space() == fe->Space())
474  {
475  int k = FElem->GetOrder();
476  int d = FElem->GetDim();
477  int l = fe->GetOrder();
478  switch (fe->Space())
479  {
480  case FunctionSpace::Pk:
481  return ((k-1)*(d-1)+(l-1));
482  case FunctionSpace::Qk:
483  return (k*(d-1)+(l-1));
484  default:
485  MFEM_ABORT("unsupported finite element");
486  }
487  }
488  MFEM_ABORT("incompatible finite elements");
489  return 0;
490 }
491 
493  Vector &trans)
494 {
495  shape.SetSize(FElem->GetDof());
496  trans.SetSize(PointMat.Height());
497 
498  FElem -> CalcShape(ip, shape);
499  PointMat.Mult(shape, trans);
500 }
501 
503  DenseMatrix &tr)
504 {
505  int dof, n, dim, i, j, k;
506 
507  dim = PointMat.Height();
508  dof = FElem->GetDof();
509  n = ir.GetNPoints();
510 
511  shape.SetSize(dof);
512  tr.SetSize(dim, n);
513 
514  for (j = 0; j < n; j++)
515  {
516  FElem -> CalcShape (ir.IntPoint(j), shape);
517  for (i = 0; i < dim; i++)
518  {
519  tr(i, j) = 0.0;
520  for (k = 0; k < dof; k++)
521  {
522  tr(i, j) += PointMat(i, k) * shape(k);
523  }
524  }
525  }
526 }
527 
529  DenseMatrix &result)
530 {
531  MFEM_ASSERT(matrix.Height() == GetDimension(), "invalid input");
532  result.SetSize(PointMat.Height(), matrix.Width());
533 
534  IntegrationPoint ip;
535  Vector col;
536 
537  for (int j = 0; j < matrix.Width(); j++)
538  {
539  ip.Set(matrix.GetColumn(j), matrix.Height());
540 
541  result.GetColumnReference(j, col);
542  Transform(ip, col);
543  }
544 }
545 
547  IntegrationPoint &ip2)
548 {
549  double vec[3];
550  Vector v (vec, Transf.GetPointMat().Height());
551 
552  Transf.Transform (ip1, v);
553  ip2.Set(vec, v.Size());
554 }
555 
557  IntegrationRule &ir2)
558 {
559  int i, n;
560 
561  n = ir1.GetNPoints();
562  for (i = 0; i < n; i++)
563  {
564  Transform (ir1.IntPoint(i), ir2.IntPoint(i));
565  }
566 }
567 
569 {
571 
572  if (mask & 4)
573  {
574  Loc1.Transform(*face_ip, eip1);
575  if (Elem1)
576  {
577  Elem1->SetIntPoint(&eip1);
578  }
579  }
580  if (mask & 8)
581  {
582  Loc2.Transform(*face_ip, eip2);
583  if (Elem2)
584  {
585  Elem2->SetIntPoint(&eip2);
586  }
587  }
588 }
589 
592 {
593  MFEM_VERIFY(mask & HAVE_ELEM1 && Elem1 != NULL, "The ElementTransformation "
594  "for the element has not been configured for side 1.");
595  return *Elem1;
596 }
597 
600 {
601  MFEM_VERIFY(mask & HAVE_ELEM2 && Elem2 != NULL, "The ElementTransformation "
602  "for the element has not been configured for side 2.");
603  return *Elem2;
604 }
605 
608 {
609  MFEM_VERIFY(mask & HAVE_LOC1, "The IntegrationPointTransformation "
610  "for the element has not been configured for side 1.");
611  return Loc1;
612 }
613 
616 {
617  MFEM_VERIFY(mask & HAVE_LOC2, "The IntegrationPointTransformation "
618  "for the element has not been configured for side 2.");
619  return Loc2;
620 }
621 
623  Vector &trans)
624 {
625  MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
626  "for the face has not been configured.");
628 }
629 
631  DenseMatrix &tr)
632 {
633  MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
634  "for the face has not been configured.");
636 }
637 
639  DenseMatrix &result)
640 {
641  MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
642  "for the face has not been configured.");
644 }
645 
647  std::ostream &os)
648 {
649  // Check that the face vertices are mapped to the same physical location
650  // when using the following three transformations:
651  // - the face transformation, *this
652  // - Loc1 + Elem1
653  // - Loc2 + Elem2, if present.
654 
655  const bool have_face = (mask & 16);
656  const bool have_el1 = (mask & 1) && (mask & 4);
657  const bool have_el2 = (mask & 2) && (mask & 8) && (Elem2No >= 0);
658  if (int(have_face) + int(have_el1) + int(have_el2) < 2)
659  {
660  // need at least two different transformations to perform a check
661  return 0.0;
662  }
663 
665 
666  double max_dist = 0.0;
667  Vector dist(v_ir.GetNPoints());
668  DenseMatrix coords_base, coords_el;
669  IntegrationRule v_eir(v_ir.GetNPoints());
670  if (have_face)
671  {
672  Transform(v_ir, coords_base);
673  if (print_level > 0)
674  {
675  os << "\nface vertex coordinates (from face transform):\n"
676  << "----------------------------------------------\n";
677  coords_base.PrintT(os, coords_base.Height());
678  }
679  }
680  if (have_el1)
681  {
682  Loc1.Transform(v_ir, v_eir);
683  Elem1->Transform(v_eir, coords_el);
684  if (print_level > 0)
685  {
686  os << "\nface vertex coordinates (from element 1 transform):\n"
687  << "---------------------------------------------------\n";
688  coords_el.PrintT(os, coords_el.Height());
689  }
690  if (have_face)
691  {
692  coords_el -= coords_base;
693  coords_el.Norm2(dist);
694  max_dist = std::max(max_dist, dist.Normlinf());
695  }
696  else
697  {
698  coords_base = coords_el;
699  }
700  }
701  if (have_el2)
702  {
703  Loc2.Transform(v_ir, v_eir);
704  Elem2->Transform(v_eir, coords_el);
705  if (print_level > 0)
706  {
707  os << "\nface vertex coordinates (from element 2 transform):\n"
708  << "---------------------------------------------------\n";
709  coords_el.PrintT(os, coords_el.Height());
710  }
711  coords_el -= coords_base;
712  coords_el.Norm2(dist);
713  max_dist = std::max(max_dist, dist.Normlinf());
714  }
715 
716  return max_dist;
717 }
718 
719 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Abstract class for all finite elements.
Definition: fe_base.hpp:235
IntegrationPointTransformation & GetIntPoint1Transformation()
Definition: eltrans.cpp:607
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
int NewtonSolve(const Vector &pt, IntegrationPoint &ip)
Definition: eltrans.cpp:172
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:135
void Get(double *p, const int dim) const
Definition: intrules.hpp:51
Tensor products of polynomials of order k.
Definition: fe_base.hpp:224
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:311
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
ElementTransformation & GetElement2Transformation()
Definition: eltrans.cpp:599
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void SetIntPoint(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.cpp:568
Element on side 1 is configured.
Definition: eltrans.hpp:513
Use a specific point, set with SetInitialGuess().
Definition: eltrans.hpp:200
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
const IntegrationPoint * IntPoint
Definition: eltrans.hpp:26
ElementTransformation * T
Definition: eltrans.hpp:231
class LinearPyramidFiniteElement PyramidFE
Definition: fe.cpp:44
virtual void Transform(const IntegrationPoint &, Vector &)
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition: eltrans.cpp:622
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
Use the center of the reference element.
Definition: eltrans.hpp:191
virtual int Transform(const Vector &pt, IntegrationPoint &ip)
Given a point, pt, in physical space, find its reference coordinates, ip.
Definition: eltrans.cpp:336
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:73
int GetDimension() const
Return the topological dimension of the reference element.
Definition: eltrans.hpp:165
IntegrationPointTransformation & GetIntPoint2Transformation()
Definition: eltrans.cpp:615
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2440
void NewtonPrint(int mode, double val)
Definition: eltrans.cpp:129
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
int Size() const
Returns the size of the vector.
Definition: vector.hpp:200
Point transformation for side 2 is configured.
Definition: eltrans.hpp:516
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:524
const DenseMatrix & EvalAdjugateJ()
Definition: eltrans.cpp:36
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:337
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:145
double Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:830
virtual void Transform(const IntegrationPoint &, Vector &)
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition: eltrans.cpp:492
Element on side 2 is configured.
Definition: eltrans.hpp:514
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
virtual int OrderJ() const
Return the order of the elements of the Jacobian of the transformation.
Definition: eltrans.cpp:443
ElementTransformation * Elem2
Definition: eltrans.hpp:522
void NewtonPrintPoint(const char *prefix, const Vector &pt, const char *suffix)
Definition: eltrans.cpp:158
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:265
Polynomials of order k.
Definition: fe_base.hpp:223
virtual int OrderW() const
Return the order of the determinant of the Jacobian (weight) of the transformation.
Definition: eltrans.cpp:457
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:97
double Weight() const
Definition: densemat.cpp:506
PointFiniteElement PointFE
Definition: point.cpp:30
Geometry Geometries
Definition: fe.cpp:49
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2512
The point is inside the element.
Definition: eltrans.hpp:224
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:67
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
Definition: geom.hpp:306
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:524
The point is probably outside the element.
Definition: eltrans.hpp:225
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:119
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:390
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1773
IntegrationRule RefPts
Definition: geom.hpp:282
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1099
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1292
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:654
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2548
Point transformation for side 1 is configured.
Definition: eltrans.hpp:515
const DenseMatrix & EvalInverseJ()
Definition: eltrans.cpp:57
virtual int OrderGrad(const FiniteElement *fe) const
Return the order of .
Definition: eltrans.cpp:471
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:583
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1381
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:383
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
double CheckConsistency(int print_level=0, std::ostream &out=mfem::out)
Check for self-consistency: compares the result of mapping the reference face vertices to physical co...
Definition: eltrans.cpp:646
virtual int Order() const =0
Return the order of the current element we are using for the transformation.
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:460
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_base.cpp:100
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:289
class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:36
Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
Face transformation is configured.
Definition: eltrans.hpp:517
Class for integration point with weight.
Definition: intrules.hpp:25
The algorithm failed to determine where the point is.
Definition: eltrans.hpp:226
int GetType() const
Get the Quadrature1D type of points used for subdivision.
Definition: geom.hpp:308
class LinearWedgeFiniteElement WedgeFE
Definition: fe.cpp:40
int dim
Definition: ex24.cpp:53
int FindClosestPhysPoint(const Vector &pt, const IntegrationRule &ir)
Find the IntegrationPoint mapped closest to pt.
Definition: eltrans.cpp:70
ElementTransformation & GetElement1Transformation()
Definition: eltrans.cpp:591
IsoparametricTransformation Transf
Definition: eltrans.hpp:464
ElementTransformation * Elem1
Definition: eltrans.hpp:522
BiLinear2DFiniteElement QuadrilateralFE
Vector data type.
Definition: vector.hpp:60
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
const DenseMatrix & EvalTransAdjugateJ()
Definition: eltrans.cpp:46
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:405
const IntegrationPoint * ip0
Definition: eltrans.hpp:234
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:52
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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:66
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:404
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...
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49