Skip to content

Commit

Permalink
Fix bug where projection onto internal boundaries fails to identify p…
Browse files Browse the repository at this point in the history 
…arent dofs in NC faces
  • Loading branch information
@hughcars
hughcars committed Apr 23, 2024
1 parent 31bd94e commit 2dabf82
Show file tree
Hide file tree
Showing 2 changed files with 119 additions and 36 deletions.
97 changes: 61 additions & 36 deletions fem/gridfunc.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2061,41 +2061,37 @@ void GridFunction::AccumulateAndCountBdrValues(
Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
Array<int> &values_counter)
{
int i, j, fdof, d, ind, vdim;
real_t val;
const FiniteElement *fe;
ElementTransformation *transf;
Array<int> vdofs;
Vector vc;

values_counter.SetSize(Size());
values_counter = 0;

vdim = fes->GetVDim();

const int vdim = fes->GetVDim();
HostReadWrite();

for (i = 0; i < fes->GetNBE(); i++)
for (int i = 0; i < fes->GetNBE(); i++)
{
if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }

fe = fes->GetBE(i);
fdof = fe->GetDof();
transf = fes->GetBdrElementTransformation(i);
const FiniteElement *fe = fes->GetBE(i);
const int fdof = fe->GetDof();
ElementTransformation *transf = fes->GetBdrElementTransformation(i);
const IntegrationRule &ir = fe->GetNodes();
fes->GetBdrElementVDofs(i, vdofs);

for (j = 0; j < fdof; j++)
for (int j = 0; j < fdof; j++)
{
const IntegrationPoint &ip = ir.IntPoint(j);
transf->SetIntPoint(&ip);
if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
for (d = 0; d < vdim; d++)
for (int d = 0; d < vdim; d++)
{
if (!vcoeff && !coeff[d]) { continue; }

val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
if ( (ind = vdofs[fdof*d+j]) < 0 )
real_t val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
int ind = vdofs[fdof*d+j];
if ( ind < 0 )
{
val = -val, ind = -1-ind;
}
Expand All @@ -2117,37 +2113,31 @@ void GridFunction::AccumulateAndCountBdrValues(
// iff A_ij != 0. It is sufficient to resolve just the first level of
// dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
// Cases like these arise in 3D when boundary edges are constrained by
// (depend on) internal faces/elements. We use the virtual method
// GetBoundaryClosure from NCMesh to resolve the dependencies.

if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
// (depend on) internal faces/elements, or for internal boundaries in 2 or
// 3D. We use the virtual method GetBoundaryClosure from NCMesh to resolve
// the dependencies.
if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
fes->GetMesh()->Dimension() == 3))
{
Vector vals;
Mesh *mesh = fes->GetMesh();
NCMesh *ncmesh = mesh->ncmesh;
Array<int> bdr_edges, bdr_vertices, bdr_faces;
ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges, bdr_faces);

for (i = 0; i < bdr_edges.Size(); i++)
auto mark_dofs = [&](ElementTransformation &transf, const FiniteElement &fe)
{
int edge = bdr_edges[i];
fes->GetEdgeVDofs(edge, vdofs);
if (vdofs.Size() == 0) { continue; }

transf = mesh->GetEdgeTransformation(edge);
transf->Attribute = -1; // TODO: set the boundary attribute
fe = fes->GetEdgeElement(edge);
if (!vcoeff)
{
vals.SetSize(fe->GetDof());
for (d = 0; d < vdim; d++)
vals.SetSize(fe.GetDof());
for (int d = 0; d < vdim; d++)
{
if (!coeff[d]) { continue; }

fe->Project(*coeff[d], *transf, vals);
fe.Project(*coeff[d], transf, vals);
for (int k = 0; k < vals.Size(); k++)
{
ind = vdofs[d*vals.Size()+k];
const int ind = vdofs[d*vals.Size()+k];
if (++values_counter[ind] == 1)
{
(*this)(ind) = vals(k);
Expand All @@ -2161,11 +2151,11 @@ void GridFunction::AccumulateAndCountBdrValues(
}
else // vcoeff != NULL
{
vals.SetSize(vdim*fe->GetDof());
fe->Project(*vcoeff, *transf, vals);
vals.SetSize(vdim*fe.GetDof());
fe.Project(*vcoeff, transf, vals);
for (int k = 0; k < vals.Size(); k++)
{
ind = vdofs[k];
const int ind = vdofs[k];
if (++values_counter[ind] == 1)
{
(*this)(ind) = vals(k);
Expand All @@ -2176,6 +2166,28 @@ void GridFunction::AccumulateAndCountBdrValues(
}
}
}
};

for (auto edge : bdr_edges)
{
fes->GetEdgeVDofs(edge, vdofs);
if (vdofs.Size() == 0) { continue; }

ElementTransformation *transf = mesh->GetEdgeTransformation(edge);
transf->Attribute = -1; // TODO: set the boundary attribute
const FiniteElement *fe = fes->GetEdgeElement(edge);
mark_dofs(*transf, *fe);
}

for (auto face : bdr_faces)
{
fes->GetFaceDofs(face, vdofs);
if (vdofs.Size() == 0) { continue; }

ElementTransformation *transf = mesh->GetFaceTransformation(face);
transf->Attribute = -1; // TODO: set the boundary attribute
const FiniteElement *fe = fes->GetFaceElement(face);
mark_dofs(*transf, *fe);
}
}
}
Expand Down Expand Up @@ -2228,16 +2240,16 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
accumulate_dofs(dofs, lvec, *this, values_counter);
}

if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
fes->GetMesh()->Dimension() == 3))
{
Mesh *mesh = fes->GetMesh();
NCMesh *ncmesh = mesh->ncmesh;
Array<int> bdr_edges, bdr_vertices, bdr_faces;
ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges, bdr_faces);

for (int i = 0; i < bdr_edges.Size(); i++)
for (auto edge : bdr_edges)
{
int edge = bdr_edges[i];
fes->GetEdgeDofs(edge, dofs);
if (dofs.Size() == 0) { continue; }

Expand All @@ -2248,6 +2260,19 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
fe->Project(vcoeff, *T, lvec);
accumulate_dofs(dofs, lvec, *this, values_counter);
}

for (auto face : bdr_faces)
{
fes->GetFaceDofs(face, dofs);
if (dofs.Size() == 0) { continue; }

T = mesh->GetFaceTransformation(face);
T->Attribute = -1; // TODO: set the boundary attribute
fe = fes->GetFaceElement(face);
lvec.SetSize(fe->GetDof());
fe->Project(vcoeff, *T, lvec);
accumulate_dofs(dofs, lvec, *this, values_counter);
}
}
}

Expand Down
58 changes: 58 additions & 0 deletions tests/unit/mesh/test_ncmesh.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1691,6 +1691,64 @@ TEST_CASE("Parallel RP=I", "[Parallel], [NCMesh]")
}
}

TEST_CASE("InternalBoundaryProjectBdrCoefficient", "[Parallel], [NCMesh]")
{
auto test_project_H1 = [](ParMesh &mesh, int order, double coef)
{
H1_FECollection fe_collection(order, mesh.SpaceDimension());
ParFiniteElementSpace fe_space(&mesh, &fe_collection);
ParGridFunction x(&fe_space);
x = -coef;
ConstantCoefficient c(coef);

// Check projecting on the internal face sets essential dof.
Array<int> ess_bdr(mesh.bdr_attributes.Max());
ess_bdr = 0;
ess_bdr.Last() = 1; // internal boundary
x.ProjectBdrCoefficient(c, ess_bdr);

Array<int> ess_tdofs_list;
fe_space.GetEssentialTrueDofs(ess_bdr, ess_tdofs_list);

const auto &tvec = x.GetTrueVector();
for (auto ess_dof : ess_tdofs_list)
{
CHECK(tvec[ess_dof] == Approx(coef).epsilon(1e-8));
}

int iess = 0;
for (int i = 0; i < tvec.Size(); i++)
{
if (iess < ess_tdofs_list.Size() && i == ess_tdofs_list[iess])
{
iess++;
continue;
}
CHECK(tvec[i] == Approx(-coef).epsilon(1e-8));
}
};

SECTION("Hex")
{
auto smesh = DividingPlaneMesh(false, true);
smesh.EnsureNCMesh();
ParMesh pmesh(MPI_COMM_WORLD, smesh);
test_project_H1(pmesh, 1, 0.25);
test_project_H1(pmesh, 2, 0.25);
test_project_H1(pmesh, 3, 0.25);
}

SECTION("Tet")
{
auto smesh = DividingPlaneMesh(true, true);
smesh.EnsureNCMesh();
ParMesh pmesh(MPI_COMM_WORLD, smesh);
test_project_H1(pmesh, 1, 0.25);
test_project_H1(pmesh, 2, 0.25);
test_project_H1(pmesh, 3, 0.25);
}
}

#endif // MFEM_USE_MPI

TEST_CASE("ReferenceCubeInternalBoundaries", "[NCMesh]")
Expand Down

0 comments on commit 2dabf82

Please sign in to comment.