phynexis-deps/mfem
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

AbsMult for ElasticityInteg

Browse Source
This commit is contained in:
Gabriel Pinochet-Soto
2025-03-28 22:13:22 -07:00
parent 7c077e656d
commit 3b49d70f35
3 changed files with 88 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/integ/bilininteg_elasticity_kernels.cpp
+20
View File
@@ -88,6 +88,26 @@ void ElasticityAddMultPA(const int dim, const int nDofs,
}
}
void ElasticityAddAbsMultPA(const int dim, const int nDofs,
const FiniteElementSpace &fespace, const CoefficientVector &lambda,
const CoefficientVector &mu, const GeometricFactors &geom,
const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y)
{
switch (dim)
{
case 2:
ElasticityAddMultPA_<2>(nDofs, fespace, lambda, mu, geom, maps, x, QVec, y,
true);
break;
case 3:
ElasticityAddMultPA_<3>(nDofs, fespace, lambda, mu, geom, maps, x, QVec, y,
true);
break;
default:
MFEM_ABORT("Only dimensions 2 and 3 supported.");
}
}
void ElasticityAssembleDiagonalPA(const int dim, const int nDofs,
const CoefficientVector &lambda,
const CoefficientVector &mu, const GeometricFactors &geom,
fem/integ/bilininteg_elasticity_kernels.hpp
+39 -9
View File
@@ -67,6 +67,11 @@ void ElasticityAddMultPA(const int dim, const int nDofs,
const CoefficientVector &mu, const GeometricFactors &geom,
const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y);
void ElasticityAddAbsMultPA(const int dim, const int nDofs,
const FiniteElementSpace &fespace, const CoefficientVector &lambda,
const CoefficientVector &mu, const GeometricFactors &geom,
const DofToQuad &maps, const Vector &x, QuadratureFunction &QVec, Vector &y);
/// @brief Elasticity component kernel for AddMultPA.
///
/// Performs y += Ax. Implemented for byNODES ordering only, and does not use
@@ -144,7 +149,7 @@ template<int dim, int i_block = -1, int j_block = -1>
void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
const CoefficientVector &lambda, const CoefficientVector &mu,
const GeometricFactors &geom, const DofToQuad &maps, const Vector &x,
QuadratureFunction &QVec, Vector &y)
QuadratureFunction &QVec, Vector &y, const bool useAbs = false)
{
static_assert((i_block < 0) == (j_block < 0),
"i_block and j_block must both be non-negative or strictly negative.");
@@ -163,8 +168,8 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
ir);
E_To_Q_Map->SetOutputLayout(QVectorLayout::byNODES);
// interpolate physical derivatives to quadrature points.
E_To_Q_Map->PhysDerivatives(x, QVec);
if (!useAbs) { E_To_Q_Map->PhysDerivatives(x, QVec); }
else { E_To_Q_Map->AbsPhysDerivatives(x, QVec); }
const int numPoints = ir.GetNPoints();
const int numEls = fespace.GetNE();
const auto lamDev = Reshape(lambda.Read(), numPoints, numEls);
@@ -206,7 +211,8 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
const int iIndex = isComponent ? 0 : i;
div += gradx(iIndex,i);
}
const real_t w = ipWeights[p] /det(invJ);
const real_t w = useAbs ? std::abs (ipWeights[p] /det(invJ)) :
ipWeights[p] /det(invJ);
for (int m = 0; m < d; m++)
{
for (int q = qLower; q < qUpper; q++)
@@ -220,8 +226,16 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
{
for (int a = 0; a < d; a++)
{
contraction += 2*((a == q)*invJ(m,j_block) + (j_block==q)*invJ(m,a))*(gradx(0,
a));
if (!useAbs)
{
contraction += 2*((a == q)*invJ(m,j_block)
+ (j_block==q)*invJ(m,a))*(gradx(0, a));
}
else
{
contraction += 2*((a == q)*std::abs(invJ(m,j_block))
+ (j_block==q)*std::abs(invJ(m,a)))*(gradx(0, a));
}
}
}
else
@@ -230,15 +244,31 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
{
for (int b = 0; b < d; b++)
{
contraction += ((a == q)*invJ(m,b) + (b==q)*invJ(m,a))
*(gradx(a,b) + gradx(b, a));
if (!useAbs)
{
contraction += ((a == q)*invJ(m,b) + (b==q)*invJ(m,a))
*(gradx(a,b) + gradx(b, a));
}
else
{
contraction += ((a == q)*std::abs(invJ(m,b)) + (b==q)*std::abs(invJ(m,a)))
*(gradx(a,b) + gradx(b, a));
}
}
}
}
// lambda*div(u)*div(v) + 2*mu*sym(grad(u))*sym(grad(v))
// contraction = 4*sym(grad(u))sym(grad(v))
const int qIndex = isComponent ? 0 : q;
Q(p,m,qIndex,e) = w*(lamDev(p, e)*invJ(m,q)*div + 0.5*muDev(p, e)*contraction);
if (!useAbs)
{
Q(p,m,qIndex,e) = w*(lamDev(p, e)*invJ(m,q)*div + 0.5*muDev(p, e)*contraction);
}
else
{
Q(p,m,qIndex,e) = w*(std::abs(lamDev(p, e)*invJ(m,q))*div
+ 0.5*std::abs(muDev(p, e))*contraction);
}
}
}
}
fem/integ/bilininteg_elasticity_pa.cpp
+29
View File
@@ -75,6 +75,35 @@ void ElasticityIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const
AddMultPA(x, y); // Operator is symmetric
}
void ElasticityIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
{
DofToQuad abs_maps;
abs_maps.FE = maps->FE;
abs_maps.IntRule = maps->IntRule;
abs_maps.mode = maps->mode;
abs_maps.ndof = maps->ndof;
abs_maps.nqpt = maps->nqpt;
abs_maps.B = maps->B;
abs_maps.Bt = maps->Bt;
abs_maps.G = maps->G;
abs_maps.Gt = maps->Gt;
abs_maps.B.Abs();
abs_maps.G.Abs();
abs_maps.Bt.Abs();
abs_maps.Gt.Abs();
internal::ElasticityAddAbsMultPA(vdim, ndofs, *fespace, *lambda_quad, *mu_quad,
*geom, abs_maps, x, *q_vec, y);
}
void ElasticityIntegrator::AddAbsMultTransposePA(const Vector &x,
Vector &y) const
{
AddAbsMultPA(x, y); // Operator is symmetric
}
void ElasticityComponentIntegrator::AssemblePA(const FiniteElementSpace &fes)
{
fespace = &fes;