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

Add more virtual functions for different restrictions. Missing face/bdr restr.

Browse Source
This commit is contained in:
Gabriel Esteban Pinochet Soto
2024-08-16 11:09:00 -07:00
parent 0b5aa139ac
commit 29f43a71fa
5 changed files with 91 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files
fem/bilinearform_ext.cpp
+22 -6
View File
@@ -798,11 +798,11 @@ void PABilinearFormExtension::AbsMult(const Vector &x, Vector &y) const
{
if (integrators[i]->Patchwise())
{
integrators[i]->AddMultNURBSPA(x, y);
integrators[i]->AddAbsMultNURBSPA(x, y);
}
else
{
integrators[i]->AddMultPA(x, y);
integrators[i]->AddAbsMultPA(x, y);
}
}
}
@@ -811,17 +811,32 @@ void PABilinearFormExtension::AbsMult(const Vector &x, Vector &y) const
if (iSz)
{
Array<Array<int>*> &elem_markers = *a->GetDBFI_Marker();
auto el_rest = dynamic_cast<const ElementRestriction*>(elem_restrict);
MFEM_VERIFY(el_rest, "elem_restrict is not ElementRestriction*!");
el_rest->AbsMult(x, localX);
auto H1elem_restrict = dynamic_cast<const ElementRestriction*>(elem_restrict);
if (H1elem_restrict)
{
H1elem_restrict->AbsMult(x, localX);
}
else
{
elem_restrict->Mult(x,localX);
}
localY = 0.0;
for (int i = 0; i < iSz; ++i)
{
AddAbsMultWithMarkers(*integrators[i], localX, elem_markers[i],
elem_attributes, false, localY);
}
el_rest->AbsMultTranspose(localY, y);
if (H1elem_restrict)
{
H1elem_restrict->AbsMultTranspose(x, localX);
}
else
{
elem_restrict->MultTranspose(x,localX);
}
}
else
{
@@ -854,6 +869,7 @@ void PABilinearFormExtension::AbsMult(const Vector &x, Vector &y) const
int_face_restrict_lex->Mult(*x_dg, int_face_X);
if (int_face_dXdn.Size() > 0)
{
MFEM_ABORT("NormalDerivativeAbsMult is not implemented...");
int_face_restrict_lex->NormalDerivativeMult(*x_dg, int_face_dXdn);
}
if (int_face_X.Size() > 0)
fem/bilininteg.cpp
+19
View File
@@ -117,6 +117,12 @@ void BilinearFormIntegrator::AddMultNURBSPA(const Vector &, Vector &) const
" is not implemented for this class.");
}
void BilinearFormIntegrator::AddAbsMultNURBSPA(const Vector &, Vector &) const
{
MFEM_ABORT("BilinearFormIntegrator::AddAbsMultNURBSPA(...)\n"
" is not implemented for this class.");
}
void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
{
MFEM_ABORT("BilinearFormIntegrator::AddMultTransposePA(...)\n"
@@ -142,12 +148,25 @@ void BilinearFormIntegrator::AddMultMF(const Vector &, Vector &) const
" is not implemented for this class.");
}
void BilinearFormIntegrator::AddAbsMultMF(const Vector &, Vector &) const
{
MFEM_ABORT("BilinearFormIntegrator::AddAbsMultMF(...)\n"
" is not implemented for this class.");
}
void BilinearFormIntegrator::AddMultTransposeMF(const Vector &, Vector &) const
{
MFEM_ABORT("BilinearFormIntegrator::AddMultTransposeMF(...)\n"
" is not implemented for this class.");
}
void BilinearFormIntegrator::AddAbsMultTransposeMF(const Vector &,
Vector &) const
{
MFEM_ABORT("BilinearFormIntegrator::AddAbsMultTransposeMF(...)\n"
" is not implemented for this class.");
}
void BilinearFormIntegrator::AssembleDiagonalMF(Vector &)
{
MFEM_ABORT("BilinearFormIntegrator::AssembleDiagonalMF(...)\n"
fem/bilininteg.hpp
+6
View File
@@ -81,6 +81,8 @@ public:
/// Method for partially assembled action on NURBS patches.
virtual void AddMultNURBSPA(const Vector&x, Vector&y) const;
virtual void AddAbsMultNURBSPA(const Vector&x, Vector&y) const;
/// Method for partially assembled transposed action.
/** Perform the transpose action of integrator on the input @a x and add the
result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
@@ -115,6 +117,8 @@ public:
called. */
virtual void AddMultMF(const Vector &x, Vector &y) const;
virtual void AddAbsMultMF(const Vector &x, Vector &y) const;
/** Perform the transpose action of integrator on the input @a x and add the
result to the output @a y. Both @a x and @a y are E-vectors, i.e. they
represent the element-wise discontinuous version of the FE space.
@@ -123,6 +127,8 @@ public:
called. */
virtual void AddMultTransposeMF(const Vector &x, Vector &y) const;
virtual void AddAbsMultTransposeMF(const Vector &x, Vector &y) const;
/// Assemble diagonal and add it to Vector @a diag.
virtual void AssembleDiagonalMF(Vector &diag);
fem/restriction.cpp
+23
View File
@@ -678,6 +678,29 @@ void ConformingFaceRestriction::Mult(const Vector& x, Vector& y) const
});
}
void ConformingFaceRestriction::MultUnsigned(const Vector& x, Vector& y) const
{
if (nf==0) { return; }
// Assumes all elements have the same number of dofs
const int nface_dofs = face_dofs;
const int vd = vdim;
const bool t = byvdim;
auto d_indices = scatter_indices.Read();
auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
auto d_y = Reshape(y.Write(), nface_dofs, vd, nf);
mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
{
const int s_idx = d_indices[i];
const int idx = (s_idx >= 0) ? s_idx : -1 - s_idx;
const int dof = i % nface_dofs;
const int face = i / nface_dofs;
for (int c = 0; c < vd; ++c)
{
d_y(dof, c, face) = d_x(t?c:idx, t?idx:c);
}
});
}
static void ConformingFaceRestriction_AddMultTranspose(
const int ndofs,
const int face_dofs,
fem/restriction.hpp
+21
View File
@@ -181,6 +181,13 @@ public:
*/
void Mult(const Vector &x, Vector &y) const override = 0;
virtual void MultUnsigned(const Vector &x, Vector &y) const
{
MFEM_ABORT("MultUnsigned not implemented yet...");
}
void AbsMult(const Vector &x, Vector &y) const override { MultUnsigned(x,y); }
/** @brief Add the face degrees of freedom @a x to the element degrees of
freedom @a y.
@@ -229,6 +236,12 @@ public:
AddMultTranspose(x, y);
}
void AbsMultTranspose(const Vector &x, Vector &y) const override
{
y = 0.0;
AddMultTransposeUnsigned(x,y);
}
/** @brief For each face, sets @a y to the partial derivative of @a x with
respect to the reference coordinate whose direction is
perpendicular to the face on the reference element.
@@ -325,6 +338,8 @@ public:
ElementDofOrdering. */
void Mult(const Vector &x, Vector &y) const override;
void MultUnsigned(const Vector &x, Vector &y) const override;
using FaceRestriction::AddMultTransposeInPlace;
/** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
@@ -348,6 +363,12 @@ public:
void AddMultTransposeUnsigned(const Vector &x, Vector &y,
const real_t a = 1.0) const override;
void AddAbsMultTranspose(const Vector &x, Vector &y) const
{
y = 0.0;
AddMultTransposeUnsigned(x,y);
}
private:
/** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
for the gathering: E-vector to L-vector.