Skip to content
Snippets Groups Projects
Commit e27992aa authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Allow multiple measures of same type -> mixed bc poisson dg

parent bea763e5
No related branches found
No related tags found
No related merge requests found
......@@ -154,20 +154,21 @@ def assembly_routine_signature():
assert False
def generate_kernel(integral):
integrand = integral.integrand()
measure = integral.integral_type()
subdomain_id = integral.subdomain_id()
subdomain_data = integral.subdomain_data()
# Now split the given integrand into accumulation expressions
from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
accterms = split_into_accumulation_terms(integrand)
# Iterate over the terms and generate a kernel
for term in accterms:
from dune.perftool.loopy.transformer import transform_accumulation_term
transform_accumulation_term(term, measure, subdomain_id)
def generate_kernel(integrals):
for integral in integrals:
integrand = integral.integrand()
measure = integral.integral_type()
subdomain_id = integral.subdomain_id()
subdomain_data = integral.subdomain_data()
# Now split the given integrand into accumulation expressions
from dune.perftool.ufl.transformations.extract_accumulation_terms import split_into_accumulation_terms
accterms = split_into_accumulation_terms(integrand)
# Iterate over the terms and generate a kernel
for term in accterms:
from dune.perftool.loopy.transformer import transform_accumulation_term
transform_accumulation_term(term, measure, subdomain_id)
# Extract the information, which is needed to create a loopy kernel.
# First extracting it, might be useful to alter it before kernel generation.
......@@ -247,10 +248,6 @@ def generate_localoperator_kernels(formdata, namedata):
# Extract the relevant attributes of the form data
form = formdata.preprocessed_form
# For the moment, I do assume that there is but one integral of each type. This might differ
# if you use different quadrature orders for different terms.
assert len(form.integrals()) == len(set(i.integral_type() for i in form.integrals()))
# Reset the generation cache
from dune.perftool.generation import delete_cache_items
delete_cache_items()
......@@ -283,23 +280,23 @@ def generate_localoperator_kernels(formdata, namedata):
with global_context(formdata=formdata, namedata=namedata):
with global_context(form_type='residual'):
# Generate the necessary residual methods
for integral in form.integrals():
# Reset the outer loop
from dune.perftool.loopy.transformer import set_outerloop
set_outerloop(None)
for measure in set(i.integral_type() for i in form.integrals()):
with global_context(integral_type=measure):
# Reset the outer loop
from dune.perftool.loopy.transformer import set_outerloop
set_outerloop(None)
with global_context(integral_type=integral.integral_type()):
enum_pattern()
pattern_baseclass()
enum_alpha()
kernel = generate_kernel(integral)
kernel = generate_kernel(form.integrals_by_type(measure))
# Maybe add numerical differentiation
if get_option("numerical_jacobian"):
include_file("dune/pdelab/localoperator/defaultimp.hh", filetag="operatorfile")
_, loptype = class_type_from_cache("operator")
which = ufl_measure_to_pdelab_measure(integral.integral_type())
which = ufl_measure_to_pdelab_measure(measure)
base_class("Dune::PDELab::NumericalJacobian{}<{}>".format(which, loptype), classtag="operator")
# Add the initializer list for that base class
......@@ -310,7 +307,7 @@ def generate_localoperator_kernels(formdata, namedata):
classtag="operator",
)
operator_kernels[(integral.integral_type(), 'residual')] = kernel
operator_kernels[(measure, 'residual')] = kernel
# Generate the necessary jacobian methods
if not get_option("numerical_jacobian"):
......@@ -319,14 +316,14 @@ def generate_localoperator_kernels(formdata, namedata):
jacform = expand_derivatives(derivative(form, form.coefficients()[0]))
with global_context(form_type="jacobian"):
for integral in jacform.integrals():
for measure in set(i.integral_type() for i in jacform.integrals()):
# Reset the outer loop
from dune.perftool.loopy.transformer import set_outerloop
set_outerloop(None)
with global_context(integral_type=integral.integral_type()):
kernel = generate_kernel(integral)
operator_kernels[(integral.integral_type(), 'jacobian')] = kernel
with global_context(integral_type=measure):
kernel = generate_kernel(jacform.integrals_by_type(measure))
operator_kernels[(measure, 'jacobian')] = kernel
# Generate dummy functions for those kernels, that vanished in the differentiation process
# We *could* solve this problem by using lambda_* terms but we do not really want that, so
......
......@@ -57,6 +57,11 @@ class Expression(Coefficient):
# Initialize a coefficient with a dummy finite element map.
Coefficient.__init__(self, _dg0)
# TODO the subdomain_data code needs a uflid, not idea how to get it here
# The standard way through class decorator fails here...
def ufl_id(self):
return 0
class FiniteElement(ufl.FiniteElement):
def __init__(self, *args, **kwargs):
......
......@@ -77,6 +77,15 @@ dune_add_system_test(TARGET poisson_dg_neumann_numdiff
INIFILE poisson_dg_neumann_numdiff.mini
)
# 8. Poisson Test Case: DG, mixed bc, symbolic differentiation
add_generated_executable(UFLFILE poisson_dg_neumann.ufl
TARGET poisson_dg_neumann_symdiff
)
dune_add_system_test(TARGET poisson_dg_neumann_symdiff
INIFILE poisson_dg_neumann_symdiff.mini
)
# Add the reference code with the PDELab localoperator that produced
# the reference vtk file
add_executable(poisson_dg_ref reference_main.cc)
......
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