diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py index 2f5455230fc1273bad61c380918c1a7778f91342..4b32d5d9b73e7a3b886f4868e5755833308cd1e7 100644 --- a/python/dune/perftool/blockstructured/basis.py +++ b/python/dune/perftool/blockstructured/basis.py @@ -2,16 +2,17 @@ from dune.perftool.generation import (backend, kernel_cached, get_backend, instruction, - temporary_variable) + temporary_variable, + globalarg, + class_member, + initializer_list) from dune.perftool.tools import get_pymbolic_basename -from dune.perftool.pdelab.driver import FEM_name_mangling -from dune.perftool.pdelab.restriction import restricted_name from dune.perftool.pdelab.basis import (declare_cache_temporary, - name_localbasis_cache) + name_localbasis_cache, + type_localbasis + ) from dune.perftool.pdelab.geometry import world_dimension from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell -from dune.perftool.blockstructured.spaces import lfs_inames -import pymbolic.primitives as prim @backend(interface="evaluate_basis", name="blockstructured") @@ -20,11 +21,9 @@ def evaluate_basis(leaf_element, name, restriction): temporary_variable(name, shape=(4,), decl_method=declare_cache_temporary(leaf_element, restriction, 'Function')) cache = name_localbasis_cache(leaf_element) qp = pymbolic_quadrature_position_in_cell(restriction) + localbasis = name_localbasis(leaf_element) instruction(inames=get_backend("quad_inames")(), - code='{} = {}.evaluateFunction({}, lfs.finiteElement().localBasis());'.format(name, - cache, - str(qp), - ), + code='{} = {}.evaluateFunction({}, {});'.format(name, cache, str(qp), localbasis), assignees=frozenset({name}), read_variables=frozenset({get_pymbolic_basename(qp)}), ) @@ -33,14 +32,26 @@ def evaluate_basis(leaf_element, name, restriction): @backend(interface="evaluate_grad", name="blockstructured") @kernel_cached def evaluate_reference_gradient(leaf_element, name, restriction): - temporary_variable(name, shape=(4,1,world_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian')) + temporary_variable(name, shape=(4, 1, world_dimension()), decl_method=declare_cache_temporary(leaf_element, restriction, 'Jacobian')) cache = name_localbasis_cache(leaf_element) qp = pymbolic_quadrature_position_in_cell(restriction) + localbasis = name_localbasis(leaf_element) instruction(inames=get_backend("quad_inames")(), - code='{} = {}.evaluateJacobian({}, lfs.finiteElement().localBasis());'.format(name, - cache, - str(qp), - ), + code='{} = {}.evaluateJacobian({}, {});'.format(name, cache, str(qp), localbasis), assignees=frozenset({name}), read_variables=frozenset({get_pymbolic_basename(qp)}), - ) \ No newline at end of file + ) + + +@class_member(classtag="operator") +def define_localbasis(leaf_element, name): + localBasis_type = type_localbasis(leaf_element) + initializer_list(name, (), classtag="operator") + return "const {} {};".format(localBasis_type, name) + + +def name_localbasis(leaf_element): + name = "microElementBasis" + globalarg(name) + define_localbasis(leaf_element, name) + return name diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py index 2195e21f70338be438b71c5ff913e1760bfd388b..fd51a2e619b51af96d240ac2cb6cd869deb7b73d 100644 --- a/python/dune/perftool/pdelab/basis.py +++ b/python/dune/perftool/pdelab/basis.py @@ -32,18 +32,32 @@ import pymbolic.primitives as prim from loopy import Reduction -def type_localbasis(element): +@class_member(classtag="operator") +def typedef_localbasis(element, name): df = "typename {}::Traits::GridView::ctype".format(type_gfs(element)) r = basetype_range() dim = world_dimension() # if isPk(element): # include_file("dune/localfunctions/lagrange/pk/pklocalbasis.hh", filetag="operatorfile") # return "Dune::PkLocalBasis<{}, {}, {}, {}>".format(df, r, dim, element._degree) + #TODO add dg support if isQk(element): include_file("dune/localfunctions/lagrange/qk/qklocalbasis.hh", filetag="operatorfile") - return "Dune::QkLocalBasis<{}, {}, {}, {}>".format(df, r, dim, element._degree) - #TODO add dg support - raise NotImplementedError("Element type not known in code generation") + basis_type = "QkLocalBasis<{}, {}, {}, {}>".format(df, r, element._degree, dim) + else: + raise NotImplementedError("Element type not known in code generation") + return "using {} = Dune::{};".format(name, basis_type) + + +def type_localbasis(element): + if isPk(element): + name = "P{}_LocalBasis".format(element._degree) + elif isQk(element): + name = "Q{}_LocalBasis".format(element._degree) + else: + raise NotImplementedError("Element type not known in code generation") + typedef_localbasis(element, name) + return name def type_localbasis_cache(element): diff --git a/python/dune/perftool/ufl/transformations/blockstructured.py b/python/dune/perftool/ufl/transformations/blockstructured.py deleted file mode 100644 index d25c95e10fd958c30f7984e41b62aa25b17fb290..0000000000000000000000000000000000000000 --- a/python/dune/perftool/ufl/transformations/blockstructured.py +++ /dev/null @@ -1,39 +0,0 @@ -from dune.perftool.options import get_option -from dune.perftool.ufl.transformations import ufl_transformation -from dune.perftool.ufl.transformations.replace import ReplaceExpression -from ufl.algorithms import MultiFunction -from ufl import as_ufl -from ufl.classes import JacobianInverse, JacobianDeterminant, Product, Division, Indexed - - -class ReplaceReferenceTransformation(MultiFunction): - def __init__(self, k): - MultiFunction.__init__(self) - self.k = k - self.visited_jit = False - - def expr(self, o): - return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands)) - -#TODO abs uses c abs -> only works for ints!!! - def abs(self,o): - if isinstance(o.ufl_operands[0], JacobianDeterminant): - return Division(o, as_ufl(self.k**2)) - else: - return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands)) - - def jacobian_determinant(self,o): - return Division(o, as_ufl(self.k**2)) - - def indexed(self, o): - expr = o.ufl_operands[0] - multiindex = o.ufl_operands[1] - if isinstance(expr, JacobianInverse): - return Product(as_ufl(self.k), Indexed(expr, multiindex)) - else: - return self.reuse_if_untouched(o, *tuple(self(op) for op in o.ufl_operands)) - - -@ufl_transformation(name="blockstructured") -def blockstructured(expr): - return ReplaceReferenceTransformation(get_option("number_of_blocks"))(expr) \ No newline at end of file