Setup jmatcoeff generation

LLNL · Apr 10, 2024 · dd6b463 · dd6b463
1 parent fa2cb20
commit dd6b463
Show file tree

Hide file tree

Showing 3 changed files with 43 additions and 8 deletions.
diff --git 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
@@ -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
@@ -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();