Skip to content

Commit

Permalink
Stiffness matrix handles kappa
Browse files Browse the repository at this point in the history
  • Loading branch information
@larsson4
larsson4 committed Apr 10, 2024
1 parent dd6b463 commit b255718
Show file tree
Hide file tree
Showing 3 changed files with 115 additions and 52 deletions.
4 changes: 2 additions & 2 deletions include/nlelast_integ.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

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

};

Expand Down
53 changes: 41 additions & 12 deletions src/nlelast_integ.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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);

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

Expand All @@ -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);
Expand All @@ -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)
Expand All @@ -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
Expand Down
110 changes: 72 additions & 38 deletions test/nlelast_test.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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,
Expand All @@ -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(
Expand All @@ -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:
Expand All @@ -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)
{
Expand Down Expand Up @@ -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);
Expand Down Expand Up @@ -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: "<<jmatcoef<<endl;


// (1,1) block works
_AssembleBlock(
// (1,1) block
_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 works
// (1,2) block
_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);
shape2, shape2, dshape2_dnM, dshape2_ps, elmat, jmat);
}

// elmat := -elmat + alpha*elmat^t + jmat
elmat *= -1.0;
PrintMatrix(elmat, "checkmat.txt");
if (kappa_is_nonzero)
{
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);
}
}

}
/* PrintMatrix(elmat, "checkmat.txt");
} */
} // namespace mfem
/**
* Simple smoke test to make sure Google Test is properly linked
Expand Down Expand Up @@ -289,7 +316,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
int dim = mesh.Dimension();
int order = 1;
double alpha = 0.0; // IIPG
double kappa = 0.0; // TODO: Enable kappa = -1.0
double kappa = -1.0; // TODO: Enable kappa = -1.0
DG_FECollection fec(order, dim, BasisType::GaussLobatto);
FiniteElementSpace fespace(&mesh, &fec, dim);

Expand All @@ -310,28 +337,28 @@ 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);
new Test_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??
//a2.AddBdrFaceIntegrator(
// new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr);
//a2.AddInteriorFaceIntegrator(
// new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa)); //Needed??
a2.AddBdrFaceIntegrator(
new DGElasticityIntegrator(lambda_c, mu_c, alpha, kappa), dir_bdr);
a2.Assemble(); */

// Create vectors to hold the values of the forms
Expand Down Expand Up @@ -366,8 +393,8 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
double norm_diff = 0.0;
cout << "Linear residual norm: " << y1.Norml2() << endl;
cout << "Nonlinear residual norm: " << y2.Norml2() << endl;

/* cout << "print y1: "<< endl;
/*
cout << "print y1: "<< endl;
for (size_t i = 0; i < y1.Size(); i++)
{
cout<<y1[i]<<endl;
Expand All @@ -378,7 +405,14 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)
{
cout<<y2[i]<<endl;
}
*/
cout << "print y1/y2: "<< endl;
for (size_t i = 0; i < y2.Size(); i++)
{
cout<<y1[i]/y2[i]<<endl;
} */


y1 -= y2;
norm_diff = y1.Norml2();
cout << "Residual difference norm: " << norm_diff << endl;
Expand Down Expand Up @@ -427,7 +461,7 @@ TEST(TempLinStiffnessMatrices, Test_NLElast)

cout << "Linear RHS norm: " << b1.Norml2() << endl;
cout << "Nonlinear RHS norm: " << b2.Norml2() << endl;

b1 -= b2;
norm_diff = b1.Norml2();

Expand Down

0 comments on commit b255718

Please sign in to comment.