diff --git a/include/nlelast_integ.hpp b/include/nlelast_integ.hpp index f4f4c383..9b2cbef4 100644 --- a/include/nlelast_integ.hpp +++ b/include/nlelast_integ.hpp @@ -24,6 +24,9 @@ namespace mfem void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; } double EvalW(const DenseMatrix &J); + + double EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); + void EvalP(const FiniteElement &el, const IntegrationPoint &ip, const DenseMatrix &PMatI, FaceElementTransformations &Trans, DenseMatrix &P); void EvalP(const FiniteElement &el, const IntegrationPoint &ip, const DenseMatrix &PMatI, ElementTransformation &Trans, DenseMatrix &P); diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp index 4f820675..1766345a 100644 --- a/src/nlelast_integ.cpp +++ b/src/nlelast_integ.cpp @@ -11,6 +11,12 @@ namespace mfem double TestLinModel::EvalW(const DenseMatrix &J) {MFEM_ABORT("TODO")}; +double TestLinModel::EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) +{ +const double wL = w * c_lambda->Eval(Ttr, ip); +const double wM = w * c_mu->Eval(Ttr, ip); +return wL + 2.0*wM; +} void TestLinModel::EvalP( const FiniteElement &el, const IntegrationPoint &ip, const DenseMatrix &PMatI, FaceElementTransformations &Trans, DenseMatrix &P) @@ -257,6 +263,14 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, Jrt.SetSize(dim); elvect.SetSize(nvdofs); elvect = 0.0; + + const bool kappa_is_nonzero = (kappa != 0.0); + if (kappa_is_nonzero) + { + jmat.SetSize(nvdofs); + jmat = 0.; + } + model->SetTransformation(Trans); shape1.SetSize(ndofs1); @@ -318,9 +332,10 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, CalcInverse(Trans.Jacobian(), Jrt); double w = ip.weight; + double wLM = 0.0; if (ndofs2) { - w /= 2.0; + w /= 2.0; el2.CalcShape(eip2, shape2); el2.CalcDShape(eip2, DSh2); @@ -330,6 +345,7 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, model->EvalP(el2, eip2, PMatI2, Trans, P2); double w2 = w / Trans.Elem2->Weight(); + wLM = model->EvalwLM(w2, *Trans.Elem2, eip2); P2 *= w2; P2.Mult(nor, tau2); } @@ -342,9 +358,12 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, model->EvalP(el1, eip1, PMatI1, Trans, P1); double w1 = w / Trans.Elem1->Weight(); + wLM += model->EvalwLM(w1, *Trans.Elem1, eip1); P1 *= w1; P1.Mult(nor, tau1); + const double jmatcoef = kappa * (nor*nor) * wLM; + // (1,1) block for (int im = 0, i = 0; im < dim; ++im) { @@ -407,6 +426,13 @@ const int dim = el1.GetDim(); elmat = 0.0; model->SetTransformation(Tr); + const bool kappa_is_nonzero = (kappa != 0.0); + if (kappa_is_nonzero) + { + jmat.SetSize(nvdofs); + jmat = 0.; + } + shape1.SetSize(ndofs1); elfun1.MakeRef(elfun_copy,0,ndofs1*dim); PMatI1.UseExternalData(elfun1.GetData(), ndofs1, dim); @@ -478,6 +504,7 @@ const int dim = el1.GetDim(); CalcInverse(Tr.Jacobian(), Jrt); double w = ip.weight; + double wLM = 0.0; if (ndofs2) { w /= 2.0; @@ -491,6 +518,7 @@ const int dim = el1.GetDim(); //model->EvalDmat(dim, ndofs1, eip2, Tr, DS2, Dmat2); model->EvalDmat(dim, ndofs2, eip2, Tr, dshape2_ps, Dmat2); double w2 = w / Tr.Elem2->Weight(); + //wLM = (wL2 + 2.0*wM2); wnor2.Set(w2,nor); } @@ -505,8 +533,9 @@ const int dim = el1.GetDim(); CalcAdjugate(Tr.Elem1->Jacobian(), adjJ1); Mult(DSh1, adjJ1, dshape1_ps); - model->EvalDmat(dim, ndofs1, eip1, Tr, dshape1_ps, Dmat1); + + const double jmatcoef = kappa * (nor*nor) * wLM; // (1,1) block //works wnor1.Set(w1,nor); @@ -526,7 +555,11 @@ const int dim = el1.GetDim(); } elmat *= -1.0; - }; + if (kappa_is_nonzero) + { + elmat += jmat; + } + }; void DGHyperelasticNLFIntegrator::AssembleBlock( const int dim, const int row_ndofs, const int col_ndofs, diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp index 429999b3..f1d4ea75 100644 --- a/test/nlelast_test.cpp +++ b/test/nlelast_test.cpp @@ -354,7 +354,6 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) //x[i] = i; } - y1.SetSize(fespace.GetTrueVSize()); y1 = 0.0; @@ -390,9 +389,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) y1/=y1.Norml2(); y2/=y2.Norml2(); - - - cout << "Scaled Linear residual norm: " << y1.Norml2() << endl; + cout << "Scaled Linear residual norm: " << y1.Norml2() << endl; cout << "Scaled Nonlinear residual norm: " << y2.Norml2() << endl; y1 -= y2; @@ -416,7 +413,6 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) cout << "Linear Stiffness matrix norm: " << a1.SpMat().MaxNorm() << endl; cout << "Stiffness matrix difference norm: " << norm_diff << endl; - LinearForm b1(&fespace); b1.AddBdrFaceIntegrator( new DGElasticityDirichletLFIntegrator( @@ -429,6 +425,9 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) init_x, lambda_c, mu_c, alpha, kappa), dir_bdr); b2.Assemble(); + cout << "Linear RHS norm: " << b1.Norml2() << endl; + cout << "Nonlinear RHS norm: " << b2.Norml2() << endl; + b1 -= b2; norm_diff = b1.Norml2();