-
Dominic Kempf authoredDominic Kempf authored
geometry.py 9.92 KiB
from dune.perftool import Restriction
from dune.perftool.pdelab import restricted_name
from dune.perftool.generation import (cached,
domain,
get_global_context_value,
iname,
preamble,
temporary_variable,
)
from dune.perftool.pdelab.quadrature import (name_quadrature_position,
name_quadrature_position_in_cell,
quadrature_preamble,
)
from ufl.algorithms import MultiFunction
from pymbolic.primitives import Variable
@iname
def _dimension_iname(context, count):
if context:
context = '_' + context
name = 'idim{}{}'.format(context, str(count))
formdata = get_global_context_value('formdata')
dim = formdata.geometric_dimension
domain(name, dim)
return name
def dimension_iname(context='', count=0):
return _dimension_iname(context, count)
def name_element_geometry_wrapper():
return 'eg'
def type_element_geometry_wrapper():
return 'EG'
def name_intersection_geometry_wrapper():
return 'ig'
def type_intersection_geometry_wrapper():
return 'IG'
def name_geometry_wrapper():
""" Selects the 'native' geometry wrapper of the kernel """
from dune.perftool.generation import get_global_context_value
it = get_global_context_value("integral_type")
if it == 'cell':
return name_element_geometry_wrapper()
if it == 'exterior_facet' or it == 'interior_facet':
return name_intersection_geometry_wrapper()
assert False
def type_geometry_wrapper():
""" Selects the 'native' geometry wrapper of the kernel """
from dune.perftool.generation import get_global_context_value
it = get_global_context_value("integral_type")
if it == 'cell':
return type_element_geometry_wrapper()
if it == 'exterior_facet' or it == 'interior_facet':
return type_intersection_geometry_wrapper()
assert False
@preamble
def define_restricted_cell(name, restriction):
ig = name_intersection_geometry_wrapper()
which = "inside" if restriction == Restriction.NEGATIVE else "outside"
return "const auto& {} = {}.{}();".format(name,
ig,
which,
)
def name_cell(restriction):
if restriction == Restriction.NONE:
eg = name_element_geometry_wrapper()
return "{}.entity()".format(eg)
else:
which = "inside" if restriction == Restriction.NEGATIVE else "outside"
name = "{}_cell".format(which)
define_restricted_cell(name, restriction)
return name
@preamble
def define_cell_geometry(name, restriction):
cell = name_cell(restriction)
return "auto {} = {}.geometry();".format(name,
cell
)
def name_cell_geometry(restriction):
name = restricted_name("cell_geo", restriction)
define_cell_geometry(name, restriction)
return name
def type_cell_geometry(restriction):
if restriction == Restriction.NONE:
eg = type_element_geometry_wrapper()
return "{}::Geometry".format(eg)
else:
ig = type_intersection_geometry_wrapper()
return "{}::Entity::Geometry".format(ig)
@preamble
def define_intersection_geometry(name):
ig = name_intersection_geometry_wrapper()
return "auto {} = {}.geometry();".format(name,
ig,
)
def name_intersection_geometry():
define_intersection_geometry("is_geo")
return "is_geo"
def name_intersection():
ig = name_intersection_geometry_wrapper()
return "{}.intersection()".format(ig)
def name_geometry():
""" Selects the 'native' geometry of the kernel """
from dune.perftool.generation import get_global_context_value
it = get_global_context_value("integral_type")
if it == 'cell':
return name_cell_geometry(Restriction.NONE)
if it == 'exterior_facet' or it == 'interior_facet':
return name_intersection_geometry()
assert False
@preamble
def define_in_cell_geometry(restriction, name):
ig = name_intersection_geometry_wrapper()
which = "In" if restriction == Restriction.NEGATIVE else "Out"
return "auto {} = {}.geometryIn{}side();".format(name,
ig,
which
)
def name_in_cell_geometry(restriction):
assert restriction is not Restriction.NONE
name = "geo_in_{}side".format("in" if restriction is Restriction.NEGATIVE else "out")
define_in_cell_geometry(restriction, name)
return name
def apply_in_cell_transformation(name, local, restriction):
geo = name_in_cell_geometry(restriction)
return quadrature_preamble("{} = {}.global({});".format(name,
geo,
local,
),
assignees=frozenset({name}),
read_variables=frozenset({local}),
)
def name_in_cell_coordinates(local, basename, restriction):
name = "{}_in_{}side".format(basename, "in" if restriction is Restriction.NEGATIVE else "out")
temporary_variable(name, shape=(name_dimension(),), shape_impl=("fv",))
apply_in_cell_transformation(name, local, restriction)
return name
def to_cell_coordinates(local, basename, restriction):
if restriction == Restriction.NONE:
return local
else:
return name_in_cell_coordinates(local, basename, restriction)
@preamble
def define_dimension(name):
geo = type_geometry_wrapper()
return 'const int {} = {}::Entity::dimension;'.format(name, geo)
def name_dimension():
define_dimension("dim")
return "dim"
@preamble
def define_intersection_dimension(name):
geo = type_intersection_geometry_wrapper()
return "const int {} = {}::Geometry::mydimension;".format(name, geo)
def name_intersection_dimension():
define_intersection_dimension("idim")
return "idim"
def evaluate_unit_outer_normal(name):
ig = name_intersection_geometry_wrapper()
qp = name_quadrature_position()
return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
assignees=frozenset({name}),
read_variables=frozenset({qp}),
)
@preamble
def declare_normal(name, shape, shape_impl):
ig = name_intersection_geometry_wrapper()
return "auto {} = {}.centerUnitOuterNormal();".format(name, ig)
def name_unit_outer_normal():
name = "outer_normal"
temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
evaluate_unit_outer_normal(name)
return "outer_normal"
def evaluate_unit_inner_normal(name):
outer = name_unit_outer_normal()
return quadrature_preamble("auto {} = -1. * {};".format(name, outer),
assignees=frozenset({name}),
read_variables=frozenset({outer}),
)
def name_unit_inner_normal():
name = "inner_normal"
temporary_variable(name, shape=(name_dimension(),), decl_method=declare_normal)
evaluate_unit_inner_normal(name)
return "inner_normal"
def type_jacobian_inverse_transposed(restriction):
geo = type_cell_geometry(restriction)
return "typename {}::JacobianInverseTransposed".format(geo)
@cached
def define_jacobian_inverse_transposed_temporary(restriction):
@preamble
def _define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
t = type_jacobian_inverse_transposed(restriction)
return "{} {};".format(t,
name,
)
return _define_jacobian_inverse_transposed_temporary
def define_jacobian_inverse_transposed(name, restriction):
dim = name_dimension()
temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary(restriction), shape=(dim, dim))
geo = name_cell_geometry(restriction)
pos = name_quadrature_position_in_cell(restriction)
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
pos,
),
assignees=frozenset({name}),
read_variables=frozenset({pos}),
)
def name_jacobian_inverse_transposed(restriction):
name = restricted_name("jit", restriction)
define_jacobian_inverse_transposed(name, restriction)
return name
def define_jacobian_determinant(name):
temporary_variable(name, shape=())
geo = name_geometry()
pos = name_quadrature_position()
code = "{} = {}.integrationElement({});".format(name,
geo,
pos,
)
return quadrature_preamble(code,
assignees=frozenset({name}),
read_variables=frozenset({pos}),
)
def name_jacobian_determinant():
name = 'detjac'
define_jacobian_determinant(name)
return name
def name_facet_jacobian_determinant():
name = 'fdetjac'
define_jacobian_determinant(name)
return name