diff --git a/python/dune/codegen/pdelab/driver/constraints.py b/python/dune/codegen/pdelab/driver/constraints.py
index aed658caaa1fef5bef7bcb0615f6383efcdf8d3a..344294a183006b072016a8b414de90742ca8d4cd 100644
--- a/python/dune/codegen/pdelab/driver/constraints.py
+++ b/python/dune/codegen/pdelab/driver/constraints.py
@@ -15,6 +15,7 @@ from dune.codegen.pdelab.driver.gridfunctionspace import (name_gfs,
type_trial_gfs,
preprocess_leaf_data,
)
+from dune.codegen.pdelab.driver.interpolate import name_grid_function_root
from ufl.classes import Expr
from ufl import FiniteElement, MixedElement, TensorElement, VectorElement, TensorProductElement
@@ -43,7 +44,7 @@ def assemble_constraints(name):
gfs = name_trial_gfs()
is_dirichlet = preprocess_leaf_data(element, "is_dirichlet")
if has_dirichlet_constraints(is_dirichlet):
- bctype_function = name_bctype_grid_function(element, is_dirichlet)
+ bctype_function = name_grid_function_root("is_dirichlet")
return "Dune::PDELab::constraints(*{}, *{}, *{});".format(bctype_function,
gfs,
name,)
@@ -53,118 +54,7 @@ def assemble_constraints(name):
#
-# Bctype lambda
-#
-
-
-def bctype_lambda(func, local=False):
- from ufl.classes import Expr
- if func is None:
- func = 0.0
- if isinstance(func, (int, float)):
- return "[&](const auto& is, const auto& xl){{ return {}; }}".format(float(func))
- elif isinstance(func, Expr):
- from dune.codegen.pdelab.driver.visitor import ufl_to_code
- return "[&](const auto& is, const auto& xl){{ {}; }}".format(ufl_to_code(func))
- raise ValueError("Expression not understood")
-
-
-#
-# Bctype function
-#
-# std::function that will be used to create the bctype grid function
-
-
-@class_member(classtag="driver_block")
-def typedef_bctype_function(name):
- leafview_type = type_leafview()
- entity = "typename {}::Intersection".format(leafview_type)
- coordinate = "typename {}::Intersection::LocalCoordinate".format(leafview_type)
- return "using {} = std::function<bool({}, {})>;".format(name, entity, coordinate)
-
-
-def type_bctype_function():
- name = "BctypeFunction"
- typedef_bctype_function(name)
- return name
-
-
-@class_member(classtag="driver_block")
-def declare_bctype_function(name, element, is_dirichlet):
- bctype_function_type = type_bctype_function()
- return "std::shared_ptr<{}> {};".format(bctype_function_type, name)
-
-
-@preamble(section="constriants", kernel="driver_block")
-def define_bctype_function(name, element, is_dirichlet):
- declare_bctype_function(name, element, is_dirichlet)
- bctype_function_type = type_bctype_function()
- bct_lambda = bctype_lambda(is_dirichlet)
- return "{} = std::make_shared<{}> ({});".format(name, bctype_function_type, bct_lambda)
-
-
-def name_bctype_function(element, is_dirichlet):
- name = get_counted_variable("bctype_function")
- define_bctype_function(name, element, is_dirichlet)
- return name
-
-
-#
-# Bctype Grid Function
-#
-# GridFunction that returns 1 if the boundary is constraint (else 0)
-
-
-@class_member(classtag="driver_block")
-def typedef_bctype_grid_function(name):
- bctype_function_type = type_bctype_function()
- return "using {} = Dune::PDELab::LocalCallableToBoundaryConditionAdapter<{}>;".format(name,
- bctype_function_type)
-
-
-def type_bctype_grid_function():
- name = "BctypeGridFunction"
- typedef_bctype_grid_function(name)
- return name
-
-
-@class_member(classtag="driver_block")
-def declare_bctype_grid_function(name):
- bctype_type = type_bctype_grid_function()
- return "std::shared_ptr<{}> {};".format(bctype_type, name)
-
-
-@preamble(section="constraints", kernel="driver_block")
-def define_bctype_grid_function(name, element, is_dirichlet):
- declare_bctype_grid_function(name)
- bctype_type = type_bctype_grid_function()
- bctype_function = name_bctype_function(element, is_dirichlet)
- return "{} = std::make_shared<{}> (*{});".format(name, bctype_type, bctype_function)
-
-
-def name_bctype_grid_function(element, is_dirichlet):
- # Note: Previously, there was a separate code branch for VectorElement here,
- # which was implemented through PDELabs Power constraints concept.
- # However this completely fails if you have different constraints for
- # the children of a VectorElement. We therefore omit this branch completely.
- if isinstance(element, MixedElement):
- k = 0
- childs = []
- for subel in element.sub_elements():
- childs.append(name_bctype_grid_function(subel, is_dirichlet[k:k + subel.value_size()]))
- k = k + subel.value_size()
- name = "{}".format("_".join(childs))
- define_composite_constraints_parameters(name, element, tuple(childs))
- return name
- else:
- assert isinstance(element, (FiniteElement, TensorProductElement))
- name = get_counted_variable("bctype_grid_function")
- define_bctype_grid_function(name, element, is_dirichlet[0])
- return name
-
-
-#
-# Composite constraints parameter
+# Composite constraints parameters
#
diff --git a/python/dune/codegen/pdelab/driver/gridoperator.py b/python/dune/codegen/pdelab/driver/gridoperator.py
index 8d0001cab2c7451d68db31332b7ef200f236f0d5..cf8504c2652abd3a1104815de46257a7746c7f3a 100644
--- a/python/dune/codegen/pdelab/driver/gridoperator.py
+++ b/python/dune/codegen/pdelab/driver/gridoperator.py
@@ -101,6 +101,11 @@ def driver_block_get_gridoperator(form_ident, name=None):
gridoperator_type = type_gridoperator(form_ident)
if not name:
name = name_gridoperator(form_ident)
+ if name == "go_mass":
+ return ["std::shared_ptr<{}> getMassGridOperator(){{".format(gridoperator_type),
+ " return {};".format(name),
+ "}"]
+
return ["std::shared_ptr<{}> getGridOperator(){{".format(gridoperator_type),
" return {};".format(name),
"}"]
@@ -154,6 +159,36 @@ def define_localoperator(name, form_ident):
def name_localoperator(form_ident):
name = "lop_{}".format(form_ident)
define_localoperator(name, form_ident)
+ driver_block_get_localoperator(form_ident, name=name)
+ return name
+
+
+@class_member(classtag="driver_block")
+def driver_block_get_localoperator(form_ident, name=None):
+ if not name:
+ name = name_localoperator(form_ident)
+ localoperator_type = type_localoperator(form_ident)
+ method_name = "getLocalOperator"
+ if form_ident == "mass":
+ method_name = "getLocalMassOperator"
+ return ["std::shared_ptr<{}> {}(){{".format(localoperator_type, method_name),
+ " return {};".format(name),
+ "}"]
+
+
+@preamble(section="postprocessing", kernel="main")
+def main_define_localoperator(name, form_ident):
+ driver_block_name = name_driver_block()
+ driver_block_get_localoperator(form_ident)
+ method_name = "getLocalOperator"
+ if form_ident == "mass":
+ method_name = "getLocalMassOperator"
+ return "auto {} = {}.{}();".format(name, driver_block_name, method_name)
+
+
+def main_name_localoperator(form_ident):
+ name = "localOperator{}".format(form_ident.capitalize())
+ main_define_localoperator(name, form_ident)
return name
diff --git a/python/dune/codegen/pdelab/driver/instationary.py b/python/dune/codegen/pdelab/driver/instationary.py
index c8edb4747fedee22a4c27612b5794aa29c32cc0f..51d20766c4ffb0a3dbe3ab3f3328aa2bce8d9de2 100644
--- a/python/dune/codegen/pdelab/driver/instationary.py
+++ b/python/dune/codegen/pdelab/driver/instationary.py
@@ -1,4 +1,5 @@
-from dune.codegen.generation import (include_file,
+from dune.codegen.generation import (class_member,
+ include_file,
preamble,
)
from dune.codegen.pdelab.driver import (get_form_ident,
@@ -12,7 +13,7 @@ from dune.codegen.pdelab.driver.gridfunctionspace import (name_trial_gfs,
)
from dune.codegen.pdelab.driver.gridoperator import (name_gridoperator,
type_gridoperator,
- name_localoperator,
+ main_name_localoperator,
)
from dune.codegen.pdelab.driver.constraints import (has_dirichlet_constraints,
name_bctype_grid_function,
@@ -23,6 +24,8 @@ from dune.codegen.pdelab.driver.interpolate import (interpolate_dirichlet_data,
)
from dune.codegen.pdelab.driver.solve import (print_matrix,
print_residual,
+ main_name_vector,
+ main_type_vector,
name_linearsolver,
name_stationarynonlinearproblemsolver,
name_vector,
@@ -44,18 +47,17 @@ def solve_instationary():
raise NotImplementedError("Instationary matrix free not implemented!")
else:
time_loop()
-
print_residual()
print_matrix()
-@preamble(section="instat")
+@preamble(section="instat", kernel="main")
def time_loop():
ini = name_initree()
- lop = name_localoperator(get_form_ident())
+ lop = main_name_localoperator(get_form_ident())
time = name_time()
element = get_trial_element()
- vector_type = type_vector(get_form_ident())
+ vector_type = main_type_vector(get_form_ident())
vector = name_vector(get_form_ident())
interpolate_dirichlet_data(vector)
gfs = name_trial_gfs()
@@ -66,7 +68,7 @@ def time_loop():
bctype = name_bctype_grid_function(element, is_dirichlet)
cc = name_constraintscontainer()
assemble_new_constraints = (" // Assemble constraints for new time step\n"
- " {}.setTime({}+dt);\n"
+ " {}->setTime({}+dt);\n"
" Dune::PDELab::constraints({}, {}, {});\n"
"\n".format(lop, time, bctype, gfs, cc)
)
@@ -79,7 +81,7 @@ def time_loop():
else:
osm = name_onestepmethod()
if has_dirichlet_constraints(is_dirichlet):
- boundary = name_boundary_grid_function_root("interpolate_expression")
+ boundary = name_grid_function_root("interpolate_expression")
apply_call = "{}.apply(time, dt, {}, {}, {}new);".format(osm, vector, boundary, vector)
else:
apply_call = "{}.apply(time, dt, {}, {}new);".format(osm, vector, vector)
@@ -117,7 +119,7 @@ def time_loop():
""]
-@preamble(section="init")
+@preamble(section="init", kernel="main")
def define_time(name):
return "double {} = 0.0;".format(name)
@@ -157,19 +159,26 @@ def type_timesteppingmethod():
return "TSM"
-@preamble(section="instat")
+@class_member(classtag="driver_block")
+def declare_timesteppingmethod(name):
+ tsm_type = type_timesteppingmethod()
+ return "std::shared_ptr<{}> {};".format(tsm_type, name)
+
+
+@preamble(section="instat", kernel="driver_block")
def define_timesteppingmethod(name):
+ declare_timesteppingmethod(name)
tsm_type = type_timesteppingmethod()
explicit = get_option('explicit_time_stepping')
if explicit:
- return "{} {};".format(tsm_type, name)
+ return "{} = std::make_shared<{}>();".format(name, tsm_type)
else:
order = get_option('time_stepping_order')
if order == 1:
ini = name_initree()
- return "{} {}({}.get<double>(\"instat.theta\",1.0));".format(tsm_type, name, ini)
+ return "{} = std::make_shared<{}>({}.get<double>(\"instat.theta\",1.0));".format(name, tsm_type, ini)
else:
- return "{} {};".format(tsm_type, name)
+ return "{} = std::make_shared<{}>();".format(name, tsm_type)
def name_timesteppingmethod():
@@ -194,12 +203,19 @@ def type_instationarygridoperator():
return "IGO"
-@preamble(section="gridoperator")
+@class_member(classtag="driver_block")
+def declare_instationarygridoperator(name):
+ igo_type = type_instationarygridoperator()
+ return "std::shared_ptr<{}> {};".format(igo_type, name)
+
+
+@preamble(section="gridoperator", kernel="driver_block")
def define_instationarygridoperator(name):
+ declare_instationarygridoperator(name)
igo_type = type_instationarygridoperator()
go = name_gridoperator(get_form_ident())
mass_go = name_gridoperator("mass")
- return "{} {}({}, {});".format(igo_type, name, go, mass_go)
+ return "{} = std::make_shared<{}>({}, {});".format(name, igo_type, go, mass_go)
def name_instationarygridoperator():
@@ -221,14 +237,21 @@ def type_onestepmethod():
return "OSM"
-@preamble(section="instat")
+@class_member(classtag="driver_block")
+def declare_onestepmethod(name):
+ ilptype = type_onestepmethod()
+ return "std::shared_ptr<{}> {};".format(ilptype, name)
+
+
+@preamble(section="instat", kernel="driver_block")
def define_onestepmethod(name):
+ declare_onestepmethod(name)
ilptype = type_onestepmethod()
tsm = name_timesteppingmethod()
igo_type = type_instationarygridoperator()
igo = name_instationarygridoperator()
snp = name_stationarynonlinearproblemsolver(igo_type, igo)
- return "{} {}({},{},{});".format(ilptype, name, tsm, igo, snp)
+ return "{} = std::make_shared<{}>({},{},{});".format(name, ilptype, tsm, igo, snp)
def name_onestepmethod():
@@ -250,13 +273,20 @@ def type_explicitonestepmethod():
return "EOSM"
-@preamble(section="instat")
+@class_member(classtag="driver_block")
+def declare_explicitonestepmethod(name):
+ eosm_type = type_explicitonestepmethod()
+ return "std::shared_ptr<{}> {};".format(eosm_type, name)
+
+
+@preamble(section="instat", kernel="driver_block")
def define_explicitonestepmethod(name):
+ declare_explicitonestepmethod(name)
eosm_type = type_explicitonestepmethod()
tsm = name_timesteppingmethod()
igo = name_instationarygridoperator()
ls = name_linearsolver()
- return "{} {}({}, {}, {});".format(eosm_type, name, tsm, igo, ls)
+ return "{} = std::make_shared<{}>({}, {}, {});".format(name, eosm_type, tsm, igo, ls)
def name_explicitonestepmethod():
diff --git a/python/dune/codegen/pdelab/driver/interpolate.py b/python/dune/codegen/pdelab/driver/interpolate.py
index 6a52f485efd93ee43c217ecddd46104b18d8840d..15efff4ba60040c0b972d1e7a4bd79794300ebec 100644
--- a/python/dune/codegen/pdelab/driver/interpolate.py
+++ b/python/dune/codegen/pdelab/driver/interpolate.py
@@ -1,3 +1,12 @@
+"""Create grid functions
+
+Create all kinds of grid functions. This includes grid functions for the
+boundary condition ('interpolate expression') the exact solution
+('exact_solution') and the grid function that decides where Dirichlet boundary
+conditions are applied ('is_dirichlet')
+"""
+
+
from dune.codegen.generation import (cached,
class_member,
get_counted_variable,
@@ -40,13 +49,15 @@ def interpolate_dirichlet_data(name):
def _grid_function_root_type(identifier):
name_dict = {"exact_solution": "ExactSolution",
- "interpolate_expression": "BoundaryGridFunction"}
+ "interpolate_expression": "BoundaryGridFunction",
+ "is_dirichlet": "BoundaryConditionType"}
return name_dict[identifier]
def _grid_function_root_name(identifier):
name_dict = {"exact_solution": "exactSolution",
- "interpolate_expression": "boundaryGridFunction"}
+ "interpolate_expression": "boundaryGridFunction",
+ "is_dirichlet": "boundaryConditionType"}
return name_dict[identifier]
@@ -56,18 +67,29 @@ def _get_grid_function_method_name(identifier):
return name_dict[identifier]
+#
+# Composite grid function
+#
+
+
@class_member(classtag="driver_block")
-def typedef_composite_grid_function(name, children):
+def typedef_composite_grid_function(name, identifier, children):
templates = ','.join('std::decay_t<decltype(*{})>'.format(c) for c in children)
- return "using {} = Dune::PDELab::CompositeGridFunction<{}>;".format(name, templates)
+ if identifier == "is_dirichlet":
+ composite_type = "Dune::PDELab::CompositeConstraintsParameters"
+ else:
+ composite_type = "Dune::PDELab::CompositeGridFunction"
+ return "using {} = {}<{}>;".format(name, composite_type, templates)
def type_composite_grid_function(identifier, children, root):
if root:
name = _grid_function_root_type(identifier)
+ elif identifier == "is_dirichlet":
+ name = "CompositeConstraintsParameters_{}".format('_'.join(c for c in children))
else:
name = "CompositeGridFunction_{}".format('_'.join(c for c in children))
- typedef_composite_grid_function(name, children)
+ typedef_composite_grid_function(name, identifier, children)
return name
@@ -84,6 +106,11 @@ def define_composite_grid_function(identifier, name, children, root=True):
return "{} = std::make_shared<{}>({});".format(name, composite_gfs_type, ', '.join('*{}'.format(c) for c in children))
+#
+# Function used to generate grid function
+#
+
+
def function_lambda(func):
assert isinstance(func, tuple)
func = func[0]
@@ -100,58 +127,92 @@ def function_lambda(func):
@class_member(classtag="driver_block")
-def typedef_function(name):
- range_type = type_range()
+def typedef_function(name, boolean=False):
leafview_type = type_leafview()
- entity = "typename {}::template Codim<0>::Entity".format(leafview_type)
- coordinate = "typename {}::template Codim<0>::Entity::Geometry::LocalCoordinate".format(leafview_type)
- return "using {} = std::function<{}({}, {})>;".format(name, range_type, entity, coordinate)
+ if boolean:
+ return_type = "bool"
+ entity = "typename {}::Intersection".format(leafview_type)
+ coordinate = "typename {}::Intersection::LocalCoordinate".format(leafview_type)
+ else:
+ return_type = type_range()
+ entity = "typename {}::template Codim<0>::Entity".format(leafview_type)
+ coordinate = "typename {}::template Codim<0>::Entity::Geometry::LocalCoordinate".format(leafview_type)
+ return "using {} = std::function<{}({}, {})>;".format(name, return_type, entity, coordinate)
-def type_function():
- name = "FunctionExpression"
- typedef_function(name)
+def type_function(boolean):
+ if boolean:
+ name = "BoolFunctionExpression"
+ boolean = True
+ else:
+ name = "FunctionExpression"
+ boolean = False
+ typedef_function(name, boolean=boolean)
return name
@class_member(classtag="driver_block")
-def declare_function(name):
- function_type = type_function()
+def declare_function(name, boolean):
+ function_type = type_function(boolean)
return "std::shared_ptr<{}> {};".format(function_type, name)
@preamble(section="vector", kernel="driver_block")
-def define_function(name, func):
- declare_function(name)
- function_type = type_function()
+def define_function(name, func, boolean):
+ declare_function(name, boolean)
+ function_type = type_function(boolean)
fct_lambda = function_lambda(func)
return "{} = std::make_shared<{}> ({});".format(name, function_type, fct_lambda)
@cached
-def name_function(func):
+def name_function(func, boolean):
name = get_counted_variable("functionExpression")
- define_function(name, func)
+ define_function(name, func, boolean)
return name
+#
+# Leaf grid function
+#
+
+
@class_member(classtag="driver_block")
-def typedef_grid_function(name):
+def typedef_grid_function(name, boolean=False):
leafview_type = type_leafview()
range_type = type_range()
- function_type = type_function()
- return "using {} = Dune::PDELab::LocalCallableToGridFunctionAdapter<{}, {}, {}, {}>;".format(name,
- leafview_type,
- range_type,
- 1,
- function_type)
+ function_type = type_function(boolean)
+ if boolean:
+ return "using {} = Dune::PDELab::LocalCallableToBoundaryConditionAdapter<{}>;".format(name,
+ function_type)
+ elif is_stationary:
+ return "using {} = Dune::PDELab::LocalCallableToGridFunctionAdapter<{}, {}, {}, {}>;".format(name,
+ leafview_type,
+ range_type,
+ 1,
+ function_type)
+ else:
+ from dune.codegen.pdelab.driver.gridoperator import type_localoperator
+ lop_type = type_localoperator(get_form_ident())
+ return "using {} = Dune::PDELab::LocalCallableToInstationaryGridFunctionAdapter" \
+ "<{}, {}, {}, {}>;".format(name,
+ leafview_type,
+ range_type,
+ 1,
+ function_type,
+ lop_type)
def type_grid_function(identifier, root):
- name = "GridFunctionLeaf"
+ if identifier == "is_dirichlet":
+ name = "BoolGridFunctionLeaf"
+ boolean = True
+ else:
+ name = "GridFunctionLeaf"
+ boolean = False
if root:
name = _grid_function_root_type(identifier)
- typedef_grid_function(name)
+ typedef_grid_function(name, boolean=boolean)
return name
@@ -165,16 +226,23 @@ def declare_grid_function(identifier, name, root):
def define_grid_function(identifier, name, func, root=True):
declare_grid_function(identifier, name, root)
gv = name_leafview()
- function_name = name_function(func)
+
+ if identifier == "is_dirichlet":
+ boolean = True
+ else:
+ boolean = False
+ function_name = name_function(func, boolean)
+
grid_function_type = type_grid_function(identifier, root)
include_file('dune/pdelab/function/callableadapter.hh', filetag='driver')
- if is_stationary():
+ if identifier == "is_dirichlet":
+ return "{} = std::make_shared<{}>(*{});".format(name, grid_function_type, function_name)
+ elif is_stationary():
return "{} = std::make_shared<{}>({}, *{});".format(name, grid_function_type, gv, function_name)
else:
- # palpo TODO
- assert False
from dune.codegen.pdelab.driver.gridoperator import name_localoperator
lop = name_localoperator(get_form_ident())
+ return "{} = std::make_shared<{}>({}, *{}, *{});".format(name, grid_function_type, gv, function_name, lop)
return "auto {} = Dune::PDELab::makeInstationaryGridFunctionFromCallable({}, {}, {});".format(name,
gv,
lambdaname,
@@ -182,7 +250,6 @@ def define_grid_function(identifier, name, func, root=True):
)
-
def name_grid_function(identifier, element, func, root=True):
assert isinstance(func, tuple)
@@ -207,6 +274,11 @@ def name_grid_function(identifier, element, func, root=True):
return name
+#
+# Name of root of grid function
+#
+
+
def name_grid_function_root(identifier):
element = get_trial_element()
func = preprocess_leaf_data(element, identifier)
@@ -214,6 +286,11 @@ def name_grid_function_root(identifier):
return name
+#
+# Use grid function outside the driver block
+#
+
+
@preamble(section="postprocessing", kernel="main")
def main_typedef_grid_function(name, identifier, treepath):
if len(treepath) == 0:
@@ -259,7 +336,6 @@ def main_define_grid_function(name, identifier, treepath):
return "auto {} = child(*{}, {});".format(name, root_name, ", ".join(i for i in indices))
-
def main_name_grid_function(identifier, treepath):
name = _grid_function_root_name(identifier)
if len(treepath) > 0: