From b255718c4810ac7dc9163ac7b2d7605d83286d18 Mon Sep 17 00:00:00 2001 From: Axel Larsson Date: Wed, 10 Apr 2024 11:12:31 -0700 Subject: [PATCH] Stiffness matrix handles kappa --- include/nlelast_integ.hpp | 4 +- src/nlelast_integ.cpp | 53 +++++++++++++----- test/nlelast_test.cpp | 110 +++++++++++++++++++++++++------------- 3 files changed, 115 insertions(+), 52 deletions(-) diff --git a/include/nlelast_integ.hpp b/include/nlelast_integ.hpp index 9b2cbef4..bf126f5c 100644 --- a/include/nlelast_integ.hpp +++ b/include/nlelast_integ.hpp @@ -27,7 +27,6 @@ namespace mfem 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); @@ -89,7 +88,8 @@ class DGHyperelasticNLFIntegrator : virtual public HyperReductionIntegrator static void AssembleBlock( 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 &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat); + const Vector &col_shape, const double jmatcoef, + const Vector &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat, DenseMatrix &jmat); }; diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp index 1766345a..cb01efc0 100644 --- a/src/nlelast_integ.cpp +++ b/src/nlelast_integ.cpp @@ -308,8 +308,8 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, Vector big_row1(dim*ndofs1); Vector big_row2(dim*ndofs2); - //for (int i = 0; i < 1; i++) - for (int i = 0; i < ir->GetNPoints(); i++) + //for (int i = 0; i < ir->GetNPoints(); i++) + for (int i = 0; i < 1; i++) { const IntegrationPoint &ip = ir->IntPoint(i); @@ -518,7 +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); + wLM = model->EvalwLM(w2, *Tr.Elem2, eip2); wnor2.Set(w2,nor); } @@ -527,6 +527,7 @@ const int dim = el1.GetDim(); Mult(DSh1, Jrt, DS1); double w1 = w / Tr.Elem1->Weight(); + wLM += model->EvalwLM(w1, *Tr.Elem1, eip1); // Temporary stuff DenseMatrix adjJ1(dim); @@ -536,35 +537,47 @@ const int dim = el1.GetDim(); model->EvalDmat(dim, ndofs1, eip1, Tr, dshape1_ps, Dmat1); const double jmatcoef = kappa * (nor*nor) * wLM; - + // (1,1) block //works wnor1.Set(w1,nor); - AssembleBlock(dim, ndofs1, ndofs1, 0, 0, shape1, wnor1, Dmat1, elmat); + AssembleBlock(dim, ndofs1, ndofs1, 0, 0, shape1, shape1, jmatcoef, wnor1,Dmat1, elmat,jmat); if (ndofs2 == 0) {continue;} shape2.Neg(); // (1,2) block works - AssembleBlock(dim, ndofs1, ndofs2, 0, dim*ndofs1, shape1, wnor2, Dmat2, elmat); + AssembleBlock(dim, ndofs1, ndofs2, 0, dim*ndofs1, shape1, shape2,jmatcoef,wnor2, Dmat2, elmat,jmat); // (2,1) block - AssembleBlock(dim, ndofs2, ndofs1, dim*ndofs1, 0, shape2, wnor1, Dmat1, elmat); + AssembleBlock(dim, ndofs2, ndofs1, dim*ndofs1, 0, shape2, shape1,jmatcoef,wnor1, Dmat1, elmat,jmat); // (2,2) block - AssembleBlock(dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, shape2, wnor2, Dmat2, elmat); + AssembleBlock(dim, ndofs2, ndofs2, dim*ndofs1, dim*ndofs1, shape2, shape2,jmatcoef,wnor2, Dmat2, elmat,jmat); } - elmat *= -1.0; - if (kappa_is_nonzero) + + // elmat := -elmat + jmat + elmat *= -1.0; + if (kappa_is_nonzero) { - elmat += jmat; + for (int i = 0; i < nvdofs; ++i) + { + for (int j = 0; j < i; ++j) + { + double mij = jmat(i,j); + elmat(i,j) += mij; + elmat(j,i) += mij; + } + elmat(i,i) += jmat(i,i); + } } }; void DGHyperelasticNLFIntegrator::AssembleBlock( 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 &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat){ + const Vector &col_shape, const double jmatcoef, + const Vector &wnor, const DenseMatrix &Dmat, DenseMatrix &elmat, DenseMatrix &jmat){ for (int n = 0, jj = col_offset; n < dim; ++n) { for (int m = 0; m < col_ndofs; ++m, ++jj) @@ -585,6 +598,22 @@ for (int n = 0, jj = col_offset; n < dim; ++n) } } } + + if (jmatcoef == 0.0) { return; } + + 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; + } + } + } }; // Domain integrator diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp index f1d4ea75..0077c14c 100644 --- a/test/nlelast_test.cpp +++ b/test/nlelast_test.cpp @@ -62,7 +62,7 @@ namespace mfem DenseMatrix &elmat, DenseMatrix &jmat); }; - void _AssembleBlock( + void _AssembleBlock( const int dim, const int row_ndofs, const int col_ndofs, const int row_offset, const int col_offset, const double jmatcoef, const Vector &col_nL, const Vector &col_nM, @@ -83,11 +83,26 @@ namespace mfem for (int idof = 0; idof < row_ndofs; ++idof, ++i) { elmat(i, j) += row_shape(idof) * tt; - //elmat(i, j) += tt; } } } } + + if (jmatcoef == 0.0) { return; } + + 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; + } + } + } }; void Test_DGElasticityIntegrator::AssembleFaceMatrix( @@ -109,11 +124,7 @@ namespace mfem const int dim = el1.GetDim(); const int ndofs1 = el1.GetDof(); - //const int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0; - - int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0; // TEMP: Prevents resizing of elmat - - + const int ndofs2 = (Trans.Elem2No >= 0) ? el2.GetDof() : 0; const int nvdofs = dim*(ndofs1 + ndofs2); // Initially 'elmat' corresponds to the term: @@ -123,9 +134,7 @@ namespace mfem // elmat := -elmat + alpha*elmat^T + jmat elmat.SetSize(nvdofs); elmat = 0.; - - //ndofs2 = 0; // TEMP: Prevents resizing of elmat - + const bool kappa_is_nonzero = (kappa != 0.0); if (kappa_is_nonzero) { @@ -160,7 +169,6 @@ 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); @@ -221,34 +229,53 @@ namespace mfem dshape1_ps.Mult(nM1, dshape1_dnM); } - const double jmatcoef = 0.0; + const double jmatcoef = kappa * (nor*nor) * wLM; + cout<<"jmatcoef is: "<