diff --git a/include/nlelast_integ.hpp b/include/nlelast_integ.hpp index 8fe660db..4aa681d5 100644 --- a/include/nlelast_integ.hpp +++ b/include/nlelast_integ.hpp @@ -11,7 +11,7 @@ namespace mfem { - class TestLinModel //: public HyperelasticModel + class LinElastMaterialModel //: public HyperelasticModel { protected: mutable double mu, lambda; @@ -19,20 +19,19 @@ namespace mfem ElementTransformation *Ttr; public: - TestLinModel(double mu_, double lambda_) + LinElastMaterialModel(double mu_, double lambda_) : mu(mu_), lambda(lambda_) { c_mu = new ConstantCoefficient(mu), c_lambda = new ConstantCoefficient(lambda); } void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; } - void SetML(ElementTransformation &Trans, const IntegrationPoint &ip); - void SetML(FaceElementTransformations &Trans, const IntegrationPoint &ip); + void SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip); + void SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip); double EvalW(const DenseMatrix &J); - double EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); + double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); void EvalP(const DenseMatrix &J, DenseMatrix &P); void EvalDmat(const int dim, const int dof, const DenseMatrix gshape, DenseMatrix &Dmat); - // void AssembleH(const int dim, const int dof, const DenseMatrix &J, const double w, DenseMatrix &elmat); }; class StVenantKirchhoffModel //: public HyperelasticModel @@ -50,7 +49,7 @@ namespace mfem void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; } double EvalW(const DenseMatrix &J); - double EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); + double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); void EvalP(const DenseMatrix &Jpt, const double lambda, const double mu, DenseMatrix &P); @@ -74,7 +73,7 @@ namespace mfem void SetTransformation(ElementTransformation &Ttr_) { Ttr = &Ttr_; } double EvalW(const DenseMatrix &J); - double EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); + double EvalDGWeight(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip); void EvalP(const DenseMatrix &Jpt, const double lambda, const double mu, DenseMatrix &P); @@ -89,7 +88,7 @@ namespace mfem { public: - DGHyperelasticNLFIntegrator(TestLinModel *m, double alpha_, double kappa_) + DGHyperelasticNLFIntegrator(LinElastMaterialModel *m, double alpha_, double kappa_) : HyperReductionIntegrator(false), model(m), lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) {} virtual void AssembleFaceVector(const FiniteElement &el1, @@ -103,7 +102,7 @@ namespace mfem const Vector &elfun, DenseMatrix &elmat); // values of all scalar basis functions for one component of u (which is a protected: - TestLinModel *model; + LinElastMaterialModel *model; Coefficient *lambda, *mu; double alpha, kappa; // vector) at the integration point in the reference space #ifndef MFEM_THREAD_SAFE @@ -148,7 +147,7 @@ namespace mfem { private: - TestLinModel *model; + LinElastMaterialModel *model; // Jrt: the Jacobian of the target-to-reference-element transformation. // Jpr: the Jacobian of the reference-to-physical-element transformation. // Jpt: the Jacobian of the target-to-physical-element transformation. @@ -162,7 +161,7 @@ namespace mfem DenseMatrix DSh, DS, Jrt, Jpr, Jpt, P, PMatI, PMatO, Dmat; public: - HyperelasticNLFIntegratorHR(TestLinModel *m) : HyperReductionIntegrator(false), model(m) + HyperelasticNLFIntegratorHR(LinElastMaterialModel *m) : HyperReductionIntegrator(false), model(m) { } @@ -186,7 +185,7 @@ namespace mfem { protected: VectorCoefficient &uD; - TestLinModel *model; + LinElastMaterialModel *model; Coefficient *lambda, *mu; double alpha, kappa; @@ -203,7 +202,7 @@ namespace mfem public: DGHyperelasticDirichletNLFIntegrator(VectorCoefficient &uD_, - TestLinModel *m, + LinElastMaterialModel *m, double alpha_, double kappa_) : uD(uD_), model(m), lambda(NULL), mu(NULL), alpha(alpha_), kappa(kappa_) {} diff --git a/src/nlelast_integ.cpp b/src/nlelast_integ.cpp index ec51ae8b..c837214d 100644 --- a/src/nlelast_integ.cpp +++ b/src/nlelast_integ.cpp @@ -8,7 +8,7 @@ using namespace std; namespace mfem { /* -double StVenantKirchhoffModel::EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) +double StVenantKirchhoffModel::EvalDGWeight(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; @@ -135,29 +135,29 @@ void StVenantKirchhoffModel::AssembleH(const DenseMatrix &Dmat, const DenseMatri } }; */ -double TestLinModel::EvalW(const DenseMatrix &J) +double LinElastMaterialModel::EvalW(const DenseMatrix &J) {MFEM_ABORT("TODO")}; -double TestLinModel::EvalwLM(const double w, ElementTransformation &Ttr, const IntegrationPoint &ip) +double LinElastMaterialModel::EvalDGWeight(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::SetML(ElementTransformation &Trans, const IntegrationPoint &ip) +void LinElastMaterialModel::SetMatParam(ElementTransformation &Trans, const IntegrationPoint &ip) { mu = c_mu->Eval(Trans, ip); lambda = c_lambda->Eval(Trans, ip); } -void TestLinModel::SetML(FaceElementTransformations &Trans, const IntegrationPoint &ip) +void LinElastMaterialModel::SetMatParam(FaceElementTransformations &Trans, const IntegrationPoint &ip) { mu = c_mu->Eval(Trans, ip); lambda = c_lambda->Eval(Trans, ip); } -void TestLinModel::EvalP(const DenseMatrix &J, DenseMatrix &P) +void LinElastMaterialModel::EvalP(const DenseMatrix &J, DenseMatrix &P) { // stress = 2*M*e(u) + L*tr(e(u))*I, where // e(u) = (1/2)*(grad(u) + grad(u)^T) @@ -187,7 +187,7 @@ void TestLinModel::EvalP(const DenseMatrix &J, DenseMatrix &P) } } -void TestLinModel::EvalDmat(const int dim, const int dof, const DenseMatrix gshape, DenseMatrix &Dmat) +void LinElastMaterialModel::EvalDmat(const int dim, const int dof, const DenseMatrix gshape, DenseMatrix &Dmat) { for (size_t i = 0; i < dim; i++) { @@ -353,13 +353,13 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, CalcAdjugate(Trans.Elem2->Jacobian(), adjJ2); Mult(DSh2, adjJ2, DS2); MultAtB(PMatI2, DS2, Jpt2); - model->SetML(Trans, eip2); + model->SetMatParam(Trans, eip2); model->EvalP(Jpt2, P2); double w2 = w / Trans.Elem2->Weight(); - wLM = model->EvalwLM(w2, *Trans.Elem2, eip2); + wLM = model->EvalDGWeight(w2, *Trans.Elem2, eip2); P2 *= w2; P2.Mult(nor, tau2); } @@ -370,13 +370,13 @@ void DGHyperelasticNLFIntegrator::AssembleFaceVector(const FiniteElement &el1, CalcAdjugate(Trans.Elem1->Jacobian(), adjJ1); Mult(DSh1, adjJ1, DS1); MultAtB(PMatI1, DS1, Jpt1); - model->SetML(Trans, eip1); + model->SetMatParam(Trans, eip1); model->EvalP(Jpt1, P1); double w1 = w / Trans.Elem1->Weight(); - wLM += model->EvalwLM(w1, *Trans.Elem1, eip1); + wLM += model->EvalDGWeight(w1, *Trans.Elem1, eip1); P1 *= w1; P1.Mult(nor, tau1); @@ -573,12 +573,12 @@ const int dim = el1.GetDim(); Mult(DSh2, adjJ2, dshape2_ps); //model->EvalDmat(dim, ndofs1, eip2, Tr, DS2, Dmat2); - model->SetML(Tr,eip2); + model->SetMatParam(Tr,eip2); model->EvalDmat(dim, ndofs2, dshape2_ps, Dmat2); //model->EvalDmat(dim, ndofs2, eip2, Tr, dshape2_ps, Dmat2); double w2 = w / Tr.Elem2->Weight(); - wLM = model->EvalwLM(w2, *Tr.Elem2, eip2); + wLM = model->EvalDGWeight(w2, *Tr.Elem2, eip2); wnor2.Set(w2,nor); } @@ -587,7 +587,7 @@ const int dim = el1.GetDim(); Mult(DSh1, Jrt, DS1); double w1 = w / Tr.Elem1->Weight(); - wLM += model->EvalwLM(w1, *Tr.Elem1, eip1); + wLM += model->EvalDGWeight(w1, *Tr.Elem1, eip1); // Temporary stuff DenseMatrix adjJ1(dim); @@ -595,7 +595,7 @@ const int dim = el1.GetDim(); Mult(DSh1, adjJ1, dshape1_ps); //model->EvalDmat(dim, ndofs1, eip1, Tr, dshape1_ps, Dmat1); - model->SetML(Tr,eip1); + model->SetMatParam(Tr,eip1); model->EvalDmat(dim, ndofs1, dshape1_ps, Dmat1); @@ -753,7 +753,7 @@ void HyperelasticNLFIntegratorHR::AssembleElementGrad(const FiniteElement &el, Mult(DSh, Jrt, DS); MultAtB(PMatI, DS, Jpt); - model->SetML(Ttr,ip); + model->SetMatParam(Ttr,ip); model->EvalDmat(dim, dof, DS, Dmat); AssembleH(dim, dof, ip.weight * Ttr.Weight(), DS, elmat); @@ -859,7 +859,7 @@ MFEM_ASSERT(Tr.Elem2No < 0, "interior boundary is not supported"); double wLM; const double w = ip.weight / Tr.Elem1->Weight(); - wLM = model->EvalwLM(w, *Tr.Elem1, eip); + wLM = model->EvalDGWeight(w, *Tr.Elem1, eip); jcoef = kappa * wLM * (nor*nor); diff --git a/test/nlelast_test.cpp b/test/nlelast_test.cpp index 6495610d..b829c9a7 100644 --- a/test/nlelast_test.cpp +++ b/test/nlelast_test.cpp @@ -346,7 +346,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast) a1.Assemble(); cout<<"a1.Height() is: "<