Skip to content

Commit

Permalink
address review comments
Browse files Browse the repository at this point in the history
  • Loading branch information
pbrubeck committed Oct 16, 2024
1 parent a5c954d commit ff0aebe
Show file tree
Hide file tree
Showing 2 changed files with 49 additions and 47 deletions.
42 changes: 3 additions & 39 deletions tests/macro/test_macro_interp_project.py
Original file line number Diff line number Diff line change
Expand Up @@ -21,6 +21,7 @@ def h1_proj(u, f, bcs=None):
# compute h1 projection of f into u's function
# space, store the result in u.
v = TestFunction(u.function_space())

d = {H2: grad, H1: grad, HCurl: curl, HDiv: div, HDivDiv: div}[u.ufl_element().sobolev_space]
F = (inner(d(u-f), d(v)) * dx
+ inner(u-f, v) * dx)
Expand Down Expand Up @@ -123,53 +124,16 @@ def test_scalar_convergence(hierarchy, dim, el, deg, convrate, op):
[(2, 'Alfeld-Sorokina', 2, 2),
(3, 'Alfeld-Sorokina', 2, 2),
(2, 'Reduced-Arnold-Qin', 2, 1),
(3, 'Guzman-Neilan', 1, 1),
(3, 'Christiansen-Hu', 1, 1),
(2, 'Bernardi-Raugel', 1, 1),
(3, 'Bernardi-Raugel', 1, 1),
(2, 'Johnson-Mercier', 1, 1),
(3, 'Johnson-Mercier', 1, 1),
(3, 'Guzman-Neilan 1st kind H1', 1, 1),
(3, 'Guzman-Neilan H1(div)', 3, 2),
])
@pytest.mark.parametrize('op', (proj, h1_proj))
def test_piola_convergence(hierarchy, dim, el, deg, convrate, op):
if op == proj:
convrate += 1
run_convergence(hierarchy[dim], el, deg, convrate, op)


# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions
def test_supersmooth_bcs(mesh):
tdim = mesh.topological_dimension()
if tdim == 3:
V = FunctionSpace(mesh, "GNH1div", 3)
else:
V = FunctionSpace(mesh, "Alfeld-Sorokina", 2)

# check that V in H1
assert V.ufl_element().sobolev_space == H1

# check that V is supersmooth
nodes = V.finat_element.fiat_equivalent.dual.nodes
deriv_nodes = [i for i, node in enumerate(nodes) if len(node.deriv_dict)]
assert len(deriv_nodes) == tdim + 1

deriv_ids = V.cell_node_list[:, deriv_nodes]
u = Function(V)

CG = FunctionSpace(mesh, "Lagrange", 2)
RT = FunctionSpace(mesh, "RT", 1)
for sub in [1, (1, 2), "on_boundary"]:
bc = DirichletBC(V, 0, sub)

# check that we have the expected number of bc nodes
nnodes = len(bc.nodes)
expected = tdim * len(DirichletBC(CG, 0, sub).nodes)
if tdim == 3:
expected += len(DirichletBC(RT, 0, sub).nodes)
assert nnodes == expected

# check that the bc does not set the derivative nodes
u.assign(111)
u.dat.data_wo[deriv_ids] = 42
bc.zero(u)
assert numpy.allclose(u.dat.data_ro[deriv_ids], 42)
54 changes: 46 additions & 8 deletions tests/macro/test_stokes_macroelements.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,6 @@
from firedrake import *
import pytest
import numpy
from firedrake import *


@pytest.fixture(params=("square", "cube"))
Expand All @@ -11,21 +12,18 @@ def mesh(request):
return UnitCubeMesh(n, n, n)


@pytest.fixture(params=("SV", "GN", "GN2", "AS"))
@pytest.fixture(params=("SV", "GN", "GN2", "GNH1div"))
def space(request, mesh):
family = request.param
dim = mesh.topological_dimension()
if family == "GN":
if family == "GN1":
V = FunctionSpace(mesh, "GN", 1)
Q = FunctionSpace(mesh, "DG", 0)
elif family == "GN2":
V = FunctionSpace(mesh, "GN2", 1)
Q = FunctionSpace(mesh, "DG", 0, variant="alfeld")
elif family == "AS":
if dim == 2:
V = FunctionSpace(mesh, "AS", 2)
else:
V = FunctionSpace(mesh, "GNH1div", dim)
elif family == "GNH1div":
V = FunctionSpace(mesh, "GNH1div", dim)
Q = FunctionSpace(mesh, "CG", 1, variant="alfeld")
elif family == "SV":
V = VectorFunctionSpace(mesh, "CG", dim, variant="alfeld")
Expand All @@ -44,3 +42,43 @@ def test_stokes_complex(mesh, space):
p.interpolate(div(u))
assert norm(div(u) - p) < 1E-10
u.dat.data_wo[k] = 0


# Test that DirichletBC does not set derivative nodes of supersmooth H1 functions
def test_supersmooth_bcs(mesh):
tdim = mesh.topological_dimension()
if tdim == 3:
V = FunctionSpace(mesh, "GNH1div", 3)
else:
V = FunctionSpace(mesh, "Alfeld-Sorokina", 2)

assert V.finat_element.complex.is_macrocell()

# check that V in H1
assert V.ufl_element().sobolev_space == H1

# check that V is supersmooth
nodes = V.finat_element.fiat_equivalent.dual.nodes
deriv_nodes = [i for i, node in enumerate(nodes) if len(node.deriv_dict)]
assert len(deriv_nodes) == tdim + 1

deriv_ids = V.cell_node_list[:, deriv_nodes]
u = Function(V)

CG = FunctionSpace(mesh, "Lagrange", 2)
RT = FunctionSpace(mesh, "RT", 1)
for sub in [1, (1, 2), "on_boundary"]:
bc = DirichletBC(V, 0, sub)

# check that we have the expected number of bc nodes
nnodes = len(bc.nodes)
expected = tdim * len(DirichletBC(CG, 0, sub).nodes)
if tdim == 3:
expected += len(DirichletBC(RT, 0, sub).nodes)
assert nnodes == expected

# check that the bc does not set the derivative nodes
u.assign(111)
u.dat.data_wo[deriv_ids] = 42
bc.zero(u)
assert numpy.allclose(u.dat.data_ro[deriv_ids], 42)

0 comments on commit ff0aebe

Please sign in to comment.