-
Notifications
You must be signed in to change notification settings - Fork 169
Description
Describe the bug
It appears that we are failing to get the facet orientations correctly on non-Lagrange bases for Q
. But this only happens in parallel if the mesh distribution is not overlapping.
Steps to Reproduce
Steps to reproduce the behavior:
Here's an MFE assembling a Form on a quadrilateral mesh with two elements using Lagrange (variant="spectral") and non-Lagrange (variant="integral") with the dfault overlap (FACET
) and without overlap (NONE
).
from firedrake import *
from firedrake.petsc import PETSc
degree = 3
overlaps = {"NONE": (DistributedMeshOverlapType.NONE, 0),
"FACET": (DistributedMeshOverlapType.FACET, 1),}
variants = ("spectral", "integral")
for variant in variants:
for overlap in overlaps:
distribution = {"overlap_type": overlaps[overlap]}
mesh = UnitSquareMesh(2, 1, quadrilateral=True, distribution_parameters=distribution)
V = FunctionSpace(mesh, "Lagrange", degree, variant=variant)
c = assemble(TestFunction(V) * dx)
# Zero out interior dofs
own = mesh.cell_set.size
interior = V.finat_element.entity_dofs()[2][0]
c.dat.data[V.cell_node_list[:own, interior]] = 0
# Zero out boundary dofs
bc = DirichletBC(V, 0, "on_boundary")
bc.zero(c)
c.dat.data[abs(c.dat.data) < 1E-10] = 0
with c.dat.vec as vec:
PETSc.Sys.Print(variant, distribution, vec.norm())
vec.view()
Output of mpiexec -n 2 python script.py
:
spectral {'overlap_type': (<DistributedMeshOverlapType.NONE: 1>, 0)} 0.04910463758239919
spectral {'overlap_type': (<DistributedMeshOverlapType.FACET: 2>, 1)} 0.04910463758239919
integral {'overlap_type': (<DistributedMeshOverlapType.NONE: 1>, 0)} 0.12500000000000014
integral {'overlap_type': (<DistributedMeshOverlapType.FACET: 2>, 1)} 0.1767766952966371
Expected behavior
The norm of the assembled Cofunction should be independent of the mesh distribution, but it should still depend on the variant.
Error message
No error message.
Additional Info
The MFE zeros out all entries except for the interface dofs on the interior of the domain. For variant="integral"
, the overlapping case gives the correct solution [0., -0.176777], whereas it appears that the non-overlapping case gives [-0.088388 -0.0883883]. This suggests that we are failing to map the contributions coming from each cell with the right entity permutation.