Skip to content
Snippets Groups Projects
Commit 00130d5a authored by René Heß's avatar René Heß
Browse files

[skip ci] Create contstraints parameters through interpolate.py 

Constraints parameter classes are basically just grid functions returning
bools. Going through interpolate.py reduces lots of reduntant code.
parent 5ef2b1de
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
python/dune/codegen/pdelab/driver/constraints.py
+ 3
− 113
View file @ 00130d5a
......@@ -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
#
......
python/dune/codegen/pdelab/driver/gridoperator.py
+ 35
− 0
View file @ 00130d5a
......@@ -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
......
python/dune/codegen/pdelab/driver/instationary.py
+ 49
− 19
View file @ 00130d5a
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():
......
python/dune/codegen/pdelab/driver/interpolate.py
+ 111
− 35
View file @ 00130d5a
"""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:
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Datenschutz und Impressum