Skip to content

Commit

Permalink
Improve high-order node transfer for mesh splitting
Browse files Browse the repository at this point in the history
  • Loading branch information
@sebastiangrimberg
sebastiangrimberg committed Jun 21, 2024
1 parent e2810a0 commit 440532b
Show file tree
Hide file tree
Showing 3 changed files with 234 additions and 59 deletions.
102 changes: 68 additions & 34 deletions palace/fem/interpolator.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,14 @@
namespace palace
{

namespace
{

constexpr auto GSLIB_BB_TOL = 0.01; // MFEM defaults, slightly reduced bounding box
constexpr auto GSLIB_NEWTON_TOL = 1.0e-12;

} // namespace

#if defined(MFEM_USE_GSLIB)
InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh &mesh)
: op(mesh.GetComm())
Expand All @@ -24,8 +32,6 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh
{
return;
}
constexpr double bb_t = 0.1; // MFEM defaults
constexpr double newton_tol = 1.0e-12;
const int dim = mesh.SpaceDimension();
MFEM_VERIFY(
mesh.Dimension() == dim,
Expand All @@ -36,16 +42,14 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh
int i = 0;
for (const auto &[idx, data] : iodata.domains.postpro.probe)
{
// Use default ordering byNODES.
xyz(i) = data.center[0];
xyz(npts + i) = data.center[1];
if (dim == 3)
for (int d = 0; d < dim; d++)
{
xyz(2 * npts + i) = data.center[2];
// Use default ordering byNODES.
xyz(d * npts + i) = data.center[d];
}
op_idx[i++] = idx;
}
op.Setup(mesh, bb_t, newton_tol, npts);
op.Setup(mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);
op.FindPoints(xyz, mfem::Ordering::byNODES);
op.SetDefaultInterpolationValue(0.0);
i = 0;
Expand All @@ -70,18 +74,26 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh
std::vector<double> InterpolationOperator::ProbeField(const mfem::ParGridFunction &U)
{
#if defined(MFEM_USE_GSLIB)
// Interpolated vector values are returned from GSLIB interpolator byNODES, which we
// transform to byVDIM for output.
// Interpolated vector values are returned from GSLIB interpolator with the same ordering
// as the source grid function, which we transform to byVDIM for output.
const int npts = op.GetCode().Size();
const int vdim = U.VectorDim();
std::vector<double> vals(npts * vdim);
mfem::Vector v(npts * vdim);
op.Interpolate(U, v);
for (int d = 0; d < vdim; d++)
if (U.FESpace()->GetOrdering() == mfem::Ordering::byVDIM)
{
mfem::Vector v(vals.data(), npts * vdim);
op.Interpolate(U, v);
}
else
{
for (int i = 0; i < npts; i++)
mfem::Vector v(npts * vdim);
op.Interpolate(U, v);
for (int d = 0; d < vdim; d++)
{
vals[i * vdim + d] = v(d * npts + i);
for (int i = 0; i < npts; i++)
{
vals[i * vdim + d] = v(d * npts + i);
}
}
}
return vals;
Expand Down Expand Up @@ -147,16 +159,14 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
for (int i = 0; i < dest_mesh.GetNE(); i++)
{
const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i);
mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
T.Transform(fe.GetNodes(), pointmat);
for (int j = 0; j < pointmat.Width(); j++)
{
// Fill with ordering byNODES.
xyz(offset + j) = pointmat(0, j);
xyz(npts + offset + j) = pointmat(1, j);
if (dim == 3)
for (int d = 0; d < dim; d++)
{
xyz(2 * npts + offset + j) = pointmat(2, j);
// Use default ordering byNODES.
xyz(d * npts + offset + j) = pointmat(d, j);
}
}
offset += pointmat.Width();
Expand All @@ -171,13 +181,11 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
mfem::FindPointsGSLIB op(comm);
constexpr double bb_t = 0.1; // MFEM defaults
constexpr double newton_tol = 1.0e-12;
op.Setup(src_mesh, bb_t, newton_tol, npts);
op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);

// Perform the interpolation and fill the target GridFunction (see MFEM's field-interp
// miniapp).
const int vdim = V.VectorDim();
const int vdim = U.VectorDim();
mfem::Vector vals(npts * vdim);
op.SetDefaultInterpolationValue(0.0);
op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
Expand All @@ -197,24 +205,24 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
const int elem_npts = fe.GetNodes().GetNPoints();
elem_vals.SetSize(elem_npts * vdim);
for (int j = 0; j < elem_npts; j++)
for (int d = 0; d < vdim; d++)
{
for (int d = 0; d < vdim; d++)
for (int j = 0; j < elem_npts; j++)
{
// Arrange element values byNODES to align with GetElementVDofs.
int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
? d * npts + offset + j
: (offset + j) * dim + d;
: (offset + j) * vdim + d;
elem_vals(d * elem_npts + j) = vals(idx);
}
}
offset += elem_npts;
const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
if (dof_trans)
{
dof_trans->TransformPrimal(elem_vals);
}
V.SetSubVector(vdofs, elem_vals);
offset += elem_npts;
}
}
else
Expand All @@ -234,21 +242,20 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
for (int i = 0; i < dest_mesh.GetNE(); i++)
{
const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i);
mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i);
mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i);
const int elem_npts = fe.GetNodes().GetNPoints();
elem_vals.SetSize(elem_npts * vdim);
for (int j = 0; j < elem_npts; j++)
for (int d = 0; d < vdim; d++)
{
for (int d = 0; d < vdim; d++)
for (int j = 0; j < elem_npts; j++)
{
// Arrange element values byVDIM for ProjectFromNodes.
int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES)
? d * npts + offset + j
: (offset + j) * dim + d;
: (offset + j) * vdim + d;
elem_vals(j * vdim + d) = vals(idx);
}
}
offset += elem_npts;
const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs);
v.SetSize(vdofs.Size());
fe.ProjectFromNodes(elem_vals, T, v);
Expand All @@ -257,13 +264,40 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V)
dof_trans->TransformPrimal(v);
}
V.SetSubVector(vdofs, v);
offset += elem_npts;
}
}
#else
MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
#endif
}

void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U,
mfem::Vector &vals, mfem::Ordering::Type ordering)
{
#if defined(MFEM_USE_GSLIB)
// Set up the interpolator.
auto &src_mesh = *U.FESpace()->GetMesh();
MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!");
const int dim = src_mesh.SpaceDimension();
const int npts = xyz.Size() / dim;
auto *src_pmesh = dynamic_cast<mfem::ParMesh *>(&src_mesh);
MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF;
mfem::FindPointsGSLIB op(comm);
op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts);

// Perform the interpolation, with the ordering of the returned values matching the
// ordering of the source grid function.
const int vdim = U.VectorDim();
MFEM_VERIFY(vals.Size() == npts * vdim, "Incorrect size for interpolated values vector!");
op.SetDefaultInterpolationValue(0.0);
op.SetL2AvgType(mfem::FindPointsGSLIB::NONE);
op.Interpolate(xyz, U, vals, ordering);
#else
MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!");
#endif
}

} // namespace fem

} // namespace palace
6 changes: 6 additions & 0 deletions palace/fem/interpolator.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,6 +42,12 @@ namespace fem
// Similar to MFEM's field-interp miniapp.
void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V);

// Interpolate a function at a specific list of points, specified using the provided
// ordering. The output vector values are always arranged byVDIM.
void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U,
mfem::Vector &V,
mfem::Ordering::Type ordering = mfem::Ordering::byNODES);

} // namespace fem

} // namespace palace
Expand Down
Loading

0 comments on commit 440532b

Please sign in to comment.