From 6af883f260a646dabc8c296c7b786d4a88d65210 Mon Sep 17 00:00:00 2001 From: Axel Larsson Date: Wed, 10 Apr 2024 13:22:06 -0700 Subject: [PATCH] Kappa term in integral --- include/nlelast_integ.hpp | 5 +++ src/nlelast_integ.cpp | 66 ++++++++++++++++++++++++++++++++++++--- test/nlelast_test.cpp | 19 +++++------ 3 files changed, 77 insertions(+), 13 deletions(-) diff --git a/include/nlelast_integ.hpp b/include/nlelast_integ.hpp index bf126f5c..5d5b3b2d 100644 --- a/include/nlelast_integ.hpp +++ b/include/nlelast_integ.hpp @@ -91,6 +91,11 @@ class DGHyperelasticNLFIntegrator : virtual public HyperReductionIntegrator const Vector &col_shape, const double jmatcoef, const Vector &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat, DenseMatrix &jmat); + static void AssembleJmat( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef,DenseMatrix &jmat); + }; // Domain integrator for nonlinear elastic DG. diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp index cb01efc0..0e1434c6 100644 --- a/src/nlelast_integ.cpp +++ b/src/nlelast_integ.cpp @@ -244,6 +244,25 @@ void _PrintVector(const Vector &vec, return; } +void DGHyperelasticNLFIntegrator::AssembleJmat( + const int dim, const int row_ndofs, const int col_ndofs, + const int row_offset, const int col_offset, const Vector &row_shape, + const Vector &col_shape, const double jmatcoef,DenseMatrix &jmat){ + for (int d = 0; d < dim; ++d) + { + const int jo = col_offset + d*col_ndofs; + const int io = row_offset + d*row_ndofs; + for (int jdof = 0, j = jo; jdof < col_ndofs; ++jdof, ++j) + { + const double sj = jmatcoef * col_shape(jdof); + for (int i = max(io,j), idof = i - io; idof < row_ndofs; ++idof, ++i) + { + jmat(i, j) += row_shape(idof) * sj; + } + } + } +}; + // Boundary integrator void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, @@ -308,8 +327,8 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, Vector big_row1(dim*ndofs1); Vector big_row2(dim*ndofs2); - //for (int i = 0; i < ir->GetNPoints(); i++) - for (int i = 0; i < 1; i++) + //for (int i = 0; i < 1; i++) + for (int i = 0; i < ir->GetNPoints(); i++) { const IntegrationPoint &ip = ir->IntPoint(i); @@ -371,10 +390,49 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, { elvect(i) += shape1(idof) * tau1(im); } + } + if (ndofs2 != 0) { + shape2.Neg(); } + if (kappa_is_nonzero) + { + jmat = 0.; + AssembleJmat( + dim, ndofs1, ndofs1, 0, 0, shape1, + shape1, jmatcoef, jmat); + if (ndofs2 != 0) { + AssembleJmat( + dim, ndofs1, ndofs2, 0, ndofs1*dim, shape1, + shape2, jmatcoef, jmat); + AssembleJmat( + dim, ndofs2, ndofs1, ndofs1*dim, 0, shape2, + shape1, jmatcoef, jmat); + AssembleJmat( + dim, ndofs2, ndofs2, ndofs1*dim, ndofs1*dim, shape2, + shape2, jmatcoef, jmat); + } + + //Flatten jmat + for (int i = 0; i < nvdofs; ++i) + { + for (int j = 0; j < i; ++j) + { + jmat(j,i) = jmat(i,j); + } + } + // Apply jmat + for (size_t i = 0; i < nvdofs; i++) + { + for (size_t j = 0; j < nvdofs; j++) + { + elvect(i) -= jmat(i,j) * elfun(j); + } + } + } + if (ndofs2 == 0) {continue;} - shape2.Neg(); + // (1,2) block for (int im = 0, i = 0; im < dim; ++im) @@ -555,7 +613,7 @@ const int dim = el1.GetDim(); AssembleBlock(dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, shape2, shape2,jmatcoef,wnor2, Dmat2, elmat,jmat); } - + // elmat := -elmat + jmat elmat *= -1.0; if (kappa_is_nonzero) diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp index 0077c14c..60647395 100644 --- a/test/nlelast_test.cpp +++ b/test/nlelast_test.cpp @@ -169,6 +169,7 @@ namespace mfem ir = &IntRules.Get(Trans.GetGeometryType(), order); } + //for (int pind = 0; pind < 1; ++pind) for (int pind = 0; pind < ir->GetNPoints(); ++pind) { const IntegrationPoint &ip = ir->IntPoint(pind); @@ -272,7 +273,7 @@ namespace mfem elmat(i,i) += jmat(i,i); } } - + PrintMatrix(elmat, "checkmat.txt"); } /* PrintMatrix(elmat, "checkmat.txt"); } */ @@ -316,7 +317,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) int dim = mesh.Dimension(); int order = 1; double alpha = 0.0; // IIPG - double kappa = -1.0; // TODO: Enable kappa = -1.0 + double kappa = -1.0; DG_FECollection fec(order, dim, BasisType::GaussLobatto); FiniteElementSpace fespace(&mesh, &fec, dim); @@ -337,19 +338,19 @@ 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 Test_DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr); + new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr); a1.Assemble(); cout<<"a1.Height() is: "<