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