diff --git a/python/dune/codegen/blockstructured/basis.py b/python/dune/codegen/blockstructured/basis.py
index e5aa4d20ebef0eb320fff559c2d2fb8b54fe909b..69097e24dc9763f9b9b7d5cb285bfcda792f1bbc 100644
--- a/python/dune/codegen/blockstructured/basis.py
+++ b/python/dune/codegen/blockstructured/basis.py
@@ -1,3 +1,5 @@
+from loopy import Reduction
+
from dune.codegen.generation import (backend,
basis_mixin,
kernel_cached,
@@ -8,7 +10,8 @@ from dune.codegen.generation import (backend,
class_member,
initializer_list,
include_file,)
-from dune.codegen.tools import get_pymbolic_basename
+from dune.codegen.pdelab.argument import name_coefficientcontainer
+from dune.codegen.tools import get_pymbolic_basename, get_pymbolic_indices
from dune.codegen.loopy.target import type_floatingpoint
from dune.codegen.pdelab.basis import (GenericBasisMixin,
declare_cache_temporary,
@@ -19,7 +22,7 @@ from dune.codegen.pdelab.driver import (isPk,
isQk,
)
from dune.codegen.pdelab.geometry import world_dimension
-from dune.codegen.pdelab.spaces import type_leaf_gfs
+from dune.codegen.pdelab.spaces import type_leaf_gfs, name_lfs
from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.blockstructured.spaces import lfs_inames
from dune.codegen.blockstructured.tools import tensor_index_to_sequential_index
@@ -48,7 +51,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1),
decl_method=declare_cache_temporary(element, restriction, 'Function'))
cache = name_localbasis_cache(element)
- qp = self.to_cell(self.quadrature_position())
+ qp = self.to_cell(self.quadrature_position_in_micro())
localbasis = name_localbasis(element)
instruction(inames=self.quadrature_inames(),
code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis),
@@ -70,7 +73,7 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
temporary_variable(name, shape=((element.degree() + 1)**world_dimension(), 1, world_dimension()),
decl_method=declare_cache_temporary(element, restriction, 'Jacobian'))
cache = name_localbasis_cache(element)
- qp = self.to_cell(self.quadrature_position())
+ qp = self.to_cell(self.quadrature_position_in_micro())
localbasis = name_localbasis(element)
instruction(inames=self.quadrature_inames(),
code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis),
@@ -78,6 +81,41 @@ class BlockStructuredBasisMixin(GenericBasisMixin):
read_variables=frozenset({get_pymbolic_basename(qp)}),
)
+ def implement_trialfunction_gradient(self, element, restriction, index):
+ rawname = "gradu_{}".format(index)
+ name = restricted_name(rawname, restriction)
+ container = name_coefficientcontainer(restriction)
+ self.evaluate_coefficient_gradient(element, name, container, restriction, index)
+ return prim.Variable(name)
+
+ @kernel_cached
+ def evaluate_coefficient(self, element, name, container, restriction, index):
+ sub_element = element
+ if isinstance(element, MixedElement):
+ sub_element = element.extract_component(index)[1]
+
+ from ufl import FiniteElement
+ assert isinstance(sub_element, FiniteElement)
+
+ temporary_variable(name, shape=(), managed=True)
+
+ lfs = name_lfs(element, restriction, index)
+ basis = self.implement_basis(sub_element, restriction, 0, context='trial')
+ basisindex = get_pymbolic_indices(basis)[:-1]
+
+ # TODO get rid ot this!
+ from dune.codegen.blockstructured.argument import pymbolic_coefficient
+ coeff = pymbolic_coefficient(container, lfs, sub_element, basisindex)
+
+ assignee = prim.Variable(name)
+ reduction_expr = prim.Product((coeff, basis))
+
+ instruction(expression=Reduction("sum", basisindex, reduction_expr, allow_simultaneous=True),
+ assignee=assignee,
+ forced_iname_deps=frozenset(self.quadrature_inames()),
+ forced_iname_deps_is_final=True,
+ )
+
# define FE basis explicitly in localoperator
@backend(interface="typedef_localbasis", name="blockstructured")
diff --git a/python/dune/codegen/blockstructured/geometry.py b/python/dune/codegen/blockstructured/geometry.py
index a2bfd1de9bed8363e693b44e4270f4de5916db90..5bb48aebc4c692617fbef78af29e6baab6291d03 100644
--- a/python/dune/codegen/blockstructured/geometry.py
+++ b/python/dune/codegen/blockstructured/geometry.py
@@ -34,10 +34,16 @@ class AxiparallelBlockStructuredGeometryMixin(AxiparallelGeometryMixin, BlockStr
return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
prim.Power(get_form_option("number_of_blocks"), local_dimension()))
- def jacobian_inverse(self, o):
- return prim.Quotient(AxiparallelGeometryMixin.facet_jacobian_determinant(self, o),
+ def jacobian_determinant(self, o):
+ name = 'detjac'
+ self.define_jacobian_determinant(name)
+ return prim.Quotient(prim.Variable(name),
prim.Power(get_form_option("number_of_blocks"), local_dimension()))
+ def jacobian_inverse(self, o):
+ return prim.Product((AxiparallelGeometryMixin.jacobian_inverse(self, o),
+ get_form_option("number_of_blocks")))
+
@geometry_mixin("blockstructured_equidistant")
class EquidistantBlockStructuredGeometryMixin(EquidistantGeometryMixin, AxiparallelBlockStructuredGeometryMixin):
diff --git a/python/dune/codegen/blockstructured/quadrature.py b/python/dune/codegen/blockstructured/quadrature.py
index f69ddae816d876fd69a45a4e80844d58bee43bf6..ac45eda11a113c9b9d5642d2e09915334c9d9b78 100644
--- a/python/dune/codegen/blockstructured/quadrature.py
+++ b/python/dune/codegen/blockstructured/quadrature.py
@@ -8,8 +8,10 @@ import pymbolic.primitives as prim
@quadrature_mixin("blockstructured")
class BlockstructuredQuadratureMixin(GenericQuadratureMixin):
- def quadrature_position(self, index=None):
- original = GenericQuadratureMixin.quadrature_position(self)
+ quadrature_position_in_micro = GenericQuadratureMixin.quadrature_position
+
+ def quadrature_position_in_macro(self, index=None):
+ original = quadrature_position_in_micro(self, index)
qp = prim.Variable(name_point_in_macro(original, self), )
if index is not None:
return prim.Subscript(qp, (index,))
diff --git a/python/dune/codegen/options.py b/python/dune/codegen/options.py
index 79ee0c051b3a4fe2d667c00c61f84c7a85bd702c..28962be77db2cbe3a181a6ae771367877acaae9c 100644
--- a/python/dune/codegen/options.py
+++ b/python/dune/codegen/options.py
@@ -184,6 +184,7 @@ def process_form_options(opt, form):
opt = opt.copy(accumulation_mixins="blockstructured",
geometry_mixins="blockstructured",
quadrature_mixins="blockstructured",
+ basis_mixins="blockstructured"
)
if opt.control: