diff --git a/dune/perftool/sumfact/invertgeometry.hh b/dune/perftool/sumfact/invertgeometry.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a09f5e5227fd7f2d3e28e5c63a61984eba2706ef
--- /dev/null
+++ b/dune/perftool/sumfact/invertgeometry.hh
@@ -0,0 +1,19 @@
+#ifndef DUNE_PERFTOOL_SUMFACT_INVERTGEOMETRY_HH
+#define DUNE_PERFTOOL_SUMFACT_INVERTGEOMETRY_HH
+
+
+template<typename T>
+inline T invert_and_return_determinant(const T a00, const T a10, const T a01, const T a11, T jacobian_inverse[4]){
+ T det = a00 * a11 - a10 * a01;
+ assert (det != 0.0);
+
+ jacobian_inverse[0] = a11 / det;
+ jacobian_inverse[1] = -a10 / det;
+ jacobian_inverse[2] = -a01 / det;
+ jacobian_inverse[3] = a00 / det;
+
+ return std::abs(det);
+}
+
+
+#endif
diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index 8849491e8e6c31f33f1a1c5d9a1e247122158137..5c1283b2096cf8bcfe55ef107c77d1ece2f8aa5e 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -9,7 +9,7 @@ from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position
from dune.perftool.blockstructured.spaces import lfs_inames
from dune.perftool.blockstructured.basis import (pymbolic_reference_gradient,
pymbolic_basis)
-from dune.perftool.blockstructured.geometry import (pymbolic_jacobian_inverse_transposed,
+from dune.perftool.blockstructured.geometry import (pymbolic_jacobian_inverse,
pymbolic_jacobian_determinant,
pymbolic_facet_jacobian_determinant,
to_global)
@@ -55,8 +55,8 @@ class BlockStructuredInterface(PDELabInterface):
def pymbolic_jacobian_determinant(self):
return pymbolic_jacobian_determinant()
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ return pymbolic_jacobian_inverse(i, j, restriction)
def pymbolic_facet_jacobian_determinant(self):
return pymbolic_facet_jacobian_determinant()
diff --git a/python/dune/perftool/blockstructured/geometry.py b/python/dune/perftool/blockstructured/geometry.py
index 1ffbdf0417c8261d2f0882b2384ed044af1ddd03..113f5940b0389e215c7fbaf732d812a9226dd26e 100644
--- a/python/dune/perftool/blockstructured/geometry.py
+++ b/python/dune/perftool/blockstructured/geometry.py
@@ -19,7 +19,7 @@ def pymbolic_jacobian_determinant():
# scale Jacobian according to the order of the blockstructure
-def pymbolic_jacobian_inverse_transposed(i, j, restriction):
+def pymbolic_jacobian_inverse(i, j, restriction):
name_jit = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
return prim.Product((get_form_option("number_of_blocks"),
prim.Subscript(prim.Variable(name_jit), (j, i))))
diff --git a/python/dune/perftool/pdelab/__init__.py b/python/dune/perftool/pdelab/__init__.py
index d3153ca0b8612b5391843a8689fa2c1c9d16bc56..f53b6ff6297efdde79344c49579e466ccd99d817 100644
--- a/python/dune/perftool/pdelab/__init__.py
+++ b/python/dune/perftool/pdelab/__init__.py
@@ -19,7 +19,7 @@ from dune.perftool.pdelab.geometry import (component_iname,
pymbolic_facet_area,
pymbolic_facet_jacobian_determinant,
pymbolic_jacobian_determinant,
- pymbolic_jacobian_inverse_transposed,
+ pymbolic_jacobian_inverse,
pymbolic_unit_inner_normal,
pymbolic_unit_outer_normal,
to_global,
@@ -125,8 +125,8 @@ class PDELabInterface(object):
def pymbolic_jacobian_determinant(self):
return pymbolic_jacobian_determinant()
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ return pymbolic_jacobian_inverse(i, j, restriction)
def pymbolic_unit_inner_normal(self):
return pymbolic_unit_inner_normal()
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index 39a46f211fb85e48384943ee46fd13880514b518..16f75495847532accedfd4e35a5a9c6a4bc2636c 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -204,14 +204,14 @@ def name_in_cell_geometry(restriction):
# TODO is it always necessary to add quadrature inames?
def apply_in_cell_transformation(name, local, restriction):
geo = name_in_cell_geometry(restriction)
- return get_backend("quadrature_preamble")("{} = {}.global({});".format(name,
- geo,
- str(local),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(local)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
- )
+ return quadrature_preamble("{} = {}.global({});".format(name,
+ geo,
+ str(local),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(local)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(local))}),
+ )
def pymbolic_in_cell_coordinates(local, restriction):
@@ -253,11 +253,11 @@ def local_dimension():
def evaluate_unit_outer_normal(name):
ig = name_intersection_geometry_wrapper()
qp = get_backend("quad_pos")()
- return get_backend("quadrature_preamble")("{} = {}.unitOuterNormal({});".format(name, ig, qp),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(qp)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
- )
+ return quadrature_preamble("{} = {}.unitOuterNormal({});".format(name, ig, qp),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(qp)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(qp))}),
+ )
@preamble
@@ -267,7 +267,7 @@ def declare_normal(name, kernel, decl_info):
def pymbolic_unit_outer_normal():
- name = "outer_normal"
+ name = "unit_outer_normal"
if not get_form_option("diagonal_transformation_matrix"):
temporary_variable(name, shape=(world_dimension(),), decl_method=declare_normal)
evaluate_unit_outer_normal(name)
@@ -340,14 +340,14 @@ def define_jacobian_inverse_transposed(name, restriction):
geo = name_cell_geometry(restriction)
pos = get_backend("qp_in_cell", selector=option_switch(("blockstructured", "sumfact")))(restriction)
- return get_backend("quadrature_preamble")("{} = {}.jacobianInverseTransposed({});".format(name,
- geo,
- str(pos),
- ),
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
+ return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
+ geo,
+ str(pos),
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
@backend(interface="name_jit", name="default")
@@ -357,10 +357,11 @@ def name_jacobian_inverse_transposed(restriction):
return name
-def pymbolic_jacobian_inverse_transposed(i, j, restriction):
+def pymbolic_jacobian_inverse(i, j, restriction):
+ # Calculate jacobian_inverse_transposed
name = get_backend(interface="name_jit", selector=option_switch("constant_transformation_matrix"))(restriction)
- # Dune only has JacobianInverseTransposed as a first class citizen,
- # so we need to switch the indices around!
+
+ # Return jacobian inverse -> (j, i) instead of (i, j)
return prim.Subscript(prim.Variable(name), (j, i))
@@ -397,11 +398,11 @@ def define_jacobian_determinant(name):
str(pos),
)
- return get_backend("quadrature_preamble")(code,
- assignees=frozenset({name}),
- read_variables=frozenset({get_pymbolic_basename(pos)}),
- depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
- )
+ return quadrature_preamble(code,
+ assignees=frozenset({name}),
+ read_variables=frozenset({get_pymbolic_basename(pos)}),
+ depends_on=frozenset({Writes(get_pymbolic_basename(pos))}),
+ )
@backend(interface="fdetjac", name="constant_transformation_matrix")
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index 6100879930e8eca6d88054538df5084a365d76f3..031d97b019df0cab8355d4295b5229584e077b9b 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -61,7 +61,6 @@ def quadrature_inames():
return (quadrature_iname(),)
-@backend(interface="quadrature_preamble")
def quadrature_preamble(code, **kw):
kw['tags'] = kw.get('tags', frozenset({})).union(frozenset({"quad"}))
return instruction(inames=get_backend(interface="quad_inames")(), code=code, **kw)
diff --git a/python/dune/perftool/sumfact/__init__.py b/python/dune/perftool/sumfact/__init__.py
index d5eed5544f52c76be02cce4e59455d5a53770ba1..d3192b6f278bc7e0d050a500da86a09d905a4a7e 100644
--- a/python/dune/perftool/sumfact/__init__.py
+++ b/python/dune/perftool/sumfact/__init__.py
@@ -78,37 +78,25 @@ class SumFactInterface(PDELabInterface):
return ret
def pymbolic_unit_outer_normal(self):
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
- ret, indices = pymbolic_unit_outer_normal(self.visitor.indices)
- self.visitor.indices = indices
- return ret
+ from dune.perftool.sumfact.geometry import pymbolic_unit_outer_normal
+ return pymbolic_unit_outer_normal(self.visitor)
def pymbolic_unit_inner_normal(self):
from dune.perftool.sumfact.geometry import pymbolic_unit_inner_normal
- ret, indices = pymbolic_unit_inner_normal(self.visitor.indices)
- self.visitor.indices = indices
- return ret
+ return pymbolic_unit_inner_normal(self.visitor)
def pymbolic_facet_jacobian_determinant(self):
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
- return pymbolic_facet_jacobian_determinant()
+ from dune.perftool.sumfact.geometry import pymbolic_facet_jacobian_determinant
+ return pymbolic_facet_jacobian_determinant(self.visitor)
def pymbolic_facet_area(self):
from dune.perftool.sumfact.geometry import pymbolic_facet_area
- return pymbolic_facet_area()
+ return pymbolic_facet_area(self.visitor)
- def pymbolic_jacobian_inverse_transposed(self, i, j, restriction):
- from dune.perftool.pdelab.geometry import pymbolic_jacobian_inverse_transposed
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- return pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ def pymbolic_jacobian_inverse(self, i, j, restriction):
+ from dune.perftool.sumfact.geometry import pymbolic_jacobian_inverse
+ return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
def pymbolic_jacobian_determinant(self):
- from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
- from dune.perftool.generation import global_context
- with global_context(visitor=self.visitor):
- return pymbolic_jacobian_determinant()
+ from dune.perftool.sumfact.geometry import pymbolic_jacobian_determinant
+ return pymbolic_jacobian_determinant(self.visitor)
diff --git a/python/dune/perftool/sumfact/accumulation.py b/python/dune/perftool/sumfact/accumulation.py
index 4e7c2a6991b1c6bda49862de0aa47167c1b1ea2a..abad6788cb31f90276dbde32e2fb5aa2b6faa3af 100644
--- a/python/dune/perftool/sumfact/accumulation.py
+++ b/python/dune/perftool/sumfact/accumulation.py
@@ -24,7 +24,6 @@ from dune.perftool.options import (get_form_option,
)
from dune.perftool.loopy.flatten import flatten_index
from dune.perftool.loopy.target import type_floatingpoint
-from dune.perftool.sumfact.quadrature import nest_quadrature_loops
from dune.perftool.pdelab.driver import FEM_name_mangling
from dune.perftool.pdelab.localoperator import determine_accumulation_space
from dune.perftool.pdelab.restriction import restricted_name
diff --git a/python/dune/perftool/sumfact/geometry.py b/python/dune/perftool/sumfact/geometry.py
index f3e5b9fde10587ffb02f802b6984b3309f9b52b5..ecec31af151da1d9bbfda0e8a6ac46ac44b8cf94 100644
--- a/python/dune/perftool/sumfact/geometry.py
+++ b/python/dune/perftool/sumfact/geometry.py
@@ -6,38 +6,71 @@ from dune.perftool.generation import (backend,
get_backend,
get_counted_variable,
iname,
+ include_file,
instruction,
kernel_cached,
preamble,
+ silenced_warning,
temporary_variable,
+ get_global_context_value,
globalarg,
+ valuearg,
)
+from dune.perftool.options import get_option
from dune.perftool.pdelab.geometry import (local_dimension,
world_dimension,
+ name_cell_geometry,
name_geometry,
)
-from dune.perftool.sumfact.switch import get_facedir
-from dune.perftool.sumfact.symbolic import SumfactKernelInterfaceBase
+from dune.perftool.pdelab.localoperator import (name_ansatz_gfs_constructor_param,
+ lop_template_ansatz_gfs,
+ lop_template_range_field,
+ )
+from dune.perftool.pdelab.restriction import restricted_name
+from dune.perftool.sumfact.basis import construct_basis_matrix_sequence
+from dune.perftool.sumfact.quadrature import (additional_inames,
+ default_quadrature_inames)
+from dune.perftool.sumfact.realization import (name_buffer_storage,
+ realize_sum_factorization_kernel,
+ )
+from dune.perftool.sumfact.switch import get_facedir, get_facemod
+from dune.perftool.sumfact.symbolic import SumfactKernelInterfaceBase, SumfactKernel
+from dune.perftool.sumfact.tabulation import (quadrature_points_per_direction,
+ BasisTabulationMatrix,
+ )
from dune.perftool.sumfact.vectorization import attach_vectorization_info
+from dune.perftool.tools import get_pymbolic_basename
from dune.perftool.options import get_form_option, option_switch
from dune.perftool.ufl.modified_terminals import Restriction
from pytools import ImmutableRecord
+from loopy.match import Writes
+
import pymbolic.primitives as prim
import numpy as np
@iname
-def corner_iname():
- name = get_counted_variable("corneriname")
- domain(name, 2 ** local_dimension())
+def global_corner_iname(restriction):
+ name = get_counted_variable(restricted_name("global_corneriname", restriction))
+ domain(name, 2 ** world_dimension())
return name
class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
- def __init__(self, dir):
- ImmutableRecord.__init__(self, dir=dir)
+ def __init__(self, direction, restriction):
+ """Base class for sum-factorized evaluation of geometry mappings
+
+ At the moment we only do this for cells and not faces. For
+ intersections we do this corresponding reference elements of the
+ neigboring cells.
+
+ Each spatial component needs a seperate sum factorization kernel. The
+ argument 'direction' specifies the component (x-component: 0,
+ y-component: 1, z-component: 2).
+ """
+ ImmutableRecord.__init__(self, direction=direction, restriction=restriction)
@property
def stage(self):
@@ -48,18 +81,22 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
def realize(self, sf, insn_dep):
# TODO: When we use vectorization this should be the offset
+ assert(get_form_option('vectorization_strategy') == 'none')
index = 0
- from dune.perftool.sumfact.realization import name_buffer_storage
+ # Note: world_dimension, since we only do evaluation of cell geometry mappings
name = "input_{}".format(sf.buffer)
temporary_variable(name,
- shape=(2 ** local_dimension(), sf.vector_width),
+ shape=(2 ** world_dimension(), sf.vector_width),
custom_base_storage=name_buffer_storage(sf.buffer, 0),
managed=True,
)
+ ciname = global_corner_iname(self.restriction)
- ciname = corner_iname()
- geo = name_geometry()
+ if self.restriction == 0:
+ geo = name_geometry()
+ else:
+ geo = name_cell_geometry(self.restriction)
# NB: We need to realize this as a C instruction, because the corner
# method does return a non-scalar, which does not fit into the current
@@ -71,7 +108,7 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
index,
geo,
ciname,
- self.dir,
+ self.direction,
)
insn = instruction(code=code,
@@ -85,23 +122,33 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
@backend(interface="spatial_coordinate", name="default")
def pymbolic_spatial_coordinate_multilinear(do_predicates, visitor):
+ """Calculate global coordinates of quadrature points for multilinear geometry mappings
+
+ The evalualation is done using sum factorization techniques.
+
+ On facets we use the geometry mapping of the self cell to get a global
+ evaluation of the quadrature points. Avoiding the geometry mapping of the
+ face itself greatly simplifies the generation of code for unstructured
+ grids. Instead of looking at all the possible embeddings of the reference
+ element of the face into the reference element of the cell we can make the
+ grid edge consistent and choose a suitable convention for the order of
+ directions in our sum factorization kernels.
+ """
assert len(visitor.indices) == 1
- # Construct the matrix sequence for the evaluation of the global coordinate.
- # We need to manually construct this one, because on facets, we want to use the
- # geometry embedding of the facet into the global space directly without going
- # through the neighboring cell geometries. That matrix sequence will only have
- # dim-1 matrices!
- from dune.perftool.sumfact.tabulation import quadrature_points_per_direction, BasisTabulationMatrix
- quadrature_size = quadrature_points_per_direction()
- matrix_sequence = tuple(BasisTabulationMatrix(direction=i, basis_size=2) for i in range(world_dimension()) if i != get_facedir(visitor.restriction))
- inp = GeoCornersInput(visitor.indices[0])
-
- from dune.perftool.sumfact.symbolic import SumfactKernel
+ # Adjust restriction for boundary facets
+ restriction = visitor.restriction
+ if get_global_context_value("integral_type") == "exterior_facet":
+ restriction = 1
+
+ # Generate sum factorization kernel and add vectorization info
+ matrix_sequence = construct_basis_matrix_sequence(facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=(2,) * world_dimension())
+ inp = GeoCornersInput(visitor.indices[0], restriction)
sf = SumfactKernel(matrix_sequence=matrix_sequence,
interface=inp,
)
-
vsf = attach_vectorization_info(sf)
# Add a sum factorization kernel that implements the evaluation of
@@ -109,7 +156,10 @@ def pymbolic_spatial_coordinate_multilinear(do_predicates, visitor):
from dune.perftool.sumfact.realization import realize_sum_factorization_kernel
var, _ = realize_sum_factorization_kernel(vsf)
+ # The results of this function is already the right component of the
+ # spatial coordinate and we don't need to index this in the visitor
visitor.indices = None
+
return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
@@ -164,7 +214,6 @@ def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
# Urgh: *SOMEHOW* construct a face direction
from dune.perftool.pdelab.restriction import Restriction
restriction = Restriction.NONE
- from dune.perftool.generation import get_global_context_value
if get_global_context_value("integral_type") == "interior_facet":
restriction = Restriction.POSITIVE
from dune.perftool.sumfact.switch import get_facedir
@@ -192,69 +241,171 @@ def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,))
-def pymbolic_unit_outer_normal(visitor_indices):
- index, = visitor_indices
+def define_outer_normal(name, visitor):
+ """Calculate outer normal (length not necessarily 1)
+
+ This will be used to calculate the unit outer normal and the determinant of
+ the Jacobian of the geometry mapping of the face.
+
+ The outer normal is calculated by multiplying the Jacobian inverse
+ transposed of the geometry mapping of the 'self'-cell with the unit outer
+ normal \hat{n_s} of the face corresponding to the intersetcion in the
+ reference element of the 'self'-cell:
+
+ n = (\nabla \mu_{T_s})^{-T} \hat{n_s}
+
+ Instead of setting up \hat{n_s} we just look at facedir_s and facemod_s.
+ """
+ facedir_s = get_facedir(Restriction.POSITIVE)
+ facemod_s = get_facemod(Restriction.POSITIVE)
+ temporary_variable(name, shape=(world_dimension(),))
+ for i in range(world_dimension()):
+ assignee = prim.Subscript(prim.Variable(name), i)
+ jit = pymbolic_jacobian_inverse(i, facedir_s, Restriction.POSITIVE, visitor)
+
+ # Note: 2*facemod_s-1 because of
+ # -1 if facemod_s = 0
+ # +1 if facemod_s = 1
+ expression = jit * (2 * facemod_s - 1)
+
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+
+
+def name_outer_normal(visitor):
+ name = "outer_normal"
+ define_outer_normal(name, visitor)
+ return name
+
+
+def pymbolic_outer_normal(visitor):
+ name = name_outer_normal(visitor)
+ return prim.Variable(name)
+
+
+def define_norm_of_outer_normal(name, visitor):
+ outer_normal = pymbolic_outer_normal(visitor)
+
+ # Calculate the norm of outer_normal
+ temporary_variable(name, shape=())
+ expression = 0
+ for i in range(world_dimension()):
+ expression = prim.Sum((expression, prim.Subscript(outer_normal, i) * prim.Subscript(outer_normal, i)))
+ expression = prim.Call(prim.Variable("sqrt"), (expression,))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+
+ # Apparently the call of sqrt shadows the dependency on outer normal. Add manually.
+ instruction(assignee=prim.Variable(name),
+ expression=expression,
+ depends_on=frozenset({Writes(get_pymbolic_basename(outer_normal))}),
+ within_inames=frozenset(inames),
+ )
+
+
+def name_norm_of_outer_normal(visitor):
+ name = "norm_of_outer_normal"
+ define_norm_of_outer_normal(name, visitor)
+ return name
+
+
+def pymbolic_norm_of_outer_normal(visitor):
+ name = name_norm_of_outer_normal(visitor)
+ return prim.Variable(name)
+
+
+def pymbolic_unit_outer_normal(visitor):
+ assert (len(visitor.indices) is 1)
+ index = visitor.indices[0]
assert isinstance(index, int)
if get_form_option("diagonal_transformation_matrix"):
+ # In this case we just return -1, 0 or 1 depending on the current
+ # facemod and facedir. We don't need to (and cannot) index these ints
+ # so we set the indices of the visitor to None.
+ visitor.indices = None
+
+ # Use facemod_s and facedir_s
from dune.perftool.sumfact.switch import get_facedir, get_facemod
if index == get_facedir(Restriction.POSITIVE):
if get_facemod(Restriction.POSITIVE):
- return 1, None
+ return 1
else:
- return -1, None
+ return -1
else:
- return 0, None
+ return 0
else:
- from dune.perftool.pdelab.geometry import pymbolic_unit_outer_normal as _norm
- return _norm(), visitor_indices
+ # In pymbolic_outer_normal we need to create the
+ # jacobian_inverse_transposed which doesn't play nicely with the
+ # visitor having indices. Maybe this should be improved. For now let's
+ # just set them to None and restore them afterwards.
+ visitor.indices = None
+ outer_normal = pymbolic_outer_normal(visitor)
+ norm = pymbolic_norm_of_outer_normal(visitor)
+ visitor.indices = (index,)
+
+ # Create unit_outer_normal by scaling outer_normal to length 1
+ name = "unit_outer_normal"
+ temporary_variable(name, shape=(world_dimension(),))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ for i in range(world_dimension()):
+ instruction(assignee=prim.Subscript(prim.Variable(name), i),
+ expression=prim.Subscript(outer_normal, i) / norm,
+ within_inames=frozenset(inames),
+ )
+
+ return prim.Variable(name)
def pymbolic_unit_inner_normal(visitor_indices):
- index, = visitor_indices
+ # Note: The implementation of this was adjusted for unstructured grids but
+ # it was not properly tested. If you need the unit inner normal just remove
+ # the assert(False) statement but make sure to verify that this does the
+ # right thing.
+ assert(False)
+
+ assert (len(visitor.indices) is 1)
+ index = visitor.indices[0]
assert isinstance(index, int)
if get_form_option("diagonal_transformation_matrix"):
- from dune.perftool.sumfact.switch import get_facedir, get_facemod
+ # In this case we just return -1, 0 or 1 depending on the current
+ # facemod and facedir. We don't need to (and cannot) index these ints
+ # so we set the indices of the visitor to None.
+ visitor.indices = None
if index == get_facedir(Restriction.POSITIVE):
if get_facemod(Restriction.POSITIVE):
- return -1, None
+ return -1
else:
- return 1, None
+ return 1
else:
- return 0, None
- else:
- from dune.perftool.pdelab.geometry import pymbolic_unit_inner_normal as _norm
- return _norm(), visitor_indices
-
-
-def pymbolic_facet_jacobian_determinant():
- if get_form_option("constant_transformation_matrix"):
- return pymbolic_constant_facet_jacobian_determinant()
+ return 0
else:
- from dune.perftool.pdelab.geometry import pymbolic_facet_jacobian_determinant as _norm
- return _norm()
+ # In pymbolic_outer_normal we need to create the
+ # jacobian_inverse_transposed which doesn't play nicely with the
+ # visitor having indices. Maybe this should be improved. For now let's
+ # just set them to None and restore them afterwards.
+ visitor.indices = None
+ outer_normal = pymbolic_outer_normal(visitor)
+ norm = pymbolic_norm_of_outer_normal(visitor)
+ visitor.indices = (index,)
+ # Create unit_outer_normal by scaling outer_normal to length 1
+ name = "unit_outer_normal"
+ temporary_variable(name, shape=(world_dimension(),))
-def pymbolic_constant_facet_jacobian_determinant():
- facedir = get_facedir(Restriction.POSITIVE)
- assert isinstance(facedir, int)
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ for i in range(world_dimension()):
+ instruction(assignee=prim.Subscript(prim.Variable(name), i),
+ expression=-1 * prim.Subscript(outer_normal, i) / norm,
+ within_inames=frozenset(inames),
+ )
- name = "fdetjac"
- define_constant_facet_jacobian_determinant(name)
- globalarg(name, shape=(world_dimension(),))
- return prim.Subscript(prim.Variable(name), (facedir,))
-
-
-@class_member(classtag="operator")
-def define_constant_facet_jacobian_determinant(name):
- define_constant_facet_jacobian_determinant_eval(name)
- from dune.perftool.pdelab.localoperator import lop_template_ansatz_gfs
- gfst = lop_template_ansatz_gfs()
- return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+ return prim.Variable(name)
@kernel_cached(kernel="operator")
def define_constant_facet_jacobian_determinant_eval(name):
- from dune.perftool.pdelab.localoperator import name_ansatz_gfs_constructor_param, lop_template_range_field
gfs = name_ansatz_gfs_constructor_param()
rft = lop_template_range_field()
code = ["{",
@@ -269,9 +420,196 @@ def define_constant_facet_jacobian_determinant_eval(name):
)
-def pymbolic_facet_area():
+@class_member(classtag="operator")
+def define_constant_facet_jacobian_determinant(name):
+ define_constant_facet_jacobian_determinant_eval(name)
+ gfst = lop_template_ansatz_gfs()
+ return "std::array<typename {}::Traits::GridView::template Codim<0>::Geometry::ctype, {}> {};".format(gfst, world_dimension(), name)
+
+
+def pymbolic_constant_facet_jacobian_determinant():
+ facedir = get_facedir(Restriction.POSITIVE)
+ assert isinstance(facedir, int)
+
+ name = "fdetjac"
+ define_constant_facet_jacobian_determinant(name)
+ globalarg(name, shape=(world_dimension(),))
+ return prim.Subscript(prim.Variable(name), (facedir,))
+
+
+def define_facet_jacobian_determinant(name, visitor):
+ """Absolute value of determinant of jacobian of facet geometry mapping
+
+ Calculate the absolute value of the determinant of the jacobian of the
+ geometry mapping from the reference cell of the intersection to the
+ intersection:
+
+ |\det(\nabla \mu_F)|
+
+ This is done using the surface metric tensor lemma from the lecture notes
+ of Jean-Luc Guermond:
+
+ |\det(\nabla \mu_F)| = \| \tilde{n} \| |\det(\nabla \mu_{T_s})|
+
+ Here \tilde{n} is the outer normal defined by:
+
+ \tilde{n} = (\nabla \mu_{T_s})^{-T} \hat{n}_s
+
+ where \hat{n}_s is the unit outer normal of the corresponding face of the
+ reference cell.
+ """
+ # Calculate norm of outer normal and determinant of jacobian of geometry
+ # mapping of the inside cell -> set restriction to 1 for boundary facets
+ if get_global_context_value("integral_type") == "exterior_facet":
+ assert visitor.restriction is 0
+ visitor.restriction = 1
+ norm_of_outer_normal = pymbolic_norm_of_outer_normal(visitor)
+ detjac = pymbolic_jacobian_determinant(visitor)
+ if get_global_context_value("integral_type") == "exterior_facet":
+ visitor.restriction = 0
+
+ # Multiply detjac and norm_of_outer_normal
+ temporary_variable(name, shape=())
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=prim.Variable(name),
+ expression=detjac * norm_of_outer_normal,
+ within_inames=frozenset(inames),
+ )
+
+
+def name_facet_jacobian_determinant(visitor):
+ name = "fdetjac"
+ define_facet_jacobian_determinant(name, visitor)
+ return name
+
+
+def pymbolic_facet_jacobian_determinant(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ return pymbolic_constant_facet_jacobian_determinant()
+ else:
+ name = name_facet_jacobian_determinant(visitor)
+ return prim.Variable(name)
+
+
+def pymbolic_facet_area(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ return pymbolic_facet_jacobian_determinant(visitor)
+ else:
+ from dune.perftool.pdelab.geometry import pymbolic_facet_area as nonconstant_facet_area
+ return nonconstant_facet_area()
+
+
+def _name_jacobian(i, j, restriction, visitor):
+ """Return the (i, j) component of the jacobian of the geometry mapping
+
+ Evaluation of the derivative of the geometry mapping is done using sum
+ factorization.
+
+ Note: At the moment this only works for the mappings from reference cells
+ to the cell and not for the geometry mappings of intersections.
+ """
+ do_predicates = visitor.do_predicates
+
+ # Create matrix sequence with derivative in j direction
+ matrix_sequence = construct_basis_matrix_sequence(derivative=j,
+ facedir=get_facedir(restriction),
+ facemod=get_facemod(restriction),
+ basis_size=(2,) * world_dimension())
+
+ # Sum factorization input for the i'th component of the geometry mapping
+ inp = GeoCornersInput(i, restriction)
+
+ sf = SumfactKernel(matrix_sequence=matrix_sequence,
+ interface=inp,
+ )
+ vsf = attach_vectorization_info(sf)
+
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
+ var, _ = realize_sum_factorization_kernel(vsf)
+
+ assert(visitor.indices is None)
+ return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
+
+
+def define_jacobian_inverse(name, restriction, visitor):
+ """Return jacobian inverse transposed of the geometry mapping (of a cell)
+
+ At the moment this only works for geometry mappings of cells and not for
+ intersection. We only consider this case as it greatly simplifies code
+ generation for unstructured grids by making the grid edge consistent and
+ avoiding the face-twist problem.
+ """
+ # Calculate the jacobian of the geometry mapping using sum factorization
+ dim = world_dimension()
+ names_jacobian = []
+ for j in range(dim):
+ for i in range(dim):
+ name_jacobian = restricted_name("jacobian_geometry_{}{}".format(i, j), restriction)
+ temporary_variable(name_jacobian)
+ assignee = prim.Variable(name_jacobian)
+ expression = _name_jacobian(i, j, restriction, visitor)
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+ names_jacobian.append(name_jacobian)
+
+ # Calculate the inverse of the jacobian of the geometry mapping and the
+ # determinant by calling a c++ function
+ name_detjac = name_jacobian_determinant(visitor)
+ temporary_variable(name_detjac, shape=())
+ temporary_variable(name, shape=(dim, dim), managed=True)
+ include_file('dune/perftool/sumfact/invertgeometry.hh', filetag='operatorfile')
+ code = "{} = invert_and_return_determinant({}, {});".format(name_detjac,
+ ", ".join(names_jacobian),
+ name)
+ silenced_warning("read_no_write({})".format(name_detjac))
+ inames = default_quadrature_inames(visitor) + additional_inames(visitor)
+ return instruction(code=code,
+ assignees=frozenset([name, name_detjac]),
+ read_variables=frozenset(names_jacobian),
+ inames=frozenset(inames),
+ )
+
+
+def name_jacobian_inverse(restriction, visitor):
+ name = restricted_name("jacobian_inverse", restriction)
+ define_jacobian_inverse(name, restriction, visitor)
+ return name
+
+
+def pymbolic_jacobian_inverse(i, j, restriction, visitor):
+ # Switch between implementations without using the backend mechanism. In
+ # the case of constant_transformation_matrix we want to use the pdelab
+ # branch, otherwise we want to go into the sumfact implementation instead
+ # of pdelab.
if get_form_option("constant_transformation_matrix"):
- return pymbolic_facet_jacobian_determinant()
+ from dune.perftool.pdelab.geometry import name_constant_jacobian_inverse_transposed
+ name = name_constant_jacobian_inverse_transposed(restriction)
+ # Return jacobian inverse -> (j, i) instead of (i, j)
+ return prim.Subscript(prim.Variable(name), (j, i))
else:
- from dune.perftool.pdelab.geometry import pymbolic_facet_area as _norm
- return _norm()
+ name = name_jacobian_inverse(restriction, visitor)
+ return prim.Subscript(prim.Variable(name), (i, j))
+
+
+def name_jacobian_determinant(visitor):
+ name = "detjac"
+ return name
+
+
+def pymbolic_jacobian_determinant(visitor):
+ if get_form_option("constant_transformation_matrix"):
+ from dune.perftool.pdelab.geometry import pymbolic_jacobian_determinant
+ return pymbolic_jacobian_determinant()
+ else:
+ # The calculation of the jacobian determinant happens in this function
+ if visitor.measure == 'cell':
+ restriction = 0
+ else:
+ restriction = 1
+ name_jacobian_inverse(restriction, visitor)
+
+ return prim.Variable(name_jacobian_determinant(visitor))
diff --git a/python/dune/perftool/sumfact/quadrature.py b/python/dune/perftool/sumfact/quadrature.py
index 586d9fd26a823c48ec1a7a220a840a5d928e30e1..6fe70fa441979c4ae00a5f44acc8136def9f5fb0 100644
--- a/python/dune/perftool/sumfact/quadrature.py
+++ b/python/dune/perftool/sumfact/quadrature.py
@@ -39,20 +39,6 @@ import loopy as lp
import numpy as np
-def nest_quadrature_loops(kernel, inames):
- insns = []
- for insn in kernel.instructions:
- if "quad" in insn.tags:
- insns.append(insn.copy(within_inames=insn.within_inames.union(frozenset(inames)),
- tags=insn.tags - frozenset({"quad"})
- )
- )
- else:
- insns.append(insn)
-
- return kernel.copy(instructions=insns)
-
-
class BaseWeight(FunctionIdentifier):
def __init__(self, accumobj):
self.accumobj = accumobj
@@ -80,6 +66,15 @@ def pymbolic_base_weight():
return 1.0
+def default_quadrature_inames(visitor):
+ info = visitor.current_info[1]
+ if info is None:
+ element = None
+ else:
+ element = info.element
+ return quadrature_inames(element)
+
+
@kernel_cached
def quadrature_inames(element):
if element is None:
@@ -177,7 +172,8 @@ def quadrature_weight(visitor):
element = None
else:
element = info.element
- return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
+ return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"),
+ tuple(prim.Variable(i) for i in quadrature_inames(element)))
def define_quadrature_position(name, local_index):
@@ -245,7 +241,9 @@ def pymbolic_indexed_quadrature_position(local_index, visitor):
return prim.Subscript(lp.symbolic.TaggedVariable(name, "operator_precomputed"), tuple(prim.Variable(i) for i in quadrature_inames(element)))
-def _additional_inames(visitor):
+def additional_inames(visitor):
+ """Return inames for iterating over ansatz space in the case of jacobians
+ """
info = visitor.current_info[1]
if info is None:
element = None
@@ -262,59 +260,3 @@ def _additional_inames(visitor):
return lfs_inames
else:
return ()
-
-
-@backend(interface="quadrature_preamble", name="sumfact")
-def quadrature_preamble(code, **kw):
- visitor = get_global_context_value('visitor')
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
-
- inames = quadrature_inames(element)
- inames = inames + _additional_inames(visitor)
-
- return instruction(inames=inames, code=code, **kw)
-
-
-@preamble
-def define_quadrature_point(name):
- local_dim = local_dimension()
- return "Dune::FieldVector<double, {}> {};".format(local_dim, name)
-
-
-def name_quadrature_point():
- name = "quadrature_point"
- define_quadrature_point(name)
- return name
-
-
-@backend(interface="quad_pos", name="sumfact")
-def pymbolic_quadrature_position():
- visitor = get_global_context_value('visitor')
-
- info = visitor.current_info[1]
- if info is None:
- element = None
- else:
- element = info.element
-
- quad_point = name_quadrature_point()
- inames = quadrature_inames(element) + _additional_inames(visitor)
- globalarg(quad_point, shape=(local_dimension(),))
- facedir = get_facedir(visitor.restriction)
- for i in range(local_dimension()):
- expression = pymbolic_indexed_quadrature_position(i, visitor)
- instruction(assignee=prim.Subscript(prim.Variable(quad_point), i),
- expression=expression,
- within_inames=frozenset(inames))
- return Variable(quad_point)
-
-
-@backend(interface="qp_in_cell", name="sumfact")
-def pymbolic_quadrature_position_in_cell(restriction):
- from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
- return pymbolic_quadrature_position_in_cell(restriction)
diff --git a/python/dune/perftool/sumfact/symbolic.py b/python/dune/perftool/sumfact/symbolic.py
index fb283a0536318c1585a0963a33ca76105571cb84..f2fbcb76abfea8a594a1df28bebfc2142f2275aa 100644
--- a/python/dune/perftool/sumfact/symbolic.py
+++ b/python/dune/perftool/sumfact/symbolic.py
@@ -417,6 +417,23 @@ class SumfactKernel(SumfactKernelBase, ImmutableRecord, prim.Variable):
if len(self.matrix_sequence) == local_dimension():
return tuple(prim.Variable(i) for i in quad_inames)
+ # # Extremly ugly. This is the special case of being on an intersection
+ # # and generating a SF kernel for the geometry evaluation of one of the
+ # # neighboring cells. This screws the usual setup since the SF sequence
+ # # is not equal to the local dimension.
+ # #
+ # # TODO: Vectorization is probably broken for this case!
+ # count_sf_sequence_length = 0
+ # for d in range(world_dimension()):
+ # if self.matrix_sequence[d].face is not None:
+ # count_sf_sequence_length += 1
+ # if count_sf_sequence_length != local_dimension():
+ # from dune.perftool.sumfact.quadrature import global_quadrature_inames
+ # quad_inames = global_quadrature_inames(element)
+
+ # for d in range(world_dimension()):
+ # print(self.matrix_sequence[d].face)
+
# Traverse all the quadrature inames and map them to their correct direction
index = []
i = 0
diff --git a/python/dune/perftool/ufl/visitor.py b/python/dune/perftool/ufl/visitor.py
index cbede6c30fa7d6754b58b57b662ba126a1dec6b5..9e14391c5aef7656dfcc6e0fe8276ffcf3016505 100644
--- a/python/dune/perftool/ufl/visitor.py
+++ b/python/dune/perftool/ufl/visitor.py
@@ -464,7 +464,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
if isinstance(i, int) and isinstance(j, int) and i != j:
return 0
- return self.interface.pymbolic_jacobian_inverse_transposed(i, j, restriction)
+ return self.interface.pymbolic_jacobian_inverse(i, j, restriction)
def jacobian(self, o):
raise NotImplementedError("How did you get Jacobian into your form? We only support JacobianInverse right now. Report!")