Interiorfaceintegrator works

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

Hide file tree

Showing 2 changed files with 36 additions and 30 deletions.
diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp
@@ -450,6 +450,11 @@ const int dim = el1.GetDim();
    Dmat1 = 0.0;
    Dmat2 = 0.0;
 
+   DenseMatrix adjJ2(dim);
+
+   DenseMatrix dshape1_ps(ndofs1, dim);
+   DenseMatrix dshape2_ps(ndofs2, dim);
+
    for (int i = 0; i < ir->GetNPoints(); i++)
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
@@ -476,10 +481,15 @@ const int dim = el1.GetDim();
       if (ndofs2)
       {
       w /= 2.0;
+      el2.CalcShape(eip2, shape2);
       el2.CalcDShape(eip2, DSh2);
       Mult(DSh2, Jrt, DS2);
 
-      model->EvalDmat(dim, ndofs1, eip2, Tr, DS2, Dmat2);
+      CalcAdjugate(Tr.Elem2->Jacobian(), adjJ2);
+      Mult(DSh2, adjJ2, dshape2_ps);
+
+      //model->EvalDmat(dim, ndofs1, eip2, Tr, DS2, Dmat2);
+      model->EvalDmat(dim, ndofs2, eip2, Tr, dshape2_ps, Dmat2);
       double w2 = w / Tr.Elem2->Weight();
       wnor2.Set(w2,nor);
       }
@@ -491,32 +501,28 @@ const int dim = el1.GetDim();
       double w1 = w / Tr.Elem1->Weight();
 
       // Temporary stuff
-      Vector nL1(dim);
-      Vector nM1(dim);
-      DenseMatrix adjJ(dim);
-      DenseMatrix dshape1_ps(ndofs1,dim);
-      Vector dshape1_dnM(ndofs1);
-      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ);
-
-      Mult(DSh1, adjJ, dshape1_ps);
+      DenseMatrix adjJ1(dim);
+      CalcAdjugate(Tr.Elem1->Jacobian(), adjJ1);
+      Mult(DSh1, adjJ1, dshape1_ps);
+
 
       model->EvalDmat(dim, ndofs1, eip1, Tr, dshape1_ps, Dmat1);
 
-      // (1,1) block
+      // (1,1) block //works
       wnor1.Set(w1,nor);
       AssembleBlock(dim, ndofs1, ndofs1, 0, 0, shape1, wnor1, Dmat1, elmat);
 
       if (ndofs2 == 0) {continue;}
-       shape2.Neg();
+       shape2.Neg(); 
 
-       // (1,2) block
-      AssembleBlock(dim, ndofs1, ndofs2, 0, 0, shape1, wnor2, Dmat2, elmat);
+       // (1,2) block works
+      AssembleBlock(dim, ndofs1, ndofs2, 0, dim*ndofs1, shape1, wnor2, Dmat2, elmat);
 
        // (2,1) block
-      AssembleBlock(dim, ndofs2, ndofs1, 0, 0, shape2, wnor1, Dmat1, elmat);
+      AssembleBlock(dim, ndofs2, ndofs1, dim*ndofs1, 0, shape2, wnor1, Dmat1, elmat); 
 
        // (2,2) block
-      AssembleBlock(dim, ndofs2, ndofs2, 0, 0, shape2, wnor2, Dmat2, elmat);
+      AssembleBlock(dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, shape2, wnor2, Dmat2, elmat);
 
    }
    elmat *= -1.0;

diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp
@@ -223,26 +223,26 @@ namespace mfem
 
        const double jmatcoef = 0.0;
 
-       // (1,1) block
-       _AssembleBlock(
+       // (1,1) block works
+        _AssembleBlock(
           dim, ndofs1, ndofs1, 0, 0, jmatcoef, nL1, nM1,
-          shape1, shape1, dshape1_dnM, dshape1_ps, elmat, jmat);
+          shape1, shape1, dshape1_dnM, dshape1_ps, elmat, jmat); 
 
        if (ndofs2 == 0) { continue; }
 
        // In both elmat and jmat, shape2 appears only with a minus sign.
        shape2.Neg();
 
-       // (1,2) block
+       // (1,2) block works
        _AssembleBlock(
           dim, ndofs1, ndofs2, 0, dim*ndofs1, jmatcoef, nL2, nM2,
-          shape1, shape2, dshape2_dnM, dshape2_ps, elmat, jmat);
+          shape1, shape2, dshape2_dnM, dshape2_ps, elmat, jmat); 
        // (2,1) block
        _AssembleBlock(
           dim, ndofs2, ndofs1, dim*ndofs1, 0, jmatcoef, nL1, nM1,
-          shape2, shape1, dshape1_dnM, dshape1_ps, elmat, jmat);
+          shape2, shape1, dshape1_dnM, dshape1_ps, elmat, jmat); 
        // (2,2) block
-       _AssembleBlock(
+        _AssembleBlock(
           dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, jmatcoef, nL2, nM2,
           shape2, shape2, dshape2_dnM, dshape2_ps, elmat, jmat); 
     }
@@ -300,7 +300,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
    lambda = _lambda;      // Set lambda for all element attributes.
    PWConstCoefficient lambda_c(lambda);
    Vector mu(mesh.attributes.Max());
-   double _mu = 1.0;      
+   double _mu = 2.33;      
    mu = _mu;       // Set mu = 1 for all element attributes.
    PWConstCoefficient mu_c(mu);
 
@@ -310,23 +310,23 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
    dir_bdr[1] = 1; // boundary attribute 2 is Dirichlet
 
    BilinearForm a1(&fespace);
-   //a1.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c));
-   //a1.AddInteriorFaceIntegrator(
-      //new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa));
+   a1.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c));
+   a1.AddInteriorFaceIntegrator(
+      new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa));
    a1.AddBdrFaceIntegrator(
       new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr);
    a1.Assemble();
    cout<<"a1.Height() is: "<<a1.Height()<<endl;
 
    TestLinModel model(_mu, _lambda);
    NonlinearForm a2(&fespace);
-   //a2.AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model));
-   //a2.AddInteriorFaceIntegrator(
-      //new DGHyperelasticNLFIntegrator(&model, alpha, kappa));
+   a2.AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model));
+   a2.AddInteriorFaceIntegrator(
+      new DGHyperelasticNLFIntegrator(&model, alpha, kappa));
    a2.AddBdrFaceIntegrator(
       new DGHyperelasticNLFIntegrator(&model, alpha, kappa), dir_bdr);
 
-      /* BilinearForm a2(&fespace);
+/*       BilinearForm a2(&fespace);
    //a1.AddDomainIntegrator(new ElasticityIntegrator(lambda_c, mu_c));
    a2.AddInteriorFaceIntegrator(
       new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa)); //Needed??