Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fix projection onto NC internal faces #4259

Merged
merged 8 commits into from
May 28, 2024
Merged

Fix projection onto NC internal faces #4259

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
2dabf82
Fix bug where projection onto internal boundaries fails to identify p…
hughcars Apr 23, 2024
63804ab
Merge branch 'master' into hughcars/nc-internal-bdr-project-fix
hughcars May 8, 2024
8cc3111
Fix GetFaceDofs -> GetFaceVDofs
hughcars May 14, 2024
75bffa6
Address reviewer feedback:
hughcars May 15, 2024
71ad30f
Merge remote-tracking branch 'origin/master' into hughcars/nc-interna…
hughcars May 15, 2024
a723665
Fix double space, add assertion of coef > 0
hughcars May 20, 2024
c64f672
Merge remote-tracking branch 'origin/master' into hughcars/nc-interna…
hughcars May 20, 2024
770bcab
Merge branch 'master' into hughcars/nc-internal-bdr-project-fix
tzanio May 26, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
95 changes: 58 additions & 37 deletions fem/gridfunc.cpp
View file
Open in desktop
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,26 @@ void GridFunction::AccumulateAndCountBdrValues(
}
}
}
};

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

ElementTransformation *transf = mesh->GetEdgeTransformation(edge);
const FiniteElement *fe = fes->GetEdgeElement(edge);
mark_dofs(*transf, *fe);
}

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

ElementTransformation *transf = mesh->GetFaceTransformation(face);
const FiniteElement *fe = fes->GetFaceElement(face);
mark_dofs(*transf, *fe);
}
}
}
Expand Down Expand Up @@ -2228,26 +2238,37 @@ 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; }

T = mesh->GetEdgeTransformation(edge);
T->Attribute = -1; // TODO: set the boundary attribute
fe = fes->GetEdgeElement(edge);
lvec.SetSize(fe->GetDof());
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);
fe = fes->GetFaceElement(face);
lvec.SetSize(fe->GetDof());
fe->Project(vcoeff, *T, lvec);
accumulate_dofs(dofs, lvec, *this, values_counter);
}
}
}

Expand Down
76 changes: 76 additions & 0 deletions tests/unit/mesh/test_ncmesh.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2811,4 +2811,80 @@ TEST_CASE("RP=I", "[NCMesh]")
}
}


TEST_CASE("InternalBoundaryProjectBdrCoefficient", "[NCMesh]")
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should there be a parallel test? It seems the issue is fixed in parallel, considering the ex1p test you provided. Do you think the serial test is sufficient, since the changes are only in gridfunc.cpp?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The issue and fix are within GridFunction::AccumulateAndCountBdrValues that ParGridFunction::ProjectBdrCoefficient dispatches to, so fixing it for the serial case fixes it for the parallel. I did have the test for parallel, but it actually only failed for -np 1 which I think is possibly due to the extra steps in uncovering ghost boundary faces. I spent some time trying to devise a failing unit test for greater numbers of processors, but eventually concluded given the bug was rooted in the serial code, testing it directly there sufficed.

{
auto test_project_H1 = [](Mesh &mesh, int order, double coef)
dylan-copeland marked this conversation as resolved.
Show resolved Hide resolved
{
H1_FECollection fe_collection(order, mesh.SpaceDimension());
FiniteElementSpace fe_space(&mesh, &fe_collection);
GridFunction 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_vdofs_list, ess_vdofs_marker;
fe_space.GetEssentialVDofs(ess_bdr, ess_vdofs_marker);
fe_space.MarkerToList(ess_vdofs_marker, ess_vdofs_list);
for (auto ess_dof : ess_vdofs_list)
{
CHECK(x[ess_dof] == Approx(coef).epsilon(1e-8));
}

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

};

auto OneSidedNCRefine = [](Mesh &mesh)
{
// Pick one element attached to the new boundary attribute and refine.
const auto interface_attr = mesh.bdr_attributes.Max();
Array<int> el_to_ref;
for (int nbe = 0; nbe < mesh.GetNBE(); nbe++)
{
if (mesh.GetBdrAttribute(nbe) == interface_attr)
{
int f, o, e1, e2;
mesh.GetBdrElementFace(nbe, &f, &o);
mesh.GetFaceElements(f, &e1, &e2);
el_to_ref.Append(e1);
}
}
mesh.GeneralRefinement(el_to_ref);
return;
};

SECTION("Hex")
{
auto smesh = DividingPlaneMesh(false, true);
smesh.EnsureNCMesh(true);
OneSidedNCRefine(smesh);
test_project_H1(smesh, 2, 0.25);
}

SECTION("Tet")
{
auto smesh = DividingPlaneMesh(true, true);
smesh.EnsureNCMesh(true);
OneSidedNCRefine(smesh);
test_project_H1(smesh, 3, 0.25);
}
}



} // namespace mfem