Skip to content

Commit

Permalink
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
introduce MixedMesh
Browse files Browse the repository at this point in the history
ksagiyam committed Jul 25, 2024

Verified

This commit was created on GitHub.com and signed with GitHub’s verified signature. The key has expired.
1 parent 8db60e0 commit 9be8132
Showing 11 changed files with 574 additions and 242 deletions.
17 changes: 16 additions & 1 deletion gem/gem.py
Original file line number Diff line number Diff line change
@@ -33,7 +33,7 @@
'IndexSum', 'ListTensor', 'Concatenate', 'Delta',
'index_sum', 'partial_indexed', 'reshape', 'view',
'indices', 'as_gem', 'FlexiblyIndexed',
'Inverse', 'Solve', 'extract_type']
'Inverse', 'Solve', 'extract_type', 'extract_variables']


class NodeMeta(type):
@@ -1044,3 +1044,18 @@ def as_gem(expr):
def extract_type(expressions, klass):
"""Collects objects of type klass in expressions."""
return tuple(node for node in traversal(expressions) if isinstance(node, klass))


def extract_variables(expressions):
"""Collects variables in expressions."""
variables = set()
for node in traversal(expressions):
if isinstance(node, Variable):
variables.add(node)
elif isinstance(node, Indexed):
# Need to go deeper to find variables wrapped in 'VariableIndex's.
for idx in node.multiindex:
if isinstance(idx, VariableIndex):
v, = extract_variables((idx.expression, ))
variables.add(v)
return variables
224 changes: 224 additions & 0 deletions tests/test_mixed_function_space_with_mixed_mesh.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,224 @@
from collections import defaultdict
from tsfc import compile_form
from ufl import (triangle, Mesh, MixedMesh, FunctionSpace, TestFunction, TrialFunction, Coefficient,
Measure, SpatialCoordinate, inner, grad, curl, div, split, as_vector, )
from ufl.pullback import identity_pullback, contravariant_piola
from ufl.sobolevspace import H1, HDiv, L2
from finat.ufl import FiniteElement, MixedElement, VectorElement
from tsfc.ufl_utils import compute_form_data
from tsfc import kernel_args


def test_mixed_function_space_with_mixed_mesh_restrictions_base():
cell = triangle
gdim = 2
elem0 = FiniteElement("Discontinuous Lagrange", cell, 2)
elem1 = FiniteElement("Discontinuous Lagrange", cell, 3)
elem = MixedElement([elem0, elem1])
mesh0 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=100)
mesh1 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=101)
domain = MixedMesh(mesh0, mesh1)
V = FunctionSpace(domain, elem)
V0 = FunctionSpace(mesh0, elem0)
V1 = FunctionSpace(mesh1, elem1)
f = Coefficient(V, count=1000)
f0, f1 = split(f)
u0 = TrialFunction(V0)
u1 = TrialFunction(V1)
v0 = TestFunction(V0)
v1 = TestFunction(V1)
dx0 = Measure("dx", mesh0)
dx1 = Measure("dx", mesh1)
ds0 = Measure("ds", mesh0)
ds1 = Measure("ds", mesh1)
dS0 = Measure("dS", mesh0)
dS1 = Measure("dS", mesh1)
f0_split = Coefficient(V0)
f1_split = Coefficient(V1)
# a
form = inner(grad(f1('|')), as_vector([1, 0])) * ds1(777)
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 1
assert integral_data.domain_integral_type_map[mesh1] == "exterior_facet"
# b
form = inner(grad(f1('|')), grad(f1('|'))) * dS0(777)
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 2
assert integral_data.domain_integral_type_map[mesh0] == "interior_facet"
assert integral_data.domain_integral_type_map[mesh1] == "exterior_facet"
# c
form = div(f) * inner(grad(f1), grad(f1)) * inner(grad(u1), grad(v0)) * dx1
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 2
assert integral_data.domain_integral_type_map[mesh0] == "cell"
assert integral_data.domain_integral_type_map[mesh1] == "cell"


def test_mixed_function_space_with_mixed_mesh_3_cg3_bdm3_dg2_dx1():
cell = triangle
gdim = 2
elem0 = FiniteElement("Lagrange", cell, 3)
elem1 = FiniteElement("Brezzi-Douglas-Marini", cell, 3)
elem2 = FiniteElement("Discontinuous Lagrange", cell, 2)
elem = MixedElement([elem0, elem1, elem2])
mesh0 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=100)
mesh1 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=101)
mesh2 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=102)
domain = MixedMesh(mesh0, mesh1, mesh2)
V = FunctionSpace(domain, elem)
V0 = FunctionSpace(mesh0, elem0)
V1 = FunctionSpace(mesh1, elem1)
V2 = FunctionSpace(mesh2, elem2)
f = Coefficient(V, count=1000)
u0 = TrialFunction(V0)
v1 = TestFunction(V1)
f0, f1, f2 = split(f)
f0_split = Coefficient(V0)
f1_split = Coefficient(V1)
f2_split = Coefficient(V2)
x2 = SpatialCoordinate(mesh2)
dx1 = Measure("dx", mesh1)
form = inner(x2, x2) * f2 * inner(grad(u0), v1) * dx1(999)
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split, f2_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 3
assert integral_data.domain_integral_type_map[mesh0] == "cell"
assert integral_data.domain_integral_type_map[mesh1] == "cell"
assert integral_data.domain_integral_type_map[mesh2] == "cell"
kernel, = compile_form(form)
assert kernel.domain_number == 0
assert kernel.integral_type == "cell"
assert kernel.subdomain_id == (999, )
assert kernel.active_domain_numbers.coordinates == (0, 1, 2)
assert kernel.active_domain_numbers.cell_orientations == ()
assert kernel.active_domain_numbers.cell_sizes == ()
assert kernel.active_domain_numbers.exterior_facets == ()
assert kernel.active_domain_numbers.interior_facets == ()
assert kernel.coefficient_numbers == ((0, (2, )), )
assert isinstance(kernel.arguments[0], kernel_args.OutputKernelArg)
assert isinstance(kernel.arguments[1], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[2], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[3], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[4], kernel_args.CoefficientKernelArg)
assert kernel.arguments[0].loopy_arg.shape == (20, 10)
assert kernel.arguments[1].loopy_arg.shape == (3 * gdim, )
assert kernel.arguments[2].loopy_arg.shape == (3 * gdim, )
assert kernel.arguments[3].loopy_arg.shape == (3 * gdim, )
assert kernel.arguments[4].loopy_arg.shape == (6, )


def test_mixed_function_space_with_mixed_mesh_restrictions_bdm3_dg2_dS0():
cell = triangle
gdim = 2
elem0 = FiniteElement("Brezzi-Douglas-Marini", cell, 3)
elem1 = FiniteElement("Discontinuous Lagrange", cell, 2)
elem = MixedElement([elem0, elem1])
mesh0 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=100)
mesh1 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=101)
domain = MixedMesh(mesh0, mesh1)
V = FunctionSpace(domain, elem)
V0 = FunctionSpace(mesh0, elem0)
V1 = FunctionSpace(mesh1, elem1)
f = Coefficient(V, count=1000)
f0, f1 = split(f)
f0_split = Coefficient(V0)
f1_split = Coefficient(V1)
u0 = TrialFunction(V0)
u1 = TrialFunction(V1)
v0 = TestFunction(V0)
v1 = TestFunction(V1)
dx0 = Measure("dx", mesh0)
dx1 = Measure("dx", mesh1)
ds0 = Measure("ds", mesh0)
ds1 = Measure("ds", mesh1)
dS0 = Measure("dS", mesh0)
dS1 = Measure("dS", mesh1)
form = inner(curl(f1('|')), curl(f1('|'))) * inner(grad(u1('|')), v0('+')) * dS0(777)
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 2
assert integral_data.domain_integral_type_map[mesh0] == "interior_facet"
assert integral_data.domain_integral_type_map[mesh1] == "exterior_facet"
kernel, = compile_form(form)
assert kernel.domain_number == 0
assert kernel.integral_type == "interior_facet"
assert kernel.subdomain_id == (777, )
assert kernel.active_domain_numbers.coordinates == (0, 1)
assert kernel.active_domain_numbers.cell_orientations == ()
assert kernel.active_domain_numbers.cell_sizes == ()
assert kernel.active_domain_numbers.exterior_facets == (1, )
assert kernel.active_domain_numbers.interior_facets == (0, )
assert kernel.coefficient_numbers == ((0, (1, )), )
assert isinstance(kernel.arguments[0], kernel_args.OutputKernelArg)
assert isinstance(kernel.arguments[1], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[2], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[3], kernel_args.CoefficientKernelArg)
assert isinstance(kernel.arguments[4], kernel_args.ExteriorFacetKernelArg)
assert isinstance(kernel.arguments[5], kernel_args.InteriorFacetKernelArg)
assert kernel.arguments[0].loopy_arg.shape == (2 * 20, 6)
assert kernel.arguments[1].loopy_arg.shape == (2 * (3 * gdim), )
assert kernel.arguments[2].loopy_arg.shape == (3 * gdim, )
assert kernel.arguments[3].loopy_arg.shape == (6, )
assert kernel.arguments[4].loopy_arg.shape == (1, )
assert kernel.arguments[5].loopy_arg.shape == (2, )


def test_mixed_function_space_with_mixed_mesh_restrictions_dg2_dg3_ds1():
cell = triangle
gdim = 2
elem0 = FiniteElement("Discontinuous Lagrange", cell, 2)
elem1 = FiniteElement("Discontinuous Lagrange", cell, 3)
elem = MixedElement([elem0, elem1])
mesh0 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=100)
mesh1 = Mesh(VectorElement("Lagrange", cell, 1), ufl_id=101)
domain = MixedMesh(mesh0, mesh1)
V = FunctionSpace(domain, elem)
V0 = FunctionSpace(mesh0, elem0)
V1 = FunctionSpace(mesh1, elem1)
f = Coefficient(V, count=1000)
f0_split = Coefficient(V0)
f1_split = Coefficient(V1)
f0, f1 = split(f)
u0 = TrialFunction(V0)
u1 = TrialFunction(V1)
v0 = TestFunction(V0)
v1 = TestFunction(V1)
dx0 = Measure("dx", mesh0)
dx1 = Measure("dx", mesh1)
ds0 = Measure("ds", mesh0)
ds1 = Measure("ds", mesh1)
dS0 = Measure("dS", mesh0)
dS1 = Measure("dS", mesh1)
form = inner(grad(f1('|')), grad(f0('-'))) * inner(grad(u0('-')), grad(v1('|'))) * ds1(777)
form_data = compute_form_data(form, do_split_coefficients={f: [f0_split, f1_split]})
integral_data, = form_data.integral_data
assert len(integral_data.domain_integral_type_map) == 2
assert integral_data.domain_integral_type_map[mesh0] == "interior_facet"
assert integral_data.domain_integral_type_map[mesh1] == "exterior_facet"
kernel, = compile_form(form)
assert kernel.domain_number == 0
assert kernel.integral_type == "exterior_facet"
assert kernel.subdomain_id == (777, )
assert kernel.active_domain_numbers.coordinates == (0, 1)
assert kernel.active_domain_numbers.cell_orientations == ()
assert kernel.active_domain_numbers.cell_sizes == ()
assert kernel.active_domain_numbers.exterior_facets == (0, )
assert kernel.active_domain_numbers.interior_facets == (1, )
assert kernel.coefficient_numbers == ((0, (0, 1)), )
assert isinstance(kernel.arguments[0], kernel_args.OutputKernelArg)
assert isinstance(kernel.arguments[1], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[2], kernel_args.CoordinatesKernelArg)
assert isinstance(kernel.arguments[3], kernel_args.CoefficientKernelArg)
assert isinstance(kernel.arguments[4], kernel_args.CoefficientKernelArg)
assert isinstance(kernel.arguments[5], kernel_args.ExteriorFacetKernelArg)
assert isinstance(kernel.arguments[6], kernel_args.InteriorFacetKernelArg)
assert kernel.arguments[0].loopy_arg.shape == (10, 2 * 6)
assert kernel.arguments[1].loopy_arg.shape == (1 * (3 * gdim), )
assert kernel.arguments[2].loopy_arg.shape == (2 * (3 * gdim), )
assert kernel.arguments[3].loopy_arg.shape == (2 * 6, )
assert kernel.arguments[4].loopy_arg.shape == (10, )
assert kernel.arguments[5].loopy_arg.shape == (1, )
assert kernel.arguments[6].loopy_arg.shape == (2, )
5 changes: 3 additions & 2 deletions tests/test_tsfc_182.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
import pytest

from ufl import Coefficient, TestFunction, dx, inner, tetrahedron, Mesh, FunctionSpace
from ufl import Coefficient, TestFunction, dx, inner, tetrahedron, Mesh, MixedMesh, FunctionSpace
from finat.ufl import FiniteElement, MixedElement, VectorElement

from tsfc import compile_form
@@ -20,7 +20,8 @@ def test_delta_elimination(mode):

element_chi_lambda = MixedElement(element_eps_p, element_lambda)
domain = Mesh(VectorElement("Lagrange", tetrahedron, 1))
space = FunctionSpace(domain, element_chi_lambda)
domains = MixedMesh(domain, domain)
space = FunctionSpace(domains, element_chi_lambda)

chi_lambda = Coefficient(space)
delta_chi_lambda = TestFunction(space)
5 changes: 3 additions & 2 deletions tests/test_tsfc_204.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,13 @@
from tsfc import compile_form
from ufl import (Coefficient, FacetNormal,
FunctionSpace, Mesh, as_matrix,
FunctionSpace, Mesh, MixedMesh, as_matrix,
dot, dS, ds, dx, facet, grad, inner, outer, split, triangle)
from finat.ufl import BrokenElement, FiniteElement, MixedElement, VectorElement


def test_physically_mapped_facet():
mesh = Mesh(VectorElement("P", triangle, 1))
meshes = MixedMesh(mesh, mesh, mesh, mesh, mesh)

# set up variational problem
U = FiniteElement("Morley", mesh.ufl_cell(), 2)
@@ -15,7 +16,7 @@ def test_physically_mapped_facet():
Vv = VectorElement(BrokenElement(V))
Qhat = VectorElement(BrokenElement(V[facet]))
Vhat = VectorElement(V[facet])
Z = FunctionSpace(mesh, MixedElement(U, Vv, Qhat, Vhat, R))
Z = FunctionSpace(meshes, MixedElement(U, Vv, Qhat, Vhat, R))

z = Coefficient(Z)
u, d, qhat, dhat, lam = split(z)
79 changes: 48 additions & 31 deletions tsfc/driver.py
Original file line number Diff line number Diff line change
@@ -8,7 +8,7 @@
from ufl.algorithms import extract_arguments, extract_coefficients
from ufl.algorithms.analysis import has_type
from ufl.classes import Form, GeometricQuantity
from ufl.domain import extract_unique_domain
from ufl.domain import extract_unique_domain, extract_domains, sort_domains

import gem
import gem.impero_utils as impero_utils
@@ -26,9 +26,9 @@


TSFCIntegralDataInfo = collections.namedtuple("TSFCIntegralDataInfo",
["domain", "integral_type", "subdomain_id", "domain_number",
["domain", "integral_type", "subdomain_id", "domain_number", "domain_integral_type_map",
"arguments",
"coefficients", "coefficient_numbers"])
"coefficients", "coefficient_split", "coefficient_numbers"])
TSFCIntegralDataInfo.__doc__ = """
Minimal set of objects for kernel builders.
@@ -47,7 +47,7 @@
"""


def compile_form(form, prefix="form", parameters=None, interface=None, diagonal=False, log=False, **kwargs):
def compile_form(form, prefix="form", parameters=None, dont_split_numbers=(), diagonal=False, log=False, **kwargs):
"""Compiles a UFL form into a set of assembly kernels.
:arg form: UFL form
@@ -77,29 +77,31 @@ def compile_form(form, prefix="form", parameters=None, interface=None, diagonal=

# Determine whether in complex mode:
complex_mode = parameters and is_complex(parameters.get("scalar_type"))
fd = ufl_utils.compute_form_data(form, complex_mode=complex_mode)
form_data = ufl_utils.compute_form_data(form,
do_split_coefficients=tuple(c for i, c in enumerate(form.coefficients()) if i not in dont_split_numbers),
complex_mode=complex_mode)
logger.info(GREEN % "compute_form_data finished in %g seconds.", time.time() - cpu_time)

kernels = []
for integral_data in fd.integral_data:
for integral_data in form_data.integral_data:
start = time.time()
kernel = compile_integral(integral_data, fd, prefix, parameters, interface=interface, diagonal=diagonal, log=log)
if kernel is not None:
kernels.append(kernel)
if integral_data.integrals:
kernel = compile_integral(integral_data, form_data, prefix, parameters, diagonal=diagonal, log=log)
if kernel is not None:
kernels.append(kernel)
logger.info(GREEN % "compile_integral finished in %g seconds.", time.time() - start)

logger.info(GREEN % "TSFC finished in %g seconds.", time.time() - cpu_time)
return kernels


def compile_integral(integral_data, form_data, prefix, parameters, interface, *, diagonal=False, log=False):
def compile_integral(integral_data, form_data, prefix, parameters, *, diagonal=False, log=False):
"""Compiles a UFL integral into an assembly kernel.
:arg integral_data: UFL integral data
:arg form_data: UFL form data
:arg prefix: kernel name will start with this string
:arg parameters: parameters object
:arg interface: backend module for the kernel interface
:arg diagonal: Are we building a kernel for the diagonal of a rank-2 element tensor?
:arg log: bool if the Kernel should be profiled with Log events
:returns: a kernel constructed by the kernel interface
@@ -108,48 +110,54 @@ def compile_integral(integral_data, form_data, prefix, parameters, interface, *,
integral_data.integral_type == "interior_facet":
raise NotImplementedError("interior facet integration in hex meshes not currently supported")
parameters = preprocess_parameters(parameters)
if interface is None:
interface = firedrake_interface_loopy.KernelBuilder
scalar_type = parameters["scalar_type"]
integral_type = integral_data.integral_type
if integral_type.startswith("interior_facet") and diagonal:
raise NotImplementedError("Sorry, we can't assemble the diagonal of a form for interior facet integrals")
mesh = integral_data.domain
arguments = form_data.preprocessed_form.arguments()
kernel_name = f"{prefix}_{integral_type}_integral"
# Dict mapping domains to index in original_form.ufl_domains()
domain_numbering = form_data.original_form.domain_numbering()
domain_number = domain_numbering[integral_data.domain]
coefficients = [form_data.function_replace_map[c] for c in integral_data.integral_coefficients]
# This is which coefficient in the original form the
# current coefficient is.
# Consider f*v*dx + g*v*ds, the full form contains two
# coefficients, but each integral only requires one.
coefficient_numbers = tuple(form_data.original_coefficient_positions[i]
for i, (_, enabled) in enumerate(zip(form_data.reduced_coefficients, integral_data.enabled_coefficients))
if enabled)
coefficients = []
coefficient_split = {}
coefficient_numbers = []
for i, (coeff_orig, enabled) in enumerate(zip(form_data.reduced_coefficients, integral_data.enabled_coefficients)):
if enabled:
coeff = form_data.function_replace_map[coeff_orig]
coefficients.append(coeff)
if coeff in form_data.coefficient_split:
coefficient_split[coeff] = form_data.coefficient_split[coeff]
coefficient_numbers.append(form_data.original_coefficient_positions[i])
mesh = integral_data.domain
all_meshes = extract_domains(form_data.original_form)
domain_number = all_meshes.index(mesh)
integral_data_info = TSFCIntegralDataInfo(domain=integral_data.domain,
integral_type=integral_data.integral_type,
subdomain_id=integral_data.subdomain_id,
domain_number=domain_number,
domain_integral_type_map={mesh: integral_data.domain_integral_type_map[mesh] if mesh in integral_data.domain_integral_type_map else None for mesh in all_meshes},
arguments=arguments,
coefficients=coefficients,
coefficient_split=coefficient_split,
coefficient_numbers=coefficient_numbers)
builder = interface(integral_data_info,
scalar_type,
diagonal=diagonal)
builder.set_coordinates(mesh)
builder.set_cell_sizes(mesh)
builder.set_coefficients(integral_data, form_data)
builder = firedrake_interface_loopy.KernelBuilder(integral_data_info,
scalar_type,
diagonal=diagonal)
builder.set_entity_numbers(all_meshes)
builder.set_coordinates(all_meshes)
builder.set_cell_orientations(all_meshes)
builder.set_cell_sizes(all_meshes)
builder.set_coefficients()
# TODO: We do not want pass constants to kernels that do not need them
# so we should attach the constants to integral data instead
builder.set_constants(form_data.constants)
ctx = builder.create_context()
for integral in integral_data.integrals:
params = parameters.copy()
params.update(integral.metadata()) # integral metadata overrides
integrand = ufl.replace(integral.integrand(), form_data.function_replace_map)
integrand_exprs = builder.compile_integrand(integrand, params, ctx)
integrand_exprs = builder.compile_integrand(integral.integrand(), params, ctx)
integral_exprs = builder.construct_integrals(integrand_exprs, params)
builder.stash_integrals(integral_exprs, params, ctx)
return builder.construct_kernel(kernel_name, ctx, log)
@@ -170,6 +178,13 @@ def preprocess_parameters(parameters):
return parameters


def collect_domains_in_integrals(mesh, integrals):
meshes = set((mesh, ))
for integral in integrals:
meshes.update(extract_domains(integral.integrand()))
return sort_domains(meshes)


def compile_expression_dual_evaluation(expression, to_element, ufl_element, *,
domain=None, interface=None,
parameters=None, log=False):
@@ -222,6 +237,7 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *,
if domain is None:
domain = extract_unique_domain(expression)
assert domain is not None
builder._domain_integral_type_map = {domain: "cell"}

# Collect required coefficients and determine numbering
coefficients = extract_coefficients(expression)
@@ -234,7 +250,8 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *,
# Create a fake coordinate coefficient for a domain.
coords_coefficient = ufl.Coefficient(ufl.FunctionSpace(domain, domain.ufl_coordinate_element()))
builder.domain_coordinate[domain] = coords_coefficient
builder.set_cell_sizes(domain)
builder.set_cell_orientations((domain, ))
builder.set_cell_sizes((domain, ))
coefficients = [coords_coefficient] + coefficients
needs_external_coords = True
builder.set_coefficients(coefficients)
@@ -250,7 +267,7 @@ def compile_expression_dual_evaluation(expression, to_element, ufl_element, *,
ufl_cell=domain.ufl_cell(),
# FIXME: change if we ever implement
# interpolation on facets.
integral_type="cell",
domain_integral_type_map={domain: "cell"},
argument_multiindices=argument_multiindices,
index_cache={},
scalar_type=parameters["scalar_type"])
81 changes: 49 additions & 32 deletions tsfc/fem.py
Original file line number Diff line number Diff line change
@@ -25,10 +25,11 @@
PositiveRestricted, QuadratureWeight,
ReferenceCellEdgeVectors, ReferenceCellVolume,
ReferenceFacetVolume, ReferenceNormal,
SpatialCoordinate)
SingleValueRestricted, SpatialCoordinate)
from ufl.corealg.map_dag import map_expr_dag, map_expr_dags
from ufl.corealg.multifunction import MultiFunction
from ufl.domain import extract_unique_domain
from ufl.algorithms import extract_arguments

from tsfc import ufl2gem
from tsfc.finatinterface import as_fiat_cell, create_element
@@ -46,7 +47,7 @@ class ContextBase(ProxyKernelInterface):

keywords = ('ufl_cell',
'fiat_cell',
'integral_type',
'domain_integral_type_map',
'integration_dim',
'entity_ids',
'argument_multiindices',
@@ -80,7 +81,7 @@ def epsilon(self):
def complex_mode(self):
return is_complex(self.scalar_type)

def entity_selector(self, callback, restriction):
def entity_selector(self, callback, domain, restriction):
"""Selects code for the correct entity at run-time. Callback
generates code for a specified entity.
@@ -94,7 +95,7 @@ def entity_selector(self, callback, restriction):
if len(self.entity_ids) == 1:
return callback(self.entity_ids[0])
else:
f = self.entity_number(restriction)
f = self.entity_number(domain, restriction)
return gem.select_expression(list(map(callback, self.entity_ids)), f)

argument_multiindices = ()
@@ -131,7 +132,8 @@ def preprocess(self, expr, context):
:arg context: The translation context.
:returns: A new UFL expression
"""
ifacet = self.interface.integral_type.startswith("interior_facet")
domain = extract_unique_domain(self.mt.terminal)
ifacet = self.interface.domain_integral_type_map[domain].startswith("interior_facet")
return preprocess_expression(expr, complex_mode=context.complex_mode,
do_apply_restrictions=ifacet)

@@ -143,7 +145,7 @@ def config(self):
return config

def cell_size(self):
return self.interface.cell_size(self.mt.restriction)
return self.interface.cell_size(extract_unique_domain(self.mt.terminal), self.mt.restriction)

def jacobian_at(self, point):
ps = PointSingleton(point)
@@ -153,6 +155,10 @@ def jacobian_at(self, point):
expr = PositiveRestricted(expr)
elif self.mt.restriction == '-':
expr = NegativeRestricted(expr)
elif self.mt.restriction == '|':
expr = SingleValueRestricted(expr)
elif self.mt.restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
config = {"point_set": PointSingleton(point)}
config.update(self.config)
context = PointSetContext(**config)
@@ -165,6 +171,10 @@ def detJ_at(self, point):
expr = PositiveRestricted(expr)
elif self.mt.restriction == '-':
expr = NegativeRestricted(expr)
elif self.mt.restriction == '|':
expr = SingleValueRestricted(expr)
elif self.mt.restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
config = {"point_set": PointSingleton(point)}
config.update(self.config)
context = PointSetContext(**config)
@@ -208,6 +218,10 @@ def physical_edge_lengths(self):
expr = PositiveRestricted(expr)
elif self.mt.restriction == '-':
expr = NegativeRestricted(expr)
elif self.mt.restriction == '|':
expr = SingleValueRestricted(expr)
elif self.mt.restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")

cell = self.interface.fiat_cell
sd = cell.get_spatial_dimension()
@@ -232,6 +246,10 @@ def physical_points(self, point_set, entity=None):
expr = PositiveRestricted(expr)
elif self.mt.restriction == '-':
expr = NegativeRestricted(expr)
elif self.mt.restriction == '|':
expr = SingleValueRestricted(expr)
elif self.mt.restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
config = {"point_set": point_set}
config.update(self.config)
if entity is not None:
@@ -328,14 +346,15 @@ def __init__(self, context):
# Can't put these in the ufl2gem mixin, since they (unlike
# everything else) want access to the translation context.
def cell_avg(self, o):
if self.context.integral_type != "cell":
domain = extract_unique_domain(o)
integral_type = self.context.domain_integral_type_map[domain]
if integral_type != "cell":
# Need to create a cell-based quadrature rule and
# translate the expression using that (c.f. CellVolume
# below).
raise NotImplementedError("CellAvg on non-cell integrals not yet implemented")
integrand, = o.ufl_operands
domain = extract_unique_domain(o)
measure = ufl.Measure(self.context.integral_type, domain=domain)
measure = ufl.Measure(integral_type, domain=domain)
integrand, degree, argument_multiindices = entity_avg(integrand / CellVolume(domain), measure, self.context.argument_multiindices)

config = {name: getattr(self.context, name)
@@ -346,17 +365,18 @@ def cell_avg(self, o):
return expr

def facet_avg(self, o):
if self.context.integral_type == "cell":
domain = extract_unique_domain(o)
integral_type = self.context.domain_integral_type_map[domain]
if integral_type == "cell":
raise ValueError("Can't take FacetAvg in cell integral")
integrand, = o.ufl_operands
domain = extract_unique_domain(o)
measure = ufl.Measure(self.context.integral_type, domain=domain)
measure = ufl.Measure(integral_type, domain=domain)
integrand, degree, argument_multiindices = entity_avg(integrand / FacetArea(domain), measure, self.context.argument_multiindices)

config = {name: getattr(self.context, name)
for name in ["ufl_cell", "index_cache", "scalar_type",
"integration_dim", "entity_ids",
"integral_type"]}
"domain_integral_type_map"]}
config.update(quadrature_degree=degree, interface=self.context,
argument_multiindices=argument_multiindices)
expr, = compile_ufl(integrand, PointSetContext(**config), point_sum=True)
@@ -393,7 +413,7 @@ def translate_geometricquantity(terminal, mt, ctx):

@translate.register(CellOrientation)
def translate_cell_orientation(terminal, mt, ctx):
return ctx.cell_orientation(mt.restriction)
return ctx.cell_orientation(extract_unique_domain(terminal), mt.restriction)


@translate.register(ReferenceCellVolume)
@@ -403,7 +423,7 @@ def translate_reference_cell_volume(terminal, mt, ctx):

@translate.register(ReferenceFacetVolume)
def translate_reference_facet_volume(terminal, mt, ctx):
assert ctx.integral_type != "cell"
assert ctx.domain_integral_type_map[extract_unique_domain(terminal)] != "cell"
# Sum of quadrature weights is entity volume
return gem.optimise.aggressive_unroll(gem.index_sum(ctx.weight_expr,
ctx.point_indices))
@@ -417,7 +437,7 @@ def translate_cell_facet_jacobian(terminal, mt, ctx):

def callback(entity_id):
return gem.Literal(make_cell_facet_jacobian(cell, facet_dim, entity_id))
return ctx.entity_selector(callback, mt.restriction)
return ctx.entity_selector(callback, extract_unique_domain(terminal), mt.restriction)


def make_cell_facet_jacobian(cell, facet_dim, facet_i):
@@ -442,7 +462,7 @@ def translate_reference_normal(terminal, mt, ctx):
def callback(facet_i):
n = ctx.fiat_cell.compute_reference_normal(ctx.integration_dim, facet_i)
return gem.Literal(n)
return ctx.entity_selector(callback, mt.restriction)
return ctx.entity_selector(callback, extract_unique_domain(terminal), mt.restriction)


@translate.register(ReferenceCellEdgeVectors)
@@ -475,7 +495,7 @@ def callback(entity_id):
data = numpy.asarray(list(map(t, ps.points)))
return gem.Literal(data.reshape(point_shape + data.shape[1:]))

return gem.partial_indexed(ctx.entity_selector(callback, mt.restriction),
return gem.partial_indexed(ctx.entity_selector(callback, extract_unique_domain(terminal), mt.restriction),
ps.indices)


@@ -526,9 +546,10 @@ def translate_cellvolume(terminal, mt, ctx):

@translate.register(FacetArea)
def translate_facetarea(terminal, mt, ctx):
assert ctx.integral_type != 'cell'
domain = extract_unique_domain(terminal)
integrand, degree = one_times(ufl.Measure(ctx.integral_type, domain=domain))
integral_type = ctx.domain_integral_type_map[domain]
assert integral_type != 'cell'
integrand, degree = one_times(ufl.Measure(integral_type, domain=domain))

config = {name: getattr(ctx, name)
for name in ["ufl_cell", "integration_dim", "scalar_type",
@@ -628,7 +649,7 @@ def callback(entity_id):
# A numerical hack that FFC used to apply on FIAT tables still
# lives on after ditching FFC and switching to FInAT.
return ffc_rounding(square, ctx.epsilon)
table = ctx.entity_selector(callback, mt.restriction)
table = ctx.entity_selector(callback, extract_unique_domain(terminal), mt.restriction)
return gem.ComponentTensor(gem.Indexed(table, argument_multiindex + sigma), sigma)


@@ -668,7 +689,7 @@ def take_singleton(xs):
per_derivative = {alpha: take_singleton(tables)
for alpha, tables in per_derivative.items()}
else:
f = ctx.entity_number(mt.restriction)
f = ctx.entity_number(extract_unique_domain(terminal), mt.restriction)
per_derivative = {alpha: gem.select_expression(tables, f)
for alpha, tables in per_derivative.items()}

@@ -701,27 +722,23 @@ def take_singleton(xs):
return result


def compile_ufl(expression, context, interior_facet=False, point_sum=False):
def compile_ufl(expression, context, point_sum=False):
"""Translate a UFL expression to GEM.
:arg expression: The UFL expression to compile.
:arg context: translation context - either a :class:`GemPointContext`
or :class:`PointSetContext`
:arg interior_facet: If ``true``, treat expression as an interior
facet integral (default ``False``)
:arg point_sum: If ``true``, return a `gem.IndexSum` of the final
gem expression along the ``context.point_indices`` (if present).
"""

# Abs-simplification
expression = simplify_abs(expression, context.complex_mode)
if interior_facet:
expressions = []
for rs in itertools.product(("+", "-"), repeat=len(context.argument_multiindices)):
expressions.append(map_expr_dag(PickRestriction(*rs), expression))
else:
expressions = [expression]

arguments = extract_arguments(expression)
domains = [extract_unique_domain(argument) for argument in arguments]
integral_types = [context.domain_integral_type_map[domain] for domain in domains]
rs_tuples = [("+", "-") if integral_type.startswith("interior_facet") else (None, ) for integral_type in integral_types]
expressions = [map_expr_dag(PickRestriction(*rs), expression) for rs in itertools.product(*rs_tuples)]
# Translate UFL to GEM, lowering finite element specific nodes
result = map_expr_dags(context.translator, expressions)
if point_sum:
6 changes: 3 additions & 3 deletions tsfc/kernel_interface/__init__.py
Original file line number Diff line number Diff line change
@@ -22,15 +22,15 @@ def constant(self, const):
"""Return the GEM expression corresponding to the constant."""

@abstractmethod
def cell_orientation(self, restriction):
def cell_orientation(self, domain, restriction):
"""Cell orientation as a GEM expression."""

@abstractmethod
def cell_size(self, restriction):
def cell_size(self, domain, restriction):
"""Mesh cell size as a GEM expression. Shape (nvertex, ) in FIAT vertex ordering."""

@abstractmethod
def entity_number(self, restriction):
def entity_number(self, domain, restriction):
"""Facet or vertex number as a GEM index."""

@abstractmethod
112 changes: 61 additions & 51 deletions tsfc/kernel_interface/common.py
Original file line number Diff line number Diff line change
@@ -3,6 +3,10 @@
import string
from functools import reduce
from itertools import chain, product
import copy

from ufl.utils.sequences import max_degree
from ufl.domain import extract_unique_domain

import gem
import gem.impero_utils as impero_utils
@@ -15,24 +19,18 @@
from gem.optimise import remove_componenttensors as prune
from gem.utils import cached_property
from numpy import asarray
from tsfc import fem, ufl_utils
from tsfc import fem
from tsfc.finatinterface import as_fiat_cell, create_element
from tsfc.kernel_interface import KernelInterface
from tsfc.logging import logger
from ufl.utils.sequences import max_degree


class KernelBuilderBase(KernelInterface):
"""Helper class for building local assembly kernels."""

def __init__(self, scalar_type, interior_facet=False):
"""Initialise a kernel builder.
:arg interior_facet: kernel accesses two cells
"""
assert isinstance(interior_facet, bool)
def __init__(self, scalar_type):
"""Initialise a kernel builder."""
self.scalar_type = scalar_type
self.interior_facet = interior_facet

self.prepare = []
self.finalise = []
@@ -56,40 +54,49 @@ def coordinate(self, domain):
def coefficient(self, ufl_coefficient, restriction):
"""A function that maps :class:`ufl.Coefficient`s to GEM
expressions."""
if restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
kernel_arg = self.coefficient_map[ufl_coefficient]
domain = extract_unique_domain(ufl_coefficient)
if ufl_coefficient.ufl_element().family() == 'Real':
return kernel_arg
elif not self.interior_facet:
elif not self._domain_integral_type_map[domain].startswith("interior_facet"): # '|' is for exterior_facet
return kernel_arg
else:
return kernel_arg[{'+': 0, '-': 1}[restriction]]

def constant(self, const):
return self.constant_map[const]

def cell_orientation(self, restriction):
def cell_orientation(self, domain, restriction):
"""Cell orientation as a GEM expression."""
f = {None: 0, '+': 0, '-': 1}[restriction]
# Assume self._cell_orientations tuple is set up at this point.
co_int = self._cell_orientations[f]
if restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
if not hasattr(self, "_cell_orientations"):
raise RuntimeError("Haven't called set_cell_orientations")
f = {None: 0, '|': 0, '+': 0, '-': 1}[restriction]
co_int = self._cell_orientations[domain][f]
return gem.Conditional(gem.Comparison("==", co_int, gem.Literal(1)),
gem.Literal(-1),
gem.Conditional(gem.Comparison("==", co_int, gem.Zero()),
gem.Literal(1),
gem.Literal(numpy.nan)))

def cell_size(self, restriction):
def cell_size(self, domain, restriction):
if restriction == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
if not hasattr(self, "_cell_sizes"):
raise RuntimeError("Haven't called set_cell_sizes")
if self.interior_facet:
return self._cell_sizes[{'+': 0, '-': 1}[restriction]]
if self._domain_integral_type_map[domain].startswith("interior_facet"):
return self._cell_sizes[domain][{'+': 0, '-': 1}[restriction]]
else:
return self._cell_sizes
return self._cell_sizes[domain]

def entity_number(self, restriction):
def entity_number(self, domain, restriction):
"""Facet or vertex number as a GEM index."""
# Assume self._entity_number dict is set up at this point.
return self._entity_number[restriction]
if not hasattr(self, "_entity_numbers"):
raise RuntimeError("Haven't called set_entity_numbers")
return self._entity_numbers[domain][restriction]

def apply_glue(self, prepare=None, finalise=None):
"""Append glue code for operations that are not handled in the
@@ -127,21 +134,18 @@ def compile_integrand(self, integrand, params, ctx):
See :meth:`create_context` for typical calling sequence.
"""
# Split Coefficients
if self.coefficient_split:
integrand = ufl_utils.split_coefficients(integrand, self.coefficient_split)
# Compile: ufl -> gem
info = self.integral_data_info
functions = list(info.arguments) + [self.coordinate(info.domain)] + list(info.coefficients)
set_quad_rule(params, info.domain.ufl_cell(), info.integral_type, functions)
quad_rule = params["quadrature_rule"]
config = self.fem_config()
config['domain_integral_type_map'] = self._domain_integral_type_map
config['argument_multiindices'] = self.argument_multiindices
config['quadrature_rule'] = quad_rule
config['index_cache'] = ctx['index_cache']
expressions = fem.compile_ufl(integrand,
fem.PointSetContext(**config),
interior_facet=self.interior_facet)
fem.PointSetContext(**config))
ctx['quadrature_indices'].extend(quad_rule.point_set.indices)
return expressions

@@ -214,7 +218,7 @@ def compile_gem(self, ctx):
oriented, needs_cell_sizes, tabulations = self.register_requirements(expressions)

# Extract Variables that are actually used
active_variables = gem.extract_type(expressions, gem.Variable)
active_variables = gem.extract_variables(expressions)
# Construct ImperoC
assignments = list(zip(return_variables, expressions))
index_ordering = get_index_ordering(ctx['quadrature_indices'], return_variables)
@@ -238,7 +242,6 @@ def fem_config(self):
integration_dim, entity_ids = lower_integral_type(fiat_cell, integral_type)
return dict(interface=self,
ufl_cell=cell,
integral_type=integral_type,
integration_dim=integration_dim,
entity_ids=entity_ids,
scalar_type=self.fem_scalar_type)
@@ -429,9 +432,9 @@ def check_requirements(ir):
rt_tabs = {}
for node in traversal(ir):
if isinstance(node, gem.Variable):
if node.name == "cell_orientations":
if node.name == "cell_orientations_0":
cell_orientations = True
elif node.name == "cell_sizes":
elif node.name == "cell_sizes_0":
cell_sizes = True
elif node.name.startswith("rt_"):
rt_tabs[node.name] = node.shape
@@ -452,7 +455,7 @@ def prepare_constant(constant, number):
constant.ufl_shape)


def prepare_coefficient(coefficient, name, interior_facet=False):
def prepare_coefficient(coefficient, name, domain_integral_type_map):
"""Bridges the kernel interface and the GEM abstraction for
Coefficients.
@@ -463,30 +466,32 @@ def prepare_coefficient(coefficient, name, interior_facet=False):
expression - GEM expression referring to the Coefficient
values
"""
assert isinstance(interior_facet, bool)

if coefficient.ufl_element().family() == 'Real':
# Constant
value_size = coefficient.ufl_element().value_size
expression = gem.reshape(gem.Variable(name, (value_size,)),
coefficient.ufl_shape)
return expression

finat_element = create_element(coefficient.ufl_element())
shape = finat_element.index_shape
size = numpy.prod(shape, dtype=int)

if not interior_facet:
expression = gem.reshape(gem.Variable(name, (size,)), shape)
else:
domain = extract_unique_domain(coefficient)
integral_type = domain_integral_type_map[domain]
if integral_type is None:
# This means that this coefficient does not exist in the DAG,
# so corresponding gem expression will never be needed.
expression = None
elif integral_type.startswith("interior_facet"):
varexp = gem.Variable(name, (2 * size,))
plus = gem.view(varexp, slice(size))
minus = gem.view(varexp, slice(size, 2 * size))
expression = (gem.reshape(plus, shape), gem.reshape(minus, shape))
else:
expression = gem.reshape(gem.Variable(name, (size,)), shape)
return expression


def prepare_arguments(arguments, multiindices, interior_facet=False, diagonal=False):
def prepare_arguments(arguments, multiindices, domain_integral_type_map, diagonal=False):
"""Bridges the kernel interface and the GEM abstraction for
Arguments. Vector Arguments are rearranged here for interior
facet integrals.
@@ -499,8 +504,8 @@ def prepare_arguments(arguments, multiindices, interior_facet=False, diagonal=Fa
expressions - GEM expressions referring to the argument
tensor
"""
assert isinstance(interior_facet, bool)

if len(multiindices) != len(arguments):
raise ValueError(f"Got inconsistent lengths of arguments ({len(arguments)}) and multiindices ({len(multiindices)})")
if len(arguments) == 0:
# No arguments
expression = gem.Indexed(gem.Variable("A", (1,)), (0,))
@@ -526,15 +531,20 @@ def expression(restricted):
tuple(chain(*multiindices)))

u_shape = numpy.array([numpy.prod(shape, dtype=int) for shape in shapes])
if interior_facet:
c_shape = tuple(2 * u_shape)
slicez = [[slice(r * s, (r + 1) * s)
for r, s in zip(restrictions, u_shape)]
for restrictions in product((0, 1), repeat=len(arguments))]
else:
c_shape = tuple(u_shape)
slicez = [[slice(s) for s in u_shape]]

varexp = gem.Variable("A", c_shape)
c_shape = copy.deepcopy(u_shape)
rs_tuples = []
for arg_num, arg in enumerate(arguments):
integral_type = domain_integral_type_map[extract_unique_domain(arg)]
if integral_type is None:
raise RuntimeError(f"Can not determine integral_type on {arg}")
if integral_type.startswith("interior_facet"):
rs_tuples.append((0, 1))
c_shape[arg_num] *= 2
else:
rs_tuples.append((0, ))
slicez = [[slice(r * s, (r + 1) * s)
for r, s in zip(restrictions, u_shape)]
for restrictions in product(*rs_tuples)]
varexp = gem.Variable("A", tuple(c_shape))
expressions = [expression(gem.view(varexp, *slices)) for slices in slicez]
return tuple(prune(expressions))
262 changes: 146 additions & 116 deletions tsfc/kernel_interface/firedrake_loopy.py

Large diffs are not rendered by default.

8 changes: 6 additions & 2 deletions tsfc/modified_terminals.py
Original file line number Diff line number Diff line change
@@ -21,7 +21,7 @@
"""Definitions of 'modified terminals', a core concept in uflacs."""

from ufl.classes import (ReferenceValue, ReferenceGrad,
NegativeRestricted, PositiveRestricted,
NegativeRestricted, PositiveRestricted, SingleValueRestricted, ToBeRestricted,
Restricted, ConstantValue,
Jacobian, SpatialCoordinate, Zero)
from ufl.checks import is_cellwise_constant
@@ -39,7 +39,7 @@ class ModifiedTerminal(object):
terminal - the underlying Terminal object
local_derivatives - tuple of ints, each meaning derivative in that local direction
reference_value - bool, whether this is represented in reference frame
restriction - None, '+' or '-'
restriction - None, '+', '-', '|', or '?'
"""

def __init__(self, expr, terminal, local_derivatives, restriction, reference_value):
@@ -175,5 +175,9 @@ def construct_modified_terminal(mt, terminal):
expr = PositiveRestricted(expr)
elif mt.restriction == '-':
expr = NegativeRestricted(expr)
elif mt.restriction == '|':
expr = SingleValueRestricted(expr)
elif mt.restriction == '?':
expr = ToBeRestricted(expr)

return expr
17 changes: 15 additions & 2 deletions tsfc/ufl_utils.py
Original file line number Diff line number Diff line change
@@ -43,8 +43,11 @@ def compute_form_data(form,
do_apply_integral_scaling=True,
do_apply_geometry_lowering=True,
preserve_geometry_types=preserve_geometry_types,
do_apply_default_restrictions=True,
do_apply_restrictions=True,
do_estimate_degrees=True,
do_split_coefficients=None,
do_assume_single_integral_type=False,
complex_mode=False):
"""Preprocess UFL form in a format suitable for TSFC. Return
form data.
@@ -59,8 +62,11 @@ def compute_form_data(form,
do_apply_integral_scaling=do_apply_integral_scaling,
do_apply_geometry_lowering=do_apply_geometry_lowering,
preserve_geometry_types=preserve_geometry_types,
do_apply_default_restrictions=do_apply_default_restrictions,
do_apply_restrictions=do_apply_restrictions,
do_estimate_degrees=do_estimate_degrees,
do_split_coefficients=do_split_coefficients,
do_assume_single_integral_type=do_assume_single_integral_type,
complex_mode=complex_mode
)
constants = extract_firedrake_constants(form)
@@ -166,6 +172,8 @@ def _modified_terminal(self, o):

positive_restricted = _modified_terminal
negative_restricted = _modified_terminal
single_value_restricted = _modified_terminal
to_be_restricted = _modified_terminal

reference_grad = _modified_terminal
reference_value = _modified_terminal
@@ -250,8 +258,13 @@ def modified_terminal(self, o):
mt = analyse_modified_terminal(o)
t = mt.terminal
r = mt.restriction
if isinstance(t, Argument) and r != self.restrictions[t.number()]:
return Zero(o.ufl_shape, o.ufl_free_indices, o.ufl_index_dimensions)
if r == '?':
raise RuntimeError("Not expecting '?' restriction at this stage")
if isinstance(t, Argument) and r in ['+', '-']:
if r == self.restrictions[t.number()]:
return o
else:
return Zero(o.ufl_shape, o.ufl_free_indices, o.ufl_index_dimensions)
else:
return o

0 comments on commit 9be8132

Please sign in to comment.