From 3cb3511c6234f43d45a0e26a25e20222c40403b4 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Mon, 10 Dec 2018 11:40:31 +0100
Subject: [PATCH] Implement sumfact_multilinear Geometry Mixin
---
python/dune/codegen/options.py | 2 +-
python/dune/codegen/pdelab/__init__.py | 13 -
python/dune/codegen/sumfact/geometry.py | 577 ++++++------------
python/dune/codegen/sumfact/quadrature.py | 12 +-
python/dune/codegen/ufl/visitor.py | 6 +-
.../poisson/poisson_2d_unstructured.mini | 1 +
.../poisson/poisson_3d_unstructured.mini | 1 +
test/sumfact/poisson/poisson_dg_2d_gmsh.mini | 1 +
.../poisson/poisson_dg_2d_unstructured.mini | 1 +
test/sumfact/poisson/poisson_dg_3d_gmsh.mini | 1 +
.../poisson/poisson_dg_3d_unstructured.mini | 1 +
.../poisson/poisson_fastdg_3d_gmsh.mini | 1 +
12 files changed, 199 insertions(+), 418 deletions(-)
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 32ede8d4..4a5a913d 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -113,7 +113,7 @@ class CodegenFormOptionsArray(ImmutableRecord):
geometry_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for geometries. Currently implemented mixins: generic, axiparallel, equidistant, sumfact_multilinear, sumfact_axiparallel, sumfact_equidistant")
quadrature_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for quadrature. Currently implemented: generic, sumfact")
basis_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for basis function evaluation. Currently implemented: generic, sumfact")
- accumulation_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for accumulation. Currently implemented: generic")
+ accumulation_mixins = CodegenOption(default="generic", helpstr="A comma separated list of mixin identifiers to use for accumulation. Currently implemented: generic, sumfact")
enable_volume = CodegenOption(default=True, helpstr="Whether to assemble volume integrals")
enable_skeleton = CodegenOption(default=True, helpstr="Whether to assemble skeleton integrals")
enable_boundary = CodegenOption(default=True, helpstr="Whether to assemble boundary integrals")
diff --git a/python/dune/codegen/pdelab/__init__.py b/python/dune/codegen/pdelab/__init__.py
index 604cacf0..7ca90b96 100644
--- a/python/dune/codegen/pdelab/__init__.py
+++ b/python/dune/codegen/pdelab/__init__.py
@@ -20,19 +20,6 @@ class PDELabInterface(object):
# The visitor instance will be registered by its init method
self.visitor = None
- # Accumulation interfaces
- def get_accumulation_info(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import get_accumulation_info
- return get_accumulation_info(expr, visitor)
-
- def list_accumulation_infos(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import list_accumulation_infos
- return list_accumulation_infos(expr, visitor)
-
- def generate_accumulation_instruction(self, expr, visitor):
- from dune.codegen.pdelab.localoperator import generate_accumulation_instruction
- return generate_accumulation_instruction(expr, visitor)
-
#
# TODO: The following ones are actually entirely PDELab independent!
# They should be placed elsewhere and be used directly in the visitor.
diff --git a/python/dune/codegen/sumfact/geometry.py b/python/dune/codegen/sumfact/geometry.py
index e4c23caf..ade47c57 100644
--- a/python/dune/codegen/sumfact/geometry.py
+++ b/python/dune/codegen/sumfact/geometry.py
@@ -17,7 +17,8 @@ from dune.codegen.generation import (backend,
valuearg,
)
from dune.codegen.options import get_option
-from dune.codegen.pdelab.geometry import (local_dimension,
+from dune.codegen.pdelab.geometry import (enforce_boundary_restriction,
+ local_dimension,
world_dimension,
name_cell_geometry,
name_geometry,
@@ -32,8 +33,7 @@ from dune.codegen.pdelab.localoperator import (name_ansatz_gfs_constructor_param
from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.sumfact.accumulation import basis_sf_kernels
from dune.codegen.sumfact.basis import construct_basis_matrix_sequence
-from dune.codegen.sumfact.quadrature import (additional_inames,
- default_quadrature_inames)
+from dune.codegen.sumfact.quadrature import additional_inames
from dune.codegen.sumfact.switch import get_facedir, get_facemod
from dune.codegen.sumfact.symbolic import SumfactKernelInterfaceBase, SumfactKernel
from dune.codegen.tools import get_pymbolic_basename
@@ -46,15 +46,191 @@ from loopy.match import Writes
import pymbolic.primitives as prim
import numpy as np
+import loopy as lp
@geometry_mixin("sumfact_multilinear")
class SumfactMultiLinearGeometryMixin(GenericPDELabGeometryMixin):
- def pymbolic_unit_outer_normal(self):
- return pymbolic_unit_outer_normal(self)
+ def _jacobian_inverse(self):
+ return "jit"
- def pymbolic_jacobian_inverse(self, i, j, restriction):
- return pymbolic_jacobian_inverse(i, j, restriction, self.visitor)
+ def jacobian_inverse(self, o):
+ # This differs from the generic one by not falling back to the
+ # transpose of the jacobian inverse.
+ restriction = enforce_boundary_restriction(self)
+
+ assert(len(self.indices) == 2)
+ i, j = self.indices
+ self.indices = None
+
+ name = restricted_name(self._jacobian_inverse(), restriction)
+ self.define_jacobian_inverse(name, restriction)
+
+ return prim.Subscript(prim.Variable(name), (i,j))
+
+ def _jacobian_determinant(self):
+ return "detjac"
+
+ def define_jacobian_determinant(self, o):
+ restriction = Restriction.NONE if self.measure == "cell" else Restriction.POSITIVE
+ self.define_jacobian_inverse(self._jacobian_inverse(), restriction)
+ return prim.Variable(self._jacobian_determinant())
+
+ def define_jacobian_inverse(self, name, restriction):
+ """Return jacobian inverse 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, self)
+ inames = self.quadrature_inames() + additional_inames(self)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ depends_on=frozenset({lp.match.Tagged("sumfact_stage1")})
+ )
+ names_jacobian.append(name_jacobian)
+
+ # The sum factorization kernels from the geometry evaluation of the
+ # jacobians will never appear in the expression for the input of
+ # stage 3. This way the SumfactCollectMapper will never see them
+ # and they will be marked as inactive. Here we explicitly mark the
+ # as used.
+ basis_sf_kernels(expression.aggregate)
+
+ # Calculate the inverse of the jacobian of the geometry mapping and the
+ # determinant by calling a c++ function. Note: The result will be column
+ # major -> fortran style.
+ name_detjac = self._jacobian_determinant()
+ temporary_variable(name_detjac, shape=())
+ ftags = "f,f"
+ temporary_variable(name, shape=(dim, dim), dim_tags=ftags, managed=True)
+ include_file('dune/codegen/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 = self.quadrature_inames() + additional_inames(self)
+ return instruction(code=code,
+ assignees=frozenset([name, name_detjac]),
+ read_variables=frozenset(names_jacobian),
+ inames=frozenset(inames),
+ )
+
+ def facet_jacobian_determinant(self, o):
+ """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.
+ """
+ detjac = self.jacobian_determinant(None)
+ norm_of_outer_normal = normalize(self.outer_normal(), world_dimension())
+
+ return prim.Product((detjac, norm_of_outer_normal))
+
+ def outer_normal(self):
+ """ This is the *unnormalized* outer normal """
+ name = "outer_normal"
+ 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)
+ self.indices = (facedir_s, i)
+ ji = self.jacobian_inverse(None)
+
+ # Note: 2*facemod_s-1 because of
+ # -1 if facemod_s = 0
+ # +1 if facemod_s = 1
+ expression = ji * (2 * facemod_s - 1)
+
+ inames = self.quadrature_inames() + additional_inames(self)
+ instruction(assignee=assignee,
+ expression=expression,
+ within_inames=frozenset(inames),
+ )
+
+ return prim.Variable(name)
+
+ def facet_normal(self, o):
+ index, = self.indices
+ self.indices = None
+
+ outer_normal = self.outer_normal()
+ norm = normalize(outer_normal, world_dimension())
+
+ return prim.Subscript(outer_normal, (index,)) / norm
+
+ def spatial_coordinate(self, o):
+ """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(self.indices) == 1
+ restriction = enforce_boundary_restriction(self)
+
+ # 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(self.indices[0], restriction)
+ sf = SumfactKernel(matrix_sequence=matrix_sequence,
+ interface=inp,
+ )
+
+ from dune.codegen.sumfact.vectorization import attach_vectorization_info
+ vsf = attach_vectorization_info(sf)
+
+ # If this sum factorization kernel was not used in the dry run we
+ # just return 0
+ if vsf == 0:
+ self.indices = None
+ return 0
+
+ # Add a sum factorization kernel that implements the evaluation of
+ # the basis functions at quadrature points (stage 1)
+ from dune.codegen.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
+ self.indices = None
+
+ return prim.Subscript(var, vsf.quadrature_index(sf, self))
@geometry_mixin("sumfact_axiparallel")
@@ -108,10 +284,10 @@ class SumfactEqudistantGeometryMixin(EquidistantGeometryMixin, SumfactAxiParalle
index, = self.indices
# Urgh: *SOMEHOW* construct a face direction. This is not breaking in the unstructured
- # case, because
+ # case, because we do not enter this code path...
from dune.codegen.pdelab.restriction import Restriction
restriction = Restriction.NONE
- if get_global_context_value("integral_type") == "interior_facet":
+ if self.measure == "interior_facet":
restriction = Restriction.POSITIVE
from dune.codegen.sumfact.switch import get_facedir
face = get_facedir(restriction)
@@ -209,57 +385,6 @@ class GeoCornersInput(SumfactKernelInterfaceBase, ImmutableRecord):
return insn_dep.union(frozenset({insn}))
-@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
-
- # 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,
- )
-
- from dune.codegen.sumfact.vectorization import attach_vectorization_info
- vsf = attach_vectorization_info(sf)
-
- # If this sum factorization kernel was not used in the dry run we
- # just return 0
- if vsf == 0:
- visitor.indices = None
- return 0
-
- # Add a sum factorization kernel that implements the evaluation of
- # the basis functions at quadrature points (stage 1)
- from dune.codegen.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))
-
-
@preamble
def define_corner(name, low):
geo = name_geometry()
@@ -303,260 +428,6 @@ def name_meshwidth():
return name
-@backend(interface="spatial_coordinate", name="diagonal_transformation_matrix")
-def pymbolic_spatial_coordinate_axiparallel(do_predicates, visitor):
- assert len(visitor.indices) == 1
- index, = visitor.indices
-
- # Urgh: *SOMEHOW* construct a face direction
- from dune.codegen.pdelab.restriction import Restriction
- restriction = Restriction.NONE
- if get_global_context_value("integral_type") == "interior_facet":
- restriction = Restriction.POSITIVE
- from dune.codegen.sumfact.switch import get_facedir
- face = get_facedir(restriction)
-
- lowcorner = name_lowerleft_corner()
- meshwidth = name_meshwidth()
-
- # If we have to decide which boundary condition to take for this
- # intersection we always take the boundary condition of the center
- # of the intersection. We assume that there are no cells with more
- # than one boundary condition.
- if do_predicates:
- x = 0.5
- elif index == face:
- x = 0
- else:
- iindex = index
- if face is not None and index > face:
- iindex = iindex - 1
- from dune.codegen.sumfact.quadrature import pymbolic_indexed_quadrature_position
- x = pymbolic_indexed_quadrature_position(iindex, visitor)
-
- visitor.indices = None
- return prim.Subscript(prim.Variable(lowcorner), (index,)) + x * prim.Subscript(prim.Variable(meshwidth), (index,))
-
-
-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)
- ji = pymbolic_jacobian_inverse(facedir_s, i, Restriction.POSITIVE, visitor)
-
- # Note: 2*facemod_s-1 because of
- # -1 if facemod_s = 0
- # +1 if facemod_s = 1
- expression = ji * (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.codegen.sumfact.switch import get_facedir, get_facemod
- if index == get_facedir(Restriction.POSITIVE):
- if get_facemod(Restriction.POSITIVE):
- return 1
- else:
- return -1
- else:
- return 0
- else:
- # 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):
- # 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"):
- # 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
- else:
- return 1
- else:
- return 0
- else:
- # 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=-1 * prim.Subscript(outer_normal, i) / norm,
- within_inames=frozenset(inames),
- )
-
- return prim.Variable(name)
-
-
-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):
- name = name_facet_jacobian_determinant(visitor)
- return prim.Variable(name)
-
-
-def pymbolic_facet_area(visitor):
- from dune.codegen.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
@@ -599,79 +470,5 @@ def _name_jacobian(i, j, restriction, visitor):
return prim.Subscript(var, vsf.quadrature_index(sf, visitor))
-def define_jacobian_inverse(name, restriction, visitor):
- """Return jacobian inverse 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)
-
- # The sum factorization kernels from the geometry evaluation of the
- # jacobians will never appear in the expression for the input of
- # stage 3. This way the SumfactCollectMapper will never see them
- # and they will be marked as inactive. Here we explicitly mark the
- # as used.
- basis_sf_kernels(expression.aggregate)
-
- # Calculate the inverse of the jacobian of the geometry mapping and the
- # determinant by calling a c++ function. Note: The result will be column
- # major -> fortran style.
- name_detjac = name_jacobian_determinant(visitor)
- temporary_variable(name_detjac, shape=())
- ftags = ",".join(["f"] * 2)
- temporary_variable(name, shape=(dim, dim), dim_tags=ftags, managed=True)
- include_file('dune/codegen/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):
- 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):
- # 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))
+def normalize(expr, dim):
+ return prim.Call(prim.Variable("sqrt"), (prim.Sum(tuple(expr[i]*expr[i] for i in range(dim))),))
diff --git a/python/dune/codegen/sumfact/quadrature.py b/python/dune/codegen/sumfact/quadrature.py
index 92a598d9..86982443 100644
--- a/python/dune/codegen/sumfact/quadrature.py
+++ b/python/dune/codegen/sumfact/quadrature.py
@@ -141,15 +141,6 @@ def base_weight_function_mangler(target, func, dtypes):
return CallMangleInfo(func.name, (NumpyType(dtype_floatingpoint()),), ())
-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:
@@ -193,8 +184,7 @@ def additional_inames(visitor):
restriction = Restriction.NONE
if get_global_context_value("integral_type") == "interior_facet":
restriction = Restriction.POSITIVE
- from dune.codegen.sumfact.basis import lfs_inames
- lfs_inames = lfs_inames(element, restriction)
+ lfs_inames = visitor.lfs_inames(element, restriction)
return lfs_inames
else:
return ()
diff --git a/python/dune/codegen/ufl/visitor.py b/python/dune/codegen/ufl/visitor.py
index 10afd7bf..376fe9d5 100644
--- a/python/dune/codegen/ufl/visitor.py
+++ b/python/dune/codegen/ufl/visitor.py
@@ -56,14 +56,14 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
return self._call(o, do_predicates)
def accumulate(self, o):
- for info in self.interface.list_accumulation_infos(o, self):
+ for info in self.list_accumulation_infos(o):
self.current_info = info
expr = self._call(o, False)
if expr != 0:
if get_form_option("simplify"):
from dune.codegen.sympy import simplify_pymbolic_expression
expr = simplify_pymbolic_expression(expr)
- self.interface.generate_accumulation_instruction(expr, self)
+ self.generate_accumulation_instruction(expr)
def _call(self, o, do_predicates):
# Reset state variables
@@ -96,7 +96,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker):
def argument(self, o):
self.interface.initialize_function_spaces(o, self)
# Update the information on where to accumulate this
- info = self.interface.get_accumulation_info(o, self)
+ info = self.get_accumulation_info(o)
if o.number() == 0:
if info != self.current_info[0]:
self.indices = None
diff --git a/test/sumfact/poisson/poisson_2d_unstructured.mini b/test/sumfact/poisson/poisson_2d_unstructured.mini
index b7de167f..5eabe782 100644
--- a/test/sumfact/poisson/poisson_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_2d_unstructured.mini
@@ -30,6 +30,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_3d_unstructured.mini b/test/sumfact/poisson/poisson_3d_unstructured.mini
index 0d8b4690..0cbf2814 100644
--- a/test/sumfact/poisson/poisson_3d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_3d_unstructured.mini
@@ -30,6 +30,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d_gmsh.mini b/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
index 79df9ea4..054a2481 100644
--- a/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_dg_2d_gmsh.mini
@@ -29,6 +29,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
index fec1c93d..55336ae7 100644
--- a/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_dg_2d_unstructured.mini
@@ -30,6 +30,7 @@ numerical_jacobian = 1, 0 | expand num
sumfact = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d_gmsh.mini b/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
index eb320a2d..0216093f 100644
--- a/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_dg_3d_gmsh.mini
@@ -31,6 +31,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_dg_3d_unstructured.mini b/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
index 271b1f54..c42be3b9 100644
--- a/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
+++ b/test/sumfact/poisson/poisson_dg_3d_unstructured.mini
@@ -31,6 +31,7 @@ sumfact = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
diff --git a/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini b/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
index 52de26ab..a6374d18 100644
--- a/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
+++ b/test/sumfact/poisson/poisson_fastdg_3d_gmsh.mini
@@ -30,6 +30,7 @@ fastdg = 1
sumfact_regular_jacobians = 1
vectorization_quadloop = 1, 0 | expand quad
vectorization_strategy = explicit, none | expand grad
+geometry_mixins = sumfact_multilinear
[formcompiler.ufl_variants]
degree = 1, 2 | expand deg
--
GitLab