12 #include "../mesh/mesh_headers.hpp"
73 MFEM_VERIFY(
T != NULL,
"invalid ElementTransformation");
81 double minDist = std::numeric_limits<double>::max();
85 for (
int i = 0; i < npts; ++i)
87 double dist = pt.
DistanceTo(physPts.GetColumn(i));
100 MFEM_VERIFY(
T != NULL,
"invalid ElementTransformation");
105 double minDist = std::numeric_limits<double>::max();
111 for (
int i = 0; i < npts; ++i)
118 double dist = dr.Norml2();
136 case 0: os <<
", ";
break;
137 case 1: os <<
"Newton: ";
break;
138 case 2: os <<
" ";
break;
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;
152 case 1: os <<
'\n';
break;
153 case 2: os <<
" (converged)\n";
break;
154 case 3: os <<
" (actual)\n";
break;
164 os << prefix <<
" = (";
165 for (
int j = 0; j < pt.
Size(); j++)
167 os << (j > 0 ?
", " :
"") << pt(j);
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;
196 for (
int it = 0;
true; )
223 err_phys = y.Normlinf();
224 if (err_phys < phys_tol)
251 prev_xip.
Get(dxd, dim);
258 if (prev_hit_bdr && real_dx_norm <
ref_tol)
269 mfem::out <<
"Newton: *** stuck on boundary!\n";
286 prev_hit_bdr = hit_bdr;
298 default: MFEM_ABORT(
"invalid solver type");
330 mfem::out <<
"Newton: *** iteration did not converge!\n";
339 MFEM_VERIFY(
T != NULL,
"invalid ElementTransformation");
375 MFEM_ABORT(
"invalid initial guess type");
397 MFEM_ABORT(
"unknown Geometry::Type!");
400 int dof = FElem->
GetDof();
403 for (
int j = 0; j < dof; j++)
410 const DenseMatrix &IsoparametricTransformation::EvalJacobian()
416 if (dshape.
Width() > 0)
426 const DenseMatrix &IsoparametricTransformation::EvalHessian()
430 int Dim = FElem->
GetDim();
433 if (d2shape.
Width() > 0)
445 switch (FElem->
Space())
452 MFEM_ABORT(
"unsupported finite element");
459 switch (FElem->
Space())
466 MFEM_ABORT(
"unsupported finite element");
481 return ((k-1)*(d-1)+(l-1));
483 return (k*(d-1)+(l-1));
485 MFEM_ABORT(
"unsupported finite element");
488 MFEM_ABORT(
"incompatible finite elements");
498 FElem -> CalcShape(ip, shape);
499 PointMat.
Mult(shape, trans);
505 int dof, n,
dim, i, j, k;
514 for (j = 0; j < n; j++)
516 FElem -> CalcShape (ir.
IntPoint(j), shape);
517 for (i = 0; i <
dim; i++)
520 for (k = 0; k < dof; k++)
522 tr(i, j) += PointMat(i, k) * shape(k);
537 for (
int j = 0; j < matrix.
Width(); j++)
553 ip2.
Set(vec, v.Size());
562 for (i = 0; i < n; i++)
593 MFEM_VERIFY(mask &
HAVE_ELEM1 &&
Elem1 != NULL,
"The ElementTransformation "
594 "for the element has not been configured for side 1.");
601 MFEM_VERIFY(mask &
HAVE_ELEM2 &&
Elem2 != NULL,
"The ElementTransformation "
602 "for the element has not been configured for side 2.");
609 MFEM_VERIFY(mask &
HAVE_LOC1,
"The IntegrationPointTransformation "
610 "for the element has not been configured for side 1.");
617 MFEM_VERIFY(mask &
HAVE_LOC2,
"The IntegrationPointTransformation "
618 "for the element has not been configured for side 2.");
625 MFEM_VERIFY(mask &
HAVE_FACE,
"The ElementTransformation "
626 "for the face has not been configured.");
633 MFEM_VERIFY(mask &
HAVE_FACE,
"The ElementTransformation "
634 "for the face has not been configured.");
641 MFEM_VERIFY(mask &
HAVE_FACE,
"The ElementTransformation "
642 "for the face has not been configured.");
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)
666 double max_dist = 0.0;
675 os <<
"\nface vertex coordinates (from face transform):\n"
676 <<
"----------------------------------------------\n";
677 coords_base.PrintT(os, coords_base.Height());
686 os <<
"\nface vertex coordinates (from element 1 transform):\n"
687 <<
"---------------------------------------------------\n";
688 coords_el.PrintT(os, coords_el.Height());
692 coords_el -= coords_base;
693 coords_el.Norm2(dist);
694 max_dist = std::max(max_dist, dist.Normlinf());
698 coords_base = coords_el;
707 os <<
"\nface vertex coordinates (from element 2 transform):\n"
708 <<
"---------------------------------------------------\n";
709 coords_el.PrintT(os, coords_el.Height());
711 coords_el -= coords_base;
712 coords_el.Norm2(dist);
713 max_dist = std::max(max_dist, dist.Normlinf());
int GetNPoints() const
Returns the number of the points in the integration rule.
Abstract class for all finite elements.
void Set(const double *p, const int dim)
void Get(double *p, const int dim) const
int GetDim() const
Returns the reference space dimension for the finite element.
void trans(const Vector &u, Vector &x)
Class for an integration rule - an Array of IntegrationPoint.
class LinearPyramidFiniteElement PyramidFE
void SetSize(int s)
Resize the vector to size s.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void CalcAdjugate(const DenseMatrix &a, DenseMatrix &adja)
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Data type dense matrix using column-major storage.
int Size() const
Returns the size of the vector.
int Space() const
Returns the type of FunctionSpace on the element.
double Normlinf() const
Returns the l_infinity norm of the vector.
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Tensor products of polynomials of order k.
PointFiniteElement PointFE
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
void SetType(const int t)
Set the Quadrature1D type of points to use for subdivision.
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
GeometryRefiner GlobGeometryRefiner
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
void GetColumn(int c, Vector &col) const
double DistanceTo(const double *p) const
Compute the Euclidean distance to another vector.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
static bool ProjectPoint(int GeomType, const IntegrationPoint &beg, IntegrationPoint &end)
Project a point end, onto the given Geometry::Type, GeomType.
void Transpose()
(*this) = (*this)^t
int GetDof() const
Returns the number of degrees of freedom in the finite element.
void subtract(const Vector &x, const Vector &y, Vector &z)
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...
void GetColumnReference(int c, Vector &col)
class Linear3DFiniteElement TetrahedronFE
Linear2DFiniteElement TriangleFE
Class for integration point with weight.
int GetType() const
Get the Quadrature1D type of points used for subdivision.
class LinearWedgeFiniteElement WedgeFE
BiLinear2DFiniteElement QuadrilateralFE
void Mult(const double *x, double *y) const
Matrix vector multiplication.
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
TriLinear3DFiniteElement HexahedronFE
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
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