20 mfem_error (
"NonlinearFormIntegrator::GetLocalStateEnergyPA(...)\n"
21 " is not implemented for this class.");
27 mfem_error (
"NonlinearFormIntegrator::AssemblePA(...)\n"
28 " is not implemented for this class.");
34 mfem_error (
"NonlinearFormIntegrator::AssemblePA(...)\n"
35 " is not implemented for this class.");
41 mfem_error (
"NonlinearFormIntegrator::AssembleGradPA(...)\n"
42 " is not implemented for this class.");
47 mfem_error (
"NonlinearFormIntegrator::AddMultPA(...)\n"
48 " is not implemented for this class.");
53 mfem_error (
"NonlinearFormIntegrator::AddMultGradPA(...)\n"
54 " is not implemented for this class.");
59 mfem_error (
"NonlinearFormIntegrator::AssembleGradDiagonalPA(...)\n"
60 " is not implemented for this class.");
65 mfem_error (
"NonlinearFormIntegrator::AssembleMF(...)\n"
66 " is not implemented for this class.");
71 mfem_error (
"NonlinearFormIntegrator::AddMultMF(...)\n"
72 " is not implemented for this class.");
79 mfem_error(
"NonlinearFormIntegrator::AssembleElementVector"
80 " is not overloaded!");
87 mfem_error(
"NonlinearFormIntegrator::AssembleFaceVector"
88 " is not overloaded!");
95 mfem_error(
"NonlinearFormIntegrator::AssembleElementGrad"
96 " is not overloaded!");
104 mfem_error(
"NonlinearFormIntegrator::AssembleFaceGrad"
105 " is not overloaded!");
111 mfem_error(
"NonlinearFormIntegrator::GetElementEnergy"
112 " is not overloaded!");
123 mfem_error(
"BlockNonlinearFormIntegrator::AssembleElementVector"
124 " is not overloaded!");
134 mfem_error(
"BlockNonlinearFormIntegrator::AssembleFaceVector"
135 " is not overloaded!");
144 mfem_error(
"BlockNonlinearFormIntegrator::AssembleElementGrad"
145 " is not overloaded!");
155 mfem_error(
"BlockNonlinearFormIntegrator::AssembleFaceGrad"
156 " is not overloaded!");
164 mfem_error(
"BlockNonlinearFormIntegrator::GetElementEnergy"
165 " is not overloaded!");
174 return 0.5*(
Z*
Z)/J.
Det();
187 for (
int i = 0; i <
dim; i++)
220 for (
int i = 0; i < dof; i++)
221 for (
int j = 0; j <= i; j++)
224 for (
int d = 0; d <
dim; d++)
229 for (
int k = 0; k <
dim; k++)
230 for (
int l = 0; l <= k; l++)
233 A(i+k*dof,j+l*dof) +=
b;
236 A(j+k*dof,i+l*dof) +=
b;
240 A(i+l*dof,j+k*dof) +=
b;
243 A(j+l*dof,i+k*dof) +=
b;
250 for (
int i = 1; i < dof; i++)
251 for (
int j = 0; j < i; j++)
253 for (
int k = 1; k <
dim; k++)
254 for (
int l = 0; l < k; l++)
257 weight*(
C(i,l)*
G(j,k) -
C(i,k)*
G(j,l) +
258 C(j,k)*
G(i,l) -
C(j,l)*
G(i,k) +
259 t*(
G(i,k)*
G(j,l) -
G(i,l)*
G(j,k)));
261 A(i+k*dof,j+l*dof) +=
a;
262 A(j+l*dof,i+k*dof) +=
a;
264 A(i+l*dof,j+k*dof) -=
a;
265 A(j+k*dof,i+l*dof) -=
a;
294 return 0.5*(
mu*(bI1 -
dim) +
K*(sJ - 1.0)*(sJ - 1.0));
350 for (
int i = 0; i < dof; i++)
351 for (
int k = 0; k <= i; k++)
354 for (
int d = 0; d <
dim; d++)
356 s += DS(i,d)*DS(k,d);
360 for (
int d = 0; d <
dim; d++)
362 A(i+d*dof,k+d*dof) += s;
366 for (
int d = 0; d <
dim; d++)
368 A(k+d*dof,i+d*dof) += s;
375 for (
int i = 0; i < dof; i++)
376 for (
int j = 0; j <
dim; j++)
377 for (
int k = 0; k < dof; k++)
378 for (
int l = 0; l <
dim; l++)
380 A(i+j*dof,k+l*dof) +=
381 a*(
C(i,j)*
G(k,l) +
G(i,j)*
C(k,l)) +
382 b*
G(i,l)*
G(k,j) + c*
G(i,j)*
G(k,l);
459 model->
EvalP(Jpt, P);
509 mfem_error(
"IncompressibleNeoHookeanIntegrator::GetElementEnergy"
510 " has incorrect block finite element space size!");
513 int dof_u = el[0]->GetDof();
514 int dim = el[0]->GetDim();
522 int intorder = 2*el[0]->GetOrder() + 3;
534 el[0]->CalcDShape(ip, DSh_u);
554 mfem_error(
"IncompressibleNeoHookeanIntegrator::AssembleElementVector"
555 " has finite element space of incorrect block number");
558 int dof_u = el[0]->GetDof();
559 int dof_p = el[1]->GetDof();
561 int dim = el[0]->GetDim();
566 mfem_error(
"IncompressibleNeoHookeanIntegrator::AssembleElementVector"
567 " is not defined on manifold meshes");
584 int intorder = 2*el[0]->GetOrder() + 3;
596 el[0]->CalcDShape(ip, DSh_u);
597 Mult(DSh_u, J0i, DS_u);
600 el[1]->CalcShape(ip, Sh_p);
602 real_t pres = Sh_p * *elfun[1];
610 P.
Add(-1.0 * pres * dJ, FinvT);
615 elvec[1]->Add(ip.
weight * Tr.
Weight() * (dJ - 1.0), Sh_p);
626 int dof_u = el[0]->GetDof();
627 int dof_p = el[1]->GetDof();
629 int dim = el[0]->GetDim();
634 elmats(1,1)->
SetSize(dof_p, dof_p);
651 int intorder = 2*el[0]->GetOrder() + 3;
660 el[0]->CalcDShape(ip, DSh_u);
661 Mult(DSh_u, J0i, DS_u);
664 el[1]->CalcShape(ip, Sh_p);
665 real_t pres = Sh_p * *elfun[1];
673 for (
int i_u = 0; i_u < dof_u; ++i_u)
675 for (
int i_dim = 0; i_dim <
dim; ++i_dim)
677 for (
int j_u = 0; j_u < dof_u; ++j_u)
679 for (
int j_dim = 0; j_dim <
dim; ++j_dim)
685 for (
int n=0; n<
dim; ++n)
687 for (
int l=0; l<
dim; ++l)
689 (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
690 dJ * (
mu * F(i_dim, l) - pres * FinvT(i_dim,l)) *
691 FinvT(j_dim,n) * DS_u(i_u,l) * DS_u(j_u, n) *
694 if (j_dim == i_dim && n==l)
696 (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
697 dJ *
mu * DS_u(i_u, l) * DS_u(j_u,n) *
703 (*elmats(0,0))(i_u + i_dim*dof_u, j_u + j_dim*dof_u) +=
704 dJ * pres * FinvT(i_dim, n) *
705 FinvT(j_dim,l) * DS_u(i_u,l) * DS_u(j_u,n) *
715 for (
int i_p = 0; i_p < dof_p; ++i_p)
717 for (
int j_u = 0; j_u < dof_u; ++j_u)
719 for (
int dim_u = 0; dim_u <
dim; ++dim_u)
721 for (
int l=0; l<
dim; ++l)
723 dJ_FinvT_DS = dJ * FinvT(dim_u,l) * DS_u(j_u, l) * Sh_p(i_p) *
725 (*elmats(1,0))(i_p, j_u + dof_u * dim_u) += dJ_FinvT_DS;
726 (*elmats(0,1))(j_u + dof_u * dim_u, i_p) -= dJ_FinvT_DS;
750 const int nd = el.
GetDof();
771 if (Q) { w *= Q->
Eval(T, ip); }
775 gradEF.
Mult(vec1, vec2);
787 const int nd = el.
GetDof();
808 trans.SetIntPoint(&ip);
813 Mult(dshape,
trans.InverseJacobian(), dshapex);
825 trans.AdjugateJacobian().Mult(vec1, vec2);
828 dshape.
Mult(vec2, vec3);
829 MultVWt(shape, vec3, elmat_comp);
831 for (
int ii = 0; ii <
dim; ii++)
833 elmat.
AddMatrix(elmat_comp, ii * nd, ii * nd);
842 for (
int ii = 0; ii <
dim; ii++)
844 for (
int jj = 0; jj <
dim; jj++)
846 elmat.
AddMatrix(w * gradEF(ii, jj), elmat_comp, ii * nd, jj * nd);
859 const int nd = el.
GetDof();
879 trans.SetIntPoint(&ip);
888 trans.AdjugateJacobian().Mult(vec1, vec2);
891 dshape.
Mult(vec2, vec3);
892 MultVWt(shape, vec3, elmat_comp);
894 for (
int ii = 0; ii <
dim; ii++)
896 elmat.
AddMatrix(elmat_comp, ii * nd, ii * nd);
908 const int nd = el.
GetDof();
931 trans.SetIntPoint(&ip);
936 Mult(dshape,
trans.InverseJacobian(), dshapex);
942 trans.AdjugateJacobian().Mult(vec1, vec2);
945 dshape.
Mult(vec2, vec3);
946 MultVWt(shape, vec3, elmat_comp);
949 for (
int ii = 0; ii <
dim; ii++)
951 elmat.
AddMatrix(.5, elmat_comp, ii * nd, ii * nd);
952 elmat.
AddMatrix(-.5, elmat_comp_T, ii * nd, ii * nd);
Dynamic 2D array using row-major layout.
void SetSize(int m, int n)
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
int Size() const
Return the logical size of the array.
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
Data type dense matrix using column-major storage.
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
void Transpose()
(*this) = (*this)^t
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
real_t Trace() const
Trace of a square matrix.
void AddMatrix(DenseMatrix &A, int ro, int co)
Perform (ro+i,co+j)+=A(i,j) for 0<=i.
void UseExternalData(real_t *d, int h, int w)
Change the data array and the size of the DenseMatrix.
void Add(const real_t c, const DenseMatrix &A)
Adds the matrix A multiplied by the number c to the matrix.
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Abstract class for all finite elements.
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
int GetDim() const
Returns the reference space dimension for the finite element.
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...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const =0
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
virtual real_t EvalW(const DenseMatrix &Jpt) const =0
Evaluate the strain energy density function, W = W(Jpt).
ElementTransformation * Ttr
virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const =0
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void SetTransformation(ElementTransformation &Ttr_)
void AssembleElementVector(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
real_t GetElementEnergy(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun) override
Computes the integral of W(Jacobian(Trt)) over a target zone.
const IntegrationRule * GetDefaultIntegrationRule(const FiniteElement &trial_fe, const FiniteElement &test_fe, const ElementTransformation &trans) const override
Subclasses should override to choose a default integration rule.
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &Ttr, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
real_t GetElementEnergy(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun) override
Compute the local energy.
void AssembleElementVector(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array< Vector * > &elvec) override
Perform the local action of the NonlinearFormIntegrator.
void AssembleElementGrad(const Array< const FiniteElement * > &el, ElementTransformation &Tr, const Array< const Vector * > &elfun, const Array2D< DenseMatrix * > &elmats) override
Assemble the local gradient matrix.
Class for integration point with weight.
Class for an integration rule - an Array of IntegrationPoint.
int GetNPoints() const
Returns the number of the points in the integration rule.
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
const IntegrationRule * GetIntegrationRule() const
Equivalent to GetIntRule, but retained for backward compatibility with applications.
void EvalP(const DenseMatrix &J, DenseMatrix &P) const override
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const override
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
real_t EvalW(const DenseMatrix &J) const override
Evaluate the strain energy density function, W = W(Jpt).
void AssembleH(const DenseMatrix &J, const DenseMatrix &DS, const real_t weight, DenseMatrix &A) const override
Evaluate the derivative of the 1st Piola-Kirchhoff stress tensor and assemble its contribution to the...
real_t EvalW(const DenseMatrix &J) const override
Evaluate the strain energy density function, W = W(Jpt).
void EvalP(const DenseMatrix &J, DenseMatrix &P) const override
Evaluate the 1st Piola-Kirchhoff stress tensor, P = P(Jpt).
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
void AssembleElementVector(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, Vector &elvect) override
Perform the local action of the NonlinearFormIntegrator.
void AssembleElementGrad(const FiniteElement &el, ElementTransformation &trans, const Vector &elfun, DenseMatrix &elmat) override
Assemble the local gradient matrix.
static const IntegrationRule & GetRule(const FiniteElement &fe, const ElementTransformation &T)
void SetSize(int s)
Resize the vector to size s.
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
void trans(const Vector &u, Vector &x)
void mfem_error(const char *msg)
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
void CalcAdjugateTranspose(const DenseMatrix &a, DenseMatrix &adjat)
Calculate the transposed adjugate of a matrix (for NxN matrices, N=1,2,3)
void MultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
void AddMultVWt(const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += v w^t.
void MultVVt(const Vector &v, DenseMatrix &vvt)
Make a matrix from a vector V.Vt.
void AddMultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
ABt += A * B^t.
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
void CalcInverseTranspose(const DenseMatrix &a, DenseMatrix &inva)
Calculate the inverse transpose of a matrix (for NxN matrices, N=1,2,3)
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)