Skip to content

Commit a5dbb11

Browse files
ksagiyamksagiyam
committed
codim 1
1 parent d1bcfe9 commit a5dbb11

File tree

8 files changed

+265
-2
lines changed

8 files changed

+265
-2
lines changed

‎firedrake/assemble.py‎

Lines changed: 33 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1646,6 +1646,7 @@ def __init__(self, form, local_knl, subdomain_id, all_integer_subdomain_ids, dia
16461646
self._constants = _FormHandler.iter_constants(form, local_knl.kinfo)
16471647
self._active_exterior_facets = _FormHandler.iter_active_exterior_facets(form, local_knl.kinfo)
16481648
self._active_interior_facets = _FormHandler.iter_active_interior_facets(form, local_knl.kinfo)
1649+
self._active_orientations_cell = _FormHandler.iter_active_orientations_cell(form, local_knl.kinfo)
16491650
self._active_orientations_exterior_facet = _FormHandler.iter_active_orientations_exterior_facet(form, local_knl.kinfo)
16501651
self._active_orientations_interior_facet = _FormHandler.iter_active_orientations_interior_facet(form, local_knl.kinfo)
16511652

@@ -1668,6 +1669,7 @@ def build(self):
16681669
assert_empty(self._constants)
16691670
assert_empty(self._active_exterior_facets)
16701671
assert_empty(self._active_interior_facets)
1672+
assert_empty(self._active_orientations_cell)
16711673
assert_empty(self._active_orientations_exterior_facet)
16721674
assert_empty(self._active_orientations_interior_facet)
16731675

@@ -1865,6 +1867,17 @@ def _as_global_kernel_arg_interior_facet(_, self):
18651867
return op2.DatKernelArg((2,), m._global_kernel_arg)
18661868

18671869

1870+
@_as_global_kernel_arg.register(kernel_args.OrientationsCellKernelArg)
1871+
def _(_, self):
1872+
mesh = next(self._active_orientations_cell)
1873+
if mesh is self._mesh:
1874+
return op2.DatKernelArg((1,))
1875+
else:
1876+
m, integral_type = mesh.topology.trans_mesh_entity_map(self._mesh.topology, self._integral_type, self._subdomain_id, self._all_integer_subdomain_ids)
1877+
assert integral_type == "cell"
1878+
return op2.DatKernelArg((1,), m._global_kernel_arg)
1879+
1880+
18681881
@_as_global_kernel_arg.register(kernel_args.OrientationsExteriorFacetKernelArg)
18691882
def _(_, self):
18701883
mesh = next(self._active_orientations_exterior_facet)
@@ -1936,6 +1949,7 @@ def __init__(self, form, bcs, local_knl, subdomain_id,
19361949
self._constants = _FormHandler.iter_constants(form, local_knl.kinfo)
19371950
self._active_exterior_facets = _FormHandler.iter_active_exterior_facets(form, local_knl.kinfo)
19381951
self._active_interior_facets = _FormHandler.iter_active_interior_facets(form, local_knl.kinfo)
1952+
self._active_orientations_cell = _FormHandler.iter_active_orientations_cell(form, local_knl.kinfo)
19391953
self._active_orientations_exterior_facet = _FormHandler.iter_active_orientations_exterior_facet(form, local_knl.kinfo)
19401954
self._active_orientations_interior_facet = _FormHandler.iter_active_orientations_interior_facet(form, local_knl.kinfo)
19411955

@@ -2201,6 +2215,17 @@ def _as_parloop_arg_interior_facet(_, self):
22012215
return op2.DatParloopArg(mesh.interior_facets.local_facet_dat, m)
22022216

22032217

2218+
@_as_parloop_arg.register(kernel_args.OrientationsCellKernelArg)
2219+
def _(_, self):
2220+
mesh = next(self._active_orientations_cell)
2221+
if mesh is self._mesh:
2222+
m = None
2223+
else:
2224+
m, integral_type = mesh.topology.trans_mesh_entity_map(self._mesh.topology, self._integral_type, self._subdomain_id, self._all_integer_subdomain_ids)
2225+
assert integral_type == "cell"
2226+
return op2.DatParloopArg(mesh.local_cell_orientation_dat, m)
2227+
2228+
22042229
@_as_parloop_arg.register(kernel_args.OrientationsExteriorFacetKernelArg)
22052230
def _(_, self):
22062231
mesh = next(self._active_orientations_exterior_facet)
@@ -2297,6 +2322,14 @@ def iter_active_interior_facets(form, kinfo):
22972322
mesh = all_meshes[i]
22982323
yield mesh
22992324

2325+
@staticmethod
2326+
def iter_active_orientations_cell(form, kinfo):
2327+
"""Yield the form cell orientations referenced in ``kinfo``."""
2328+
all_meshes = extract_domains(form)
2329+
for i in kinfo.active_domain_numbers.orientations_cell:
2330+
mesh = all_meshes[i]
2331+
yield mesh
2332+
23002333
@staticmethod
23012334
def iter_active_orientations_exterior_facet(form, kinfo):
23022335
"""Yield the form exterior facet orientations referenced in ``kinfo``."""

‎firedrake/mesh.py‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -709,6 +709,12 @@ def entity_orientations(self):
709709
"""
710710
pass
711711

712+
@property
713+
@abc.abstractmethod
714+
def local_cell_orientation_dat(self):
715+
"""Local cell orientation dat."""
716+
pass
717+
712718
@abc.abstractmethod
713719
def _facets(self, kind):
714720
pass
@@ -1287,6 +1293,16 @@ def cell_closure(self):
12871293
def entity_orientations(self):
12881294
return dmcommon.entity_orientations(self, self.cell_closure)
12891295

1296+
@utils.cached_property
1297+
def local_cell_orientation_dat(self):
1298+
"""Local cell orientation dat."""
1299+
return op2.Dat(
1300+
op2.DataSet(self.cell_set, 1),
1301+
self.entity_orientations[:, [-1]],
1302+
gem.uint_type,
1303+
f"{self.name}_local_cell_orientation"
1304+
)
1305+
12901306
@PETSc.Log.EventDecorator()
12911307
def _facets(self, kind):
12921308
if kind not in ["interior", "exterior"]:
@@ -1822,6 +1838,11 @@ def cell_closure(self):
18221838
def entity_orientations(self):
18231839
return self._base_mesh.entity_orientations
18241840

1841+
@utils.cached_property
1842+
def local_cell_orientation_dat(self):
1843+
"""Local cell orientation dat."""
1844+
return self._base_mesh.local_cell_orientation_dat
1845+
18251846
def _facets(self, kind):
18261847
if kind not in ["interior", "exterior"]:
18271848
raise ValueError("Unknown facet type '%s'" % kind)
@@ -2085,6 +2106,11 @@ def cell_closure(self):
20852106

20862107
entity_orientations = None
20872108

2109+
@property
2110+
def local_cell_orientation_dat(self):
2111+
"""Local cell orientation dat."""
2112+
raise NotImplementedError("Not implemented for VertexOnlyMeshTopology")
2113+
20882114
def _facets(self, kind):
20892115
"""Raises an AttributeError since cells in a
20902116
`VertexOnlyMeshTopology` have no facets.

‎firedrake/slate/slac/compiler.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -200,6 +200,7 @@ def generate_loopy_kernel(slate_expr, compiler_parameters=None):
200200
cell_sizes=(0, ) if builder.bag.needs_cell_sizes else (),
201201
exterior_facets=(),
202202
interior_facets=(),
203+
orientations_cell=(),
203204
orientations_exterior_facet=(),
204205
orientations_interior_facet=(),),
205206
coefficient_numbers=coefficient_numbers,

‎tests/firedrake/submesh/test_submesh_assemble.py‎

Lines changed: 99 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -326,3 +326,102 @@ def test_submesh_assemble_cell_cell_equation_bc():
326326
[- 1. / 3., - 1. / 6., 0., 0.]])
327327
assert np.allclose(A.M[0][0].values, M00)
328328
assert np.allclose(A.M[0][1].values, M01)
329+
330+
331+
def test_submesh_assemble_cell_facet_integral_various():
332+
# CG1 DoF numbers (nprocs = 1):
333+
#
334+
# 5-------1-------2
335+
# | | |
336+
# | | | mesh
337+
# | | |
338+
# 4-------0-------3
339+
#
340+
# 0
341+
# |
342+
# | subm
343+
# |
344+
# 1
345+
#
346+
distribution_parameters={
347+
"overlap_type": (DistributedMeshOverlapType.RIDGE, 1),
348+
}
349+
subdomain_id = 777
350+
mesh = RectangleMesh(2, 1, 2., 1., quadrilateral=True, distribution_parameters=distribution_parameters)
351+
x, y = SpatialCoordinate(mesh)
352+
V1 = FunctionSpace(mesh, "HDiv Trace", 0)
353+
f1 = Function(V1).interpolate(conditional(And(x > 0.9, x < 1.1), 1., 0.))
354+
mesh = RelabeledMesh(mesh, [f1], [subdomain_id])
355+
x, y = SpatialCoordinate(mesh)
356+
subm = Submesh(mesh, mesh.topological_dimension() - 1, subdomain_id)
357+
subx, suby = SpatialCoordinate(subm)
358+
V0 = FunctionSpace(mesh, "CG", 1)
359+
V1 = FunctionSpace(subm, "CG", 1)
360+
V = V0 * V1
361+
u = TrialFunction(V)
362+
v = TestFunction(V)
363+
u0, u1 = split(u)
364+
v0, v1 = split(v)
365+
coordV0 = VectorFunctionSpace(mesh, "CG", 1)
366+
coordV1 = VectorFunctionSpace(subm, "CG", 1)
367+
coordV = coordV0 * coordV1
368+
coords = Function(coordV)
369+
coords.sub(0).assign(mesh.coordinates)
370+
coords.sub(1).assign(subm.coordinates)
371+
coords0, coords1 = split(coords)
372+
M10 = np.array(
373+
[
374+
[1. / 6., 1. / 3., 0., 0., 0., 0.],
375+
[1. / 3., 1. / 6., 0., 0., 0., 0.],
376+
]
377+
)
378+
M10w = np.array(
379+
[
380+
[1. / 12., 1. / 4., 0., 0., 0., 0.],
381+
[1. / 12., 1. / 12., 0., 0., 0., 0.],
382+
]
383+
)
384+
M10ww = np.array(
385+
[
386+
[1. / 20., 1. / 5., 0., 0., 0., 0.],
387+
[1. / 30., 1. / 20., 0., 0., 0., 0.],
388+
]
389+
)
390+
# Use subm as primal integration domain.
391+
measure = Measure(
392+
"dx", subm,
393+
intersect_measures=(
394+
Measure("dS", mesh),
395+
),
396+
)
397+
a = inner(u0('-'), v1) * measure
398+
A = assemble(a, mat_type="nest")
399+
assert np.allclose(A.M[1][0].values, M10)
400+
a = inner(u1, v0('+')) * measure
401+
A = assemble(a, mat_type="nest")
402+
assert np.allclose(A.M[0][1].values, np.transpose(M10))
403+
a = y * inner(u0('-'), v1) * measure
404+
A = assemble(a, mat_type="nest")
405+
assert np.allclose(A.M[1][0].values, M10w)
406+
a = y * suby * inner(u0('-'), v1) * measure
407+
A = assemble(a, mat_type="nest")
408+
assert np.allclose(A.M[1][0].values, M10ww)
409+
a = coords0[1] * inner(u0('-'), v1) * measure
410+
A = assemble(a, mat_type="nest")
411+
assert np.allclose(A.M[1][0].values, M10w)
412+
a = coords0[1] * coords1[1] * inner(u0('-'), v1) * measure
413+
A = assemble(a, mat_type="nest")
414+
assert np.allclose(A.M[1][0].values, M10ww)
415+
# Use mesh as primal integration domain.
416+
measure = Measure(
417+
"dS", mesh,
418+
intersect_measures=(
419+
Measure("dx", subm),
420+
),
421+
)
422+
a = inner(u0('+'), v1) * measure(subdomain_id)
423+
A = assemble(a, mat_type="nest")
424+
assert np.allclose(A.M[1][0].values, M10)
425+
a = inner(u1, v0('-')) * measure(subdomain_id)
426+
A = assemble(a, mat_type="nest")
427+
assert np.allclose(A.M[0][1].values, np.transpose(M10))

‎tests/firedrake/submesh/test_submesh_solve.py‎

Lines changed: 83 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -458,3 +458,86 @@ def test_submesh_solve_cell_cell_equation_bc(nref, degree, simplex):
458458
solve(a == L, sol, bcs=[dbc, ebc])
459459
assert sqrt(assemble(inner(sol[0] - x * y, sol[0] - x * y) * dx_outer)) < 1.e-12
460460
assert sqrt(assemble(inner(sol[1] - x * y, sol[1] - x * y) * dx_inner)) < 1.e-12
461+
462+
463+
@pytest.mark.parallel(nprocs=8)
464+
@pytest.mark.parametrize('simplex', [True, False])
465+
@pytest.mark.parametrize('direction', [0, 1, 2])
466+
def test_submesh_solve_cell_facet_3d(simplex, direction):
467+
dim = 3
468+
distribution_parameters={
469+
"overlap_type": (DistributedMeshOverlapType.RIDGE, 1),
470+
}
471+
if simplex:
472+
nref = 3
473+
mesh = BoxMesh(2 ** nref, 2 ** nref, 2 ** nref, 1., 1., 1., hexahedral=False, distribution_parameters=distribution_parameters)
474+
xyz = SpatialCoordinate(mesh)
475+
DG0 = FunctionSpace(mesh, "DG", 0)
476+
c1 = Function(DG0).interpolate(conditional(xyz[direction] < 0.318310, 1, 0))
477+
c2 = Function(DG0).interpolate(conditional(xyz[direction] > 0.318310, 1, 0))
478+
mesh = RelabeledMesh(mesh, [c1, c2], [1, 2])
479+
family = "P"
480+
else:
481+
mesh = Mesh(join(cwd, "..", "meshes", "cube_hex.msh"), distribution_parameters=distribution_parameters)
482+
xyz = SpatialCoordinate(mesh)
483+
DG0 = FunctionSpace(mesh, "DQ", 0)
484+
c1 = Function(DG0).interpolate(conditional(xyz[direction] < 0.318310, 1, 0))
485+
c2 = Function(DG0).interpolate(conditional(xyz[direction] > 0.318310, 1, 0))
486+
HDivTrace0 = FunctionSpace(mesh, "Q", 2)
487+
f1 = Function(HDivTrace0).interpolate(conditional(xyz[0] < .001, 1, 0))
488+
f2 = Function(HDivTrace0).interpolate(conditional(xyz[0] > .999, 1, 0))
489+
f3 = Function(HDivTrace0).interpolate(conditional(xyz[1] < .001, 1, 0))
490+
f4 = Function(HDivTrace0).interpolate(conditional(xyz[1] < .999, 1, 0))
491+
f5 = Function(HDivTrace0).interpolate(conditional(xyz[2] < .001, 1, 0))
492+
f6 = Function(HDivTrace0).interpolate(conditional(xyz[2] < .999, 1, 0))
493+
mesh = RelabeledMesh(mesh, [c1, c2, f1, f2, f3, f4, f5, f6], [1, 2, 1, 2, 3, 4, 5, 6])
494+
family = "Q"
495+
mesh1 = Submesh(mesh, dim, 1)
496+
x1, y1, z1 = SpatialCoordinate(mesh1)
497+
mesh2 = Submesh(mesh, dim, 2)
498+
x2, y2, z2 = SpatialCoordinate(mesh2)
499+
label_interf = 7 # max + 1
500+
mesh12 = Submesh(mesh2, dim - 1, label_interf)
501+
dx1 = Measure("dx", mesh1)
502+
dx2 = Measure("dx", mesh2)
503+
dx12 = Measure("dx", mesh12)
504+
dx12_ds1 = Measure("dx", mesh12, intersect_measures=(Measure("ds", mesh1),))
505+
ds1_dx12 = Measure("ds", mesh1, intersect_measures=(Measure("dx", mesh12),))
506+
ds2_dx12 = Measure("ds", mesh2, intersect_measures=(Measure("dx", mesh12),))
507+
# Check sanity.
508+
vol1 = assemble(Constant(1) * dx1)
509+
vol2 = assemble(Constant(1) * dx2)
510+
assert abs(vol1 + vol2 - 1.) < 1.e-14
511+
# Solve Poisson problem.
512+
V1 = FunctionSpace(mesh1, family, 3)
513+
V12 = FunctionSpace(mesh12, family, 5)
514+
V2 = FunctionSpace(mesh2, family, 4)
515+
V = V1 * V12 * V2
516+
u = TrialFunction(V)
517+
v = TestFunction(V)
518+
u1, u12, u2 = split(u)
519+
v1, v12, v2 = split(v)
520+
g1 = cos(2 * pi * x1) * cos(2 * pi * y1) * cos(2 * pi * z1)
521+
g2 = cos(2 * pi * x2) * cos(2 * pi * y2) * cos(2 * pi * z2)
522+
f1 = 12 * pi**2 * g1
523+
f2 = 12 * pi**2 * g2
524+
n1 = FacetNormal(mesh1)
525+
a = (
526+
inner(grad(u1), grad(v1)) * dx1 - inner(u12, v1) * ds1_dx12(label_interf) +
527+
inner(u12, v12) * dx12 +
528+
inner(grad(u2), grad(v2)) * dx2 + inner(u12, v2) * ds2_dx12(label_interf)
529+
)
530+
L = (
531+
inner(f1, v1) * dx1 +
532+
inner(dot(grad(g1), n1), v12) * dx12_ds1 +
533+
inner(f2, v2) * dx2
534+
)
535+
sol = Function(V)
536+
bc1 = DirichletBC(V.sub(0), g1, (2 * direction + 1,))
537+
bc2 = DirichletBC(V.sub(2), g2, (2 * direction + 2,))
538+
solve(a == L, sol, bcs=[bc1, bc2])
539+
sol1, sol12, sol2 = split(sol)
540+
error1 = assemble(inner(sol1 - g1, sol1 - g1) * dx1)
541+
error2 = assemble(inner(sol2 - g2, sol2 - g2) * dx2)
542+
assert sqrt(error1) < 4.e-4
543+
assert sqrt(error2) < 4.e-5

‎tsfc/fem.py‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
import numpy
1010
import ufl
1111
from FIAT.orientation_utils import Orientation as FIATOrientation
12-
from FIAT.reference_element import UFCHexahedron, UFCSimplex, make_affine_mapping
12+
from FIAT.reference_element import UFCHexahedron, UFCQuadrilateral, UFCSimplex, make_affine_mapping
1313
from FIAT.reference_element import TensorProductCell
1414
from finat.physically_mapped import (NeedsCoordinateMappingElement,
1515
PhysicalGeometry)
@@ -130,6 +130,9 @@ def use_canonical_quadrature_point_ordering(self):
130130
return True
131131
elif len(set(cell_integral_type_map)) > 1: # mixed cell types
132132
return True
133+
elif any(isinstance(cell, UFCHexahedron) for cell in cell_integral_type_map) and \
134+
any(isinstance(cell, UFCQuadrilateral) for cell in cell_integral_type_map):
135+
return True
133136
return False
134137

135138

‎tsfc/kernel_args.py‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -54,6 +54,10 @@ class InteriorFacetKernelArg(KernelArg):
5454
...
5555

5656

57+
class OrientationsCellKernelArg(KernelArg):
58+
...
59+
60+
5761
class OrientationsExteriorFacetKernelArg(KernelArg):
5862
...
5963

0 commit comments

Comments
 (0)