Skip to content

Commit

Permalink
GetSharedFaceTransformationsByLocalIndex can be used to replace GetSh…
Browse files Browse the repository at this point in the history 
…aredFaceTransformations and the need for passing around a local to shared face map
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Jan 5, 2024
1 parent 9b04b4f commit 38375a1
Show file tree
Hide file tree
Showing 4 changed files with 27 additions and 75 deletions.
12 changes: 5 additions & 7 deletions palace/fem/coefficient.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,9 +7,9 @@ namespace palace
{

void BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(
int i, const mfem::ParMesh &mesh, const std::unordered_map<int, int> &local_to_shared,
mfem::FaceElementTransformations &FET, mfem::IsoparametricTransformation &T1,
mfem::IsoparametricTransformation &T2, const mfem::IntegrationPoint *ip)
int i, const mfem::ParMesh &mesh, mfem::FaceElementTransformations &FET,
mfem::IsoparametricTransformation &T1, mfem::IsoparametricTransformation &T2,
const mfem::IntegrationPoint *ip)
{
// Return transformations for elements attached to the given boundary element. FET.Elem1
// always exists but FET.Elem2 may not if the element is truly a single-sided boundary.
Expand All @@ -25,8 +25,7 @@ void BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(
if (info2 >= 0 && iel2 < 0)
{
// Face is shared with another subdomain.
const int &ishared = local_to_shared.at(f);
mesh.GetSharedFaceTransformations(ishared, &FET, &T1, &T2);
mesh.GetSharedFaceTransformationsByLocalIndex(f, &FET, &T1, &T2);
}
else
{
Expand All @@ -51,8 +50,7 @@ void BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(
// too.
MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
"Unexpected element type in BdrGridFunctionCoefficient!");
GetBdrElementNeighborTransformations(T.ElementNo, mesh, local_to_shared, FET, T1, T2,
&ip);
GetBdrElementNeighborTransformations(T.ElementNo, mesh, FET, T1, T2, &ip);

// If desired, get vector pointing from center of boundary element into element 1 for
// orientations.
Expand Down
40 changes: 13 additions & 27 deletions palace/fem/coefficient.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@

#include <complex>
#include <memory>
#include <unordered_map>
#include <utility>
#include <vector>
#include <mfem.hpp>
Expand All @@ -31,7 +30,6 @@ class BdrGridFunctionCoefficient
// XX TODO: For thread-safety (multiple threads evaluating a coefficient simultaneously),
// the FET, FET.Elem1, and FET.Elem2 objects cannot be shared
const mfem::ParMesh &mesh;
const std::unordered_map<int, int> &local_to_shared;
mfem::FaceElementTransformations FET;
mfem::IsoparametricTransformation T1, T2, TF;

Expand All @@ -40,20 +38,16 @@ class BdrGridFunctionCoefficient
mfem::Vector *C1 = nullptr);

public:
BdrGridFunctionCoefficient(const mfem::ParMesh &mesh,
const std::unordered_map<int, int> &local_to_shared)
: mesh(mesh), local_to_shared(local_to_shared)
{
}
BdrGridFunctionCoefficient(const mfem::ParMesh &mesh) : mesh(mesh) {}

// For a boundary element, return the element transformation objects for the neighboring
// domain elements. FET.Elem2 may be nullptr if the boundary is a true one-sided boundary,
// but if it is shared with another subdomain then it will be populated. Expects
// ParMesh::ExchangeFaceNbrData has been called already.
static void GetBdrElementNeighborTransformations(
int i, const mfem::ParMesh &mesh, const std::unordered_map<int, int> &local_to_shared,
mfem::FaceElementTransformations &FET, mfem::IsoparametricTransformation &T1,
mfem::IsoparametricTransformation &T2, const mfem::IntegrationPoint *ip = nullptr);
int i, const mfem::ParMesh &mesh, mfem::FaceElementTransformations &FET,
mfem::IsoparametricTransformation &T1, mfem::IsoparametricTransformation &T2,
const mfem::IntegrationPoint *ip = nullptr);

// Return normal vector to the boundary element at an integration point (it is assumed
// that the element transformation has already been configured at the integration point of
Expand Down Expand Up @@ -81,10 +75,9 @@ class BdrCurrentVectorCoefficient : public mfem::VectorCoefficient,
BdrCurrentVectorCoefficient(const mfem::ParGridFunction &gf,
const MaterialOperator &mat_op)
: mfem::VectorCoefficient(mat_op.SpaceDimension()),
BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
B(gf), mat_op(mat_op), C1(gf.VectorDim()), W(gf.VectorDim()), VU(gf.VectorDim()),
VL(gf.VectorDim()), nor(gf.VectorDim())
BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()), B(gf), mat_op(mat_op),
C1(gf.VectorDim()), W(gf.VectorDim()), VU(gf.VectorDim()), VL(gf.VectorDim()),
nor(gf.VectorDim())
{
}

Expand Down Expand Up @@ -137,8 +130,7 @@ class BdrChargeCoefficient : public mfem::Coefficient, public BdrGridFunctionCoe

public:
BdrChargeCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op)
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
E(gf), mat_op(mat_op), C1(gf.VectorDim()), W(gf.VectorDim()), VU(gf.VectorDim()),
VL(gf.VectorDim()), nor(gf.VectorDim())
{
Expand Down Expand Up @@ -179,8 +171,7 @@ class BdrFluxCoefficient : public mfem::Coefficient, public BdrGridFunctionCoeff
public:
BdrFluxCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op,
const mfem::Vector &d)
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
B(gf), dir(d), V(gf.VectorDim()), VL(gf.VectorDim()), nor(gf.VectorDim())
{
}
Expand Down Expand Up @@ -275,8 +266,7 @@ class DielectricInterfaceCoefficient : public mfem::Coefficient,
DielectricInterfaceCoefficient(const mfem::ParGridFunction &gf,
const MaterialOperator &mat_op, double ti, double ei,
const mfem::Vector &s)
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
E(gf), mat_op(mat_op), ts(ti), epsilon(ei), side(s), C1(gf.VectorDim()),
V(gf.VectorDim()), nor(gf.VectorDim())
{
Expand Down Expand Up @@ -365,8 +355,7 @@ class EnergyDensityCoefficient : public mfem::Coefficient, public BdrGridFunctio

public:
EnergyDensityCoefficient(const GridFunctionType &gf, const MaterialOperator &mat_op)
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
U(gf), mat_op(mat_op), V(mat_op.SpaceDimension())
{
}
Expand Down Expand Up @@ -460,9 +449,7 @@ class BdrFieldVectorCoefficient : public mfem::VectorCoefficient,
public:
BdrFieldVectorCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op)
: mfem::VectorCoefficient(mat_op.SpaceDimension()),
BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
U(gf), mat_op(mat_op)
BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()), U(gf), mat_op(mat_op)
{
}

Expand Down Expand Up @@ -494,8 +481,7 @@ class BdrFieldCoefficient : public mfem::Coefficient, public BdrGridFunctionCoef

public:
BdrFieldCoefficient(const mfem::ParGridFunction &gf, const MaterialOperator &mat_op)
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh(),
mat_op.GetLocalToSharedFaceMap()),
: mfem::Coefficient(), BdrGridFunctionCoefficient(*gf.ParFESpace()->GetParMesh()),
U(gf), mat_op(mat_op)
{
}
Expand Down
41 changes: 7 additions & 34 deletions palace/models/materialoperator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -278,20 +278,6 @@ mfem::DenseMatrix ToDenseMatrix(const config::SymmetricMatrixData<N> &data)
return M;
}

auto BuildLocalToSharedFaceMap(const mfem::ParMesh &mesh)
{
// Construct shared face mapping for boundary coefficients. The inverse mapping is
// constructed as part of mfem::ParMesh, but we need this mapping when looping over
// all mesh faces.
std::unordered_map<int, int> l2s;
l2s.reserve(mesh.GetNSharedFaces());
for (int i = 0; i < mesh.GetNSharedFaces(); i++)
{
l2s[mesh.GetSharedFace(i)] = i;
}
return l2s;
}

auto BuildAttributeGlobalToLocal(const mfem::ParMesh &mesh)
{
// Set up sparse map from global domain attributes to local ones on this process.
Expand Down Expand Up @@ -327,25 +313,18 @@ auto BuildAttributeGlobalToLocal(const mfem::ParMesh &mesh)
}

auto GetBdrNeighborAttribute(int i, const mfem::ParMesh &mesh,
const std::unordered_map<int, int> &face_loc_to_shared,
mfem::FaceElementTransformations &FET,
mfem::IsoparametricTransformation &T1,
mfem::IsoparametricTransformation &T2)
{
// For internal boundaries, use the element which corresponds to the vacuum domain, or
// at least the one with the higher speed of light.
BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(
i, mesh, face_loc_to_shared, FET, T1, T2);
// return (FET.Elem2 && GetLightSpeedMin(FET.Elem2->Attribute) >
// GetLightSpeedMax(FET.Elem1->Attribute))
// ? FET.Elem2->Attribute
// : FET.Elem1->Attribute;
BdrGridFunctionCoefficient::GetBdrElementNeighborTransformations(i, mesh, FET, T1, T2);
return (FET.Elem2 && FET.Elem2->Attribute < FET.Elem1->Attribute) ? FET.Elem2->Attribute
: FET.Elem1->Attribute;
}

auto BuildBdrAttributeGlobalToLocal(const mfem::ParMesh &mesh,
const std::unordered_map<int, int> &face_loc_to_shared)
auto BuildBdrAttributeGlobalToLocal(const mfem::ParMesh &mesh)
{
// Set up sparse map from global boundary attributes to local ones on this process. Each
// original global boundary attribute maps to a key-value pairing of global domain
Expand All @@ -357,7 +336,7 @@ auto BuildBdrAttributeGlobalToLocal(const mfem::ParMesh &mesh,
for (int i = 0; i < mesh.GetNBE(); i++)
{
const int attr = mesh.GetBdrAttribute(i);
const int nbr_attr = GetBdrNeighborAttribute(i, mesh, face_loc_to_shared, FET, T1, T2);
const int nbr_attr = GetBdrNeighborAttribute(i, mesh, FET, T1, T2);
auto &bdr_attr_map = loc_bdr_attr[attr];
if (bdr_attr_map.find(nbr_attr) == bdr_attr_map.end())
{
Expand All @@ -372,9 +351,8 @@ auto BuildBdrAttributeGlobalToLocal(const mfem::ParMesh &mesh,
MaterialOperator::MaterialOperator(const IoData &iodata, mfem::ParMesh &mesh) : mesh(mesh)
{
mesh.ExchangeFaceNbrData();
face_loc_to_shared = BuildLocalToSharedFaceMap(mesh);
loc_attr = BuildAttributeGlobalToLocal(mesh);
loc_bdr_attr = BuildBdrAttributeGlobalToLocal(mesh, face_loc_to_shared);
loc_bdr_attr = BuildBdrAttributeGlobalToLocal(mesh);

SetUpMaterialProperties(iodata, mesh);
}
Expand Down Expand Up @@ -604,26 +582,21 @@ int MaterialOperator::GetAttributeGlobalToLocal(mfem::ElementTransformation &T)
"Invalid domain attribute " << T.Attribute << "!");
const int nbr_attr = [&]()
{
// XX TODO INCORRECT FOR H-MULTIGRID: T.ElementNo SHOULD BE USED TO FIND THE MESH
// NEIGHBOR ON THE COARSE MESH

mfem::FaceElementTransformations FET; // XX TODO: Preallocate these for all elements
mfem::IsoparametricTransformation T1, T2;
if (const auto *submesh = dynamic_cast<const mfem::ParSubMesh *>(T.mesh))
{
MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::ELEMENT,
"Unexpected element type in GetAttributeGlobalToLocal!");
return GetBdrNeighborAttribute(submesh->GetParentElementIDMap()[T.ElementNo],
*submesh->GetParent(), face_loc_to_shared, FET, T1,
T2);
*submesh->GetParent(), FET, T1, T2);
}
else
{
MFEM_ASSERT(T.ElementType == mfem::ElementTransformation::BDR_ELEMENT,
"Unexpected element type in GetAttributeGlobalToLocal!");
return GetBdrNeighborAttribute(T.ElementNo,
*static_cast<const mfem::ParMesh *>(T.mesh),
face_loc_to_shared, FET, T1, T2);
return GetBdrNeighborAttribute(
T.ElementNo, *static_cast<const mfem::ParMesh *>(T.mesh), FET, T1, T2);
}
}();
auto it = bdr_attr_map->second.find(nbr_attr);
Expand Down
9 changes: 2 additions & 7 deletions palace/models/materialoperator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,15 +35,12 @@ class MaterialOperator
// penetration depth.
mfem::Array<int> losstan_attr, conductivity_attr, london_attr;

// Shared face mapping for boundary coefficients.
std::unordered_map<int, int> face_loc_to_shared;

// Attribute mapping for (global, 1-based) domain and boundary attributes to those on this
// process (still 1-based). For boundaries, the inner map is a mapping from neighboring
// domain attribute to the resulting local boundary attribute (to discern boundary
// elements with global boundary attribute which borders more than one domain). Interior
// boundaries use as neighbor the element which corresponds to the vacuum domain, or at
// least the one with the higher speed of light.
// boundaries use as neighbor the element with the smaller domain attribute in order to
// be consistent when the interior boundary element normals are not aligned.
std::unordered_map<int, int> loc_attr;
std::unordered_map<int, std::unordered_map<int, int>> loc_bdr_attr;

Expand Down Expand Up @@ -96,8 +93,6 @@ class MaterialOperator
const auto &GetAttributeToMaterial() const { return attr_mat; }
mfem::Array<int> GetBdrAttributeToMaterial() const;

const auto &GetLocalToSharedFaceMap() const { return face_loc_to_shared; }

const auto &GetAttributeGlobalToLocal() const { return loc_attr; }

const auto &GetBdrAttributeGlobalToLocal() const { return loc_bdr_attr; }
Expand Down

0 comments on commit 38375a1

Please sign in to comment.