MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
eltrans.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
16namespace 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();
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 real_t minDist = std::numeric_limits<real_t>::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 real_t 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 real_t minDist = std::numeric_limits<real_t>::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 real_t 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
174{
175 MFEM_ASSERT(pt.Size() == T->GetSpaceDim(), "invalid point");
176
177 const real_t 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 real_t 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, (real_t)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, (real_t)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, (real_t)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, real_t(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, (real_t)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 {
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
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 RefinedGeometry &RefG = *refiner.Refine(T->GetGeometryType(), order);
359 int closest_idx = (init_guess_type == ClosestPhysNode) ?
360 FindClosestPhysPoint(pt, RefG.RefPts) :
361 FindClosestRefPoint(pt, RefG.RefPts);
362 ip0 = &RefG.RefPts.IntPoint(closest_idx);
363 }
364 break;
365 }
366
367 case GivenPoint:
368 break;
369
370 default:
371 MFEM_ABORT("invalid initial guess type");
372 }
373
374 // Call the solver ...
375 return NewtonSolve(pt, ip);
376}
377
378
380 Geometry::Type GeomType)
381{
382 switch (GeomType)
383 {
384 case Geometry::POINT : FElem = &PointFE; break;
385 case Geometry::SEGMENT : FElem = &SegmentFE; break;
386 case Geometry::TRIANGLE : FElem = &TriangleFE; break;
387 case Geometry::SQUARE : FElem = &QuadrilateralFE; break;
388 case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
389 case Geometry::CUBE : FElem = &HexahedronFE; break;
390 case Geometry::PRISM : FElem = &WedgeFE; break;
391 case Geometry::PYRAMID : FElem = &PyramidFE; break;
392 default:
393 MFEM_ABORT("unknown Geometry::Type!");
394 }
395 int dim = FElem->GetDim();
396 int dof = FElem->GetDof();
397 const IntegrationRule &nodes = FElem->GetNodes();
398 PointMat.SetSize(dim, dof);
399 for (int j = 0; j < dof; j++)
400 {
401 nodes.IntPoint(j).Get(&PointMat(0,j), dim);
402 }
403 geom = GeomType;
404}
405
406const DenseMatrix &IsoparametricTransformation::EvalJacobian()
407{
408 MFEM_ASSERT((EvalState & JACOBIAN_MASK) == 0, "");
409
410 dshape.SetSize(FElem->GetDof(), FElem->GetDim());
411 dFdx.SetSize(PointMat.Height(), dshape.Width());
412 if (dshape.Width() > 0)
413 {
414 FElem->CalcDShape(*IntPoint, dshape);
415 Mult(PointMat, dshape, dFdx);
416 }
418
419 return dFdx;
420}
421
422const DenseMatrix &IsoparametricTransformation::EvalHessian()
423{
424 MFEM_ASSERT((EvalState & HESSIAN_MASK) == 0, "");
425
426 int Dim = FElem->GetDim();
427 d2shape.SetSize(FElem->GetDof(), (Dim*(Dim+1))/2);
428 d2Fdx2.SetSize(PointMat.Height(), d2shape.Width());
429 if (d2shape.Width() > 0)
430 {
431 FElem->CalcHessian(*IntPoint, d2shape);
432 Mult(PointMat, d2shape, d2Fdx2);
433 }
435
436 return d2Fdx2;
437}
438
440{
441 switch (FElem->Space())
442 {
444 return (FElem->GetOrder()-1);
446 return (FElem->GetOrder());
447 default:
448 MFEM_ABORT("unsupported finite element");
449 }
450 return 0;
451}
452
454{
455 switch (FElem->Space())
456 {
458 return (FElem->GetOrder() - 1) * FElem->GetDim();
460 return (FElem->GetOrder() * FElem->GetDim() - 1);
461 default:
462 MFEM_ABORT("unsupported finite element");
463 }
464 return 0;
465}
466
468{
469 if (FElem->Space() == fe->Space())
470 {
471 int k = FElem->GetOrder();
472 int d = FElem->GetDim();
473 int l = fe->GetOrder();
474 switch (fe->Space())
475 {
477 return ((k-1)*(d-1)+(l-1));
479 return (k*(d-1)+(l-1));
480 default:
481 MFEM_ABORT("unsupported finite element");
482 }
483 }
484 MFEM_ABORT("incompatible finite elements");
485 return 0;
486}
487
489 Vector &trans)
490{
491 MFEM_ASSERT(FElem != nullptr, "Must provide a valid FiniteElement object!");
492 shape.SetSize(FElem->GetDof());
493 trans.SetSize(PointMat.Height());
494
495 FElem -> CalcShape(ip, shape);
496 PointMat.Mult(shape, trans);
497}
498
500 DenseMatrix &tr)
501{
502 int dof, n, dim, i, j, k;
503
504 dim = PointMat.Height();
505 dof = FElem->GetDof();
506 n = ir.GetNPoints();
507
508 shape.SetSize(dof);
509 tr.SetSize(dim, n);
510
511 for (j = 0; j < n; j++)
512 {
513 FElem -> CalcShape (ir.IntPoint(j), shape);
514 for (i = 0; i < dim; i++)
515 {
516 tr(i, j) = 0.0;
517 for (k = 0; k < dof; k++)
518 {
519 tr(i, j) += PointMat(i, k) * shape(k);
520 }
521 }
522 }
523}
524
526 DenseMatrix &result)
527{
528 MFEM_ASSERT(matrix.Height() == GetDimension(), "invalid input");
529 result.SetSize(PointMat.Height(), matrix.Width());
530
532 Vector col;
533
534 for (int j = 0; j < matrix.Width(); j++)
535 {
536 ip.Set(matrix.GetColumn(j), matrix.Height());
537
538 result.GetColumnReference(j, col);
539 Transform(ip, col);
540 }
541}
542
544 IntegrationPoint &ip2)
545{
546 real_t vec[3];
547 Vector v (vec, Transf.GetPointMat().Height());
548
549 Transf.Transform (ip1, v);
550 ip2.Set(vec, v.Size());
551}
552
554 IntegrationRule &ir2)
555{
556 int i, n;
557
558 n = ir1.GetNPoints();
559 for (i = 0; i < n; i++)
560 {
561 Transform (ir1.IntPoint(i), ir2.IntPoint(i));
562 }
563}
564
566{
568
569 if (mask & 4)
570 {
571 Loc1.Transform(*face_ip, eip1);
572 if (Elem1)
573 {
574 Elem1->SetIntPoint(&eip1);
575 }
576 }
577 if (mask & 8)
578 {
579 Loc2.Transform(*face_ip, eip2);
580 if (Elem2)
581 {
582 Elem2->SetIntPoint(&eip2);
583 }
584 }
585}
586
589{
590 MFEM_VERIFY(mask & HAVE_ELEM1 && Elem1 != NULL, "The ElementTransformation "
591 "for the element has not been configured for side 1.");
592 return *Elem1;
593}
594
597{
598 MFEM_VERIFY(mask & HAVE_ELEM2 && Elem2 != NULL, "The ElementTransformation "
599 "for the element has not been configured for side 2.");
600 return *Elem2;
601}
602
605{
606 MFEM_VERIFY(mask & HAVE_LOC1, "The IntegrationPointTransformation "
607 "for the element has not been configured for side 1.");
608 return Loc1;
609}
610
613{
614 MFEM_VERIFY(mask & HAVE_LOC2, "The IntegrationPointTransformation "
615 "for the element has not been configured for side 2.");
616 return Loc2;
617}
618
620 Vector &trans)
621{
622 MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
623 "for the face has not been configured.");
625}
626
628 DenseMatrix &tr)
629{
630 MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
631 "for the face has not been configured.");
633}
634
636 DenseMatrix &result)
637{
638 MFEM_VERIFY(mask & HAVE_FACE, "The ElementTransformation "
639 "for the face has not been configured.");
641}
642
644 std::ostream &os)
645{
646 // Check that the face vertices are mapped to the same physical location
647 // when using the following three transformations:
648 // - the face transformation, *this
649 // - Loc1 + Elem1
650 // - Loc2 + Elem2, if present.
651
652 const bool have_face = (mask & 16);
653 const bool have_el1 = (mask & 1) && (mask & 4);
654 const bool have_el2 = (mask & 2) && (mask & 8) && (Elem2No >= 0);
655 if (int(have_face) + int(have_el1) + int(have_el2) < 2)
656 {
657 // need at least two different transformations to perform a check
658 return 0.0;
659 }
660
662
663 real_t max_dist = 0.0;
664 Vector dist(v_ir.GetNPoints());
665 DenseMatrix coords_base, coords_el;
666 IntegrationRule v_eir(v_ir.GetNPoints());
667 if (have_face)
668 {
669 Transform(v_ir, coords_base);
670 if (print_level > 0)
671 {
672 os << "\nface vertex coordinates (from face transform):\n"
673 << "----------------------------------------------\n";
674 coords_base.PrintT(os, coords_base.Height());
675 }
676 }
677 if (have_el1)
678 {
679 Loc1.Transform(v_ir, v_eir);
680 Elem1->Transform(v_eir, coords_el);
681 if (print_level > 0)
682 {
683 os << "\nface vertex coordinates (from element 1 transform):\n"
684 << "---------------------------------------------------\n";
685 coords_el.PrintT(os, coords_el.Height());
686 }
687 if (have_face)
688 {
689 coords_el -= coords_base;
690 coords_el.Norm2(dist);
691 max_dist = std::max(max_dist, dist.Normlinf());
692 }
693 else
694 {
695 coords_base = coords_el;
696 }
697 }
698 if (have_el2)
699 {
700 Loc2.Transform(v_ir, v_eir);
701 Elem2->Transform(v_eir, coords_el);
702 if (print_level > 0)
703 {
704 os << "\nface vertex coordinates (from element 2 transform):\n"
705 << "---------------------------------------------------\n";
706 coords_el.PrintT(os, coords_el.Height());
707 }
708 coords_el -= coords_base;
709 coords_el.Norm2(dist);
710 max_dist = std::max(max_dist, dist.Normlinf());
711 }
712
713 return max_dist;
714}
715
716}
Data type dense matrix using column-major storage.
Definition: densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition: densemat.cpp:220
void Transpose()
(*this) = (*this)^t
Definition: densemat.cpp:1532
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
real_t Weight() const
Definition: densemat.cpp:592
virtual void PrintT(std::ostream &out=mfem::out, int width_=4) const
Prints the transpose matrix to stream out.
Definition: densemat.cpp:2379
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1443
void Norm2(real_t *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:898
virtual int Order() const =0
Return the order of the current element we are using for the transformation.
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:145
const DenseMatrix & AdjugateJacobian()
Return the adjugate of the Jacobian matrix of the transformation at the currently set IntegrationPoin...
Definition: eltrans.hpp:135
const DenseMatrix & EvalAdjugateJ()
Definition: eltrans.cpp:36
virtual int GetSpaceDim() const =0
Get the dimension of the target (physical) space.
const IntegrationPoint * IntPoint
Definition: eltrans.hpp:26
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition: eltrans.hpp:119
int GetDimension() const
Return the topological dimension of the reference element.
Definition: eltrans.hpp:165
const DenseMatrix & EvalInverseJ()
Definition: eltrans.cpp:57
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
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 ...
ElementTransformation * Elem2
Definition: eltrans.hpp:525
void SetIntPoint(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.cpp:565
ElementTransformation * Elem1
Definition: eltrans.hpp:525
@ HAVE_ELEM2
Element on side 2 is configured.
Definition: eltrans.hpp:517
@ HAVE_LOC1
Point transformation for side 1 is configured.
Definition: eltrans.hpp:518
@ HAVE_ELEM1
Element on side 1 is configured.
Definition: eltrans.hpp:516
@ HAVE_FACE
Face transformation is configured.
Definition: eltrans.hpp:520
@ HAVE_LOC2
Point transformation for side 2 is configured.
Definition: eltrans.hpp:519
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:527
IntegrationPointTransformation & GetIntPoint1Transformation()
Definition: eltrans.cpp:604
ElementTransformation & GetElement1Transformation()
Definition: eltrans.cpp:588
virtual void Transform(const IntegrationPoint &, Vector &)
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition: eltrans.cpp:619
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:527
IntegrationPointTransformation & GetIntPoint2Transformation()
Definition: eltrans.cpp:612
real_t 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:643
ElementTransformation & GetElement2Transformation()
Definition: eltrans.cpp:596
Abstract class for all finite elements.
Definition: fe_base.hpp:239
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:101
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:333
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:316
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe_base.hpp:395
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...
int Space() const
Returns the type of FunctionSpace on the element.
Definition: fe_base.hpp:343
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:329
@ Pk
Polynomials of order k.
Definition: fe_base.hpp:226
@ Qk
Tensor products of polynomials of order k.
Definition: fe_base.hpp:227
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:1136
const IntegrationPoint & GetCenter(int GeomType) const
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
Definition: geom.cpp:293
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:435
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
Definition: geom.cpp:616
IsoparametricTransformation Transf
Definition: eltrans.hpp:467
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:543
Class for integration point with weight.
Definition: intrules.hpp:35
void Get(real_t *p, const int dim) const
Definition: intrules.hpp:60
void Set(const real_t *p, const int dim)
Definition: intrules.hpp:46
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:259
virtual int Transform(const Vector &pt, IntegrationPoint &ip)
Given a point, pt, in physical space, find its reference coordinates, ip.
Definition: eltrans.cpp:336
@ GivenPoint
Use a specific point, set with SetInitialGuess().
Definition: eltrans.hpp:201
@ Center
Use the center of the reference element.
Definition: eltrans.hpp:192
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
const IntegrationPoint * ip0
Definition: eltrans.hpp:235
void NewtonPrintPoint(const char *prefix, const Vector &pt, const char *suffix)
Definition: eltrans.cpp:158
@ Outside
The point is probably outside the element.
Definition: eltrans.hpp:226
@ Unknown
The algorithm failed to determine where the point is.
Definition: eltrans.hpp:227
@ Inside
The point is inside the element.
Definition: eltrans.hpp:225
int FindClosestPhysPoint(const Vector &pt, const IntegrationRule &ir)
Find the IntegrationPoint mapped closest to pt.
Definition: eltrans.cpp:70
int NewtonSolve(const Vector &pt, IntegrationPoint &ip)
Definition: eltrans.cpp:172
void NewtonPrint(int mode, real_t val)
Definition: eltrans.cpp:129
ElementTransformation * T
Definition: eltrans.hpp:232
virtual int OrderJ() const
Return the order of the elements of the Jacobian of the transformation.
Definition: eltrans.cpp:439
virtual int OrderGrad(const FiniteElement *fe) const
Return the order of .
Definition: eltrans.cpp:467
void SetIdentityTransformation(Geometry::Type GeomType)
Set the FiniteElement Geometry for the reference elements being used.
Definition: eltrans.cpp:379
const DenseMatrix & GetPointMat() const
Return the stored point matrix.
Definition: eltrans.hpp:405
virtual int OrderW() const
Return the order of the determinant of the Jacobian (weight) of the transformation.
Definition: eltrans.cpp:453
virtual void Transform(const IntegrationPoint &, Vector &)
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition: eltrans.cpp:488
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
IntegrationRule RefPts
Definition: geom.hpp:317
Vector data type.
Definition: vector.hpp:80
real_t Normlinf() const
Returns the l_infinity norm of the vector.
Definition: vector.cpp:850
real_t Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:832
int Size() const
Returns the size of the vector.
Definition: vector.hpp:218
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:538
real_t DistanceTo(const real_t *p) const
Compute the Euclidean distance to another vector.
Definition: vector.hpp:690
int dim
Definition: ex24.cpp:53
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:52
PointFiniteElement PointFE
Definition: point.cpp:42
TriLinear3DFiniteElement HexahedronFE
Definition: hexahedron.cpp:60
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
Definition: densemat.cpp:2679
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
Geometry Geometries
Definition: fe.cpp:49
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2715
MFEM_EXPORT class Linear3DFiniteElement TetrahedronFE
Definition: fe.cpp:36
MFEM_EXPORT class LinearWedgeFiniteElement WedgeFE
Definition: fe.cpp:40
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
Definition: densemat.cpp:2607
BiLinear2DFiniteElement QuadrilateralFE
void subtract(const Vector &x, const Vector &y, Vector &z)
Definition: vector.cpp:472
class LinearPyramidFiniteElement PyramidFE
Definition: fe.cpp:44
float real_t
Definition: config.hpp:43
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32