Skip to content
Snippets Groups Projects
Commit 0cebef06 authored by Marcel Koch's avatar Marcel Koch Committed by Dominic Kempf
Browse files

small fixes

parent fd3fe07c
No related branches found
No related tags found
No related merge requests found
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")
......
......@@ -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):
......
......@@ -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,))
......
......@@ -188,6 +188,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:
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment