Skip to content

Commit

Permalink
Renamed variablesRenamed material model name and some equations
Browse files Browse the repository at this point in the history
  • Loading branch information
@larsson4
larsson4 committed Apr 17, 2024
1 parent 2b67820 commit f384ea8
Show file tree
Hide file tree
Showing 3 changed files with 31 additions and 32 deletions.
27 changes: 13 additions & 14 deletions include/nlelast_integ.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,28 +11,27 @@
namespace mfem
{

class TestLinModel //: public HyperelasticModel
class LinElastMaterialModel //: public HyperelasticModel
{
protected:
mutable double mu, lambda;
Coefficient *c_mu, *c_lambda;
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
Expand All @@ -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);

Expand All @@ -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);

Expand All @@ -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,
Expand All @@ -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
Expand Down Expand Up @@ -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.
Expand All @@ -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)
{
}

Expand All @@ -186,7 +185,7 @@ namespace mfem
{
protected:
VectorCoefficient &uD;
TestLinModel *model;
LinElastMaterialModel *model;
Coefficient *lambda, *mu;
double alpha, kappa;

Expand All @@ -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_) {}

Expand Down
34 changes: 17 additions & 17 deletions src/nlelast_integ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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++)
{
Expand Down Expand Up @@ -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);
}
Expand All @@ -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);

Expand Down Expand Up @@ -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);
}

Expand All @@ -587,15 +587,15 @@ 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);
CalcAdjugate(Tr.Elem1->Jacobian(), adjJ1);
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);

Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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);


Expand Down
2 changes: 1 addition & 1 deletion test/nlelast_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -346,7 +346,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
a1.Assemble();
cout<<"a1.Height() is: "<<a1.Height()<<endl;

TestLinModel model(_mu, _lambda);
LinElastMaterialModel model(_mu, _lambda);
NonlinearForm a2(&fespace);
a2.AddDomainIntegrator(new HyperelasticNLFIntegratorHR(&model));
a2.AddInteriorFaceIntegrator(
Expand Down

0 comments on commit f384ea8

Please sign in to comment.