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

Hagen Poiseuille example is working fully!

parent c6f171d3
No related branches found
No related tags found
No related merge requests found
......@@ -672,6 +672,8 @@ def define_vector(name):
@preamble
def define_boundary_lambda(boundary, name):
if boundary is None:
return "auto {} = [&](const auto& x){{ return 0.0; }};".format(name)
if boundary.is_global:
return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr)
else:
......@@ -694,28 +696,53 @@ def define_boundary_function(boundary, name):
)
@preamble
def define_compositegfs_parameterfunction(name, *children):
return "Dune::PDELab::CompositeGridFunction<{}> {}({});".format(', '.join('decltype({})'.format(c) for c in children),
name,
', '.join(children))
@symbol
def name_boundary_function(boundary):
define_boundary_function(boundary, "b")
# TODO add namespace data here
return "b"
def name_boundary_function(expr):
from ufl import MixedElement, VectorElement
if isinstance(expr, VectorElement):
element, (_, boundary) = get_constraints_metadata(expr)
from dune.perftool.generation import get_global_context_value
name = get_global_context_value("namedata").get(id(boundary), "zero")
define_boundary_function(boundary, name)
return name
if isinstance(expr, MixedElement):
children = tuple(name_boundary_function(el) for el in expr.sub_elements())
name = '_'.join(children)
define_compositegfs_parameterfunction(name, *children)
return name
# This is a scalar leaf of the ansatz function tree
element, (_, boundary) = get_constraints_metadata(expr)
from dune.perftool.generation import get_global_context_value
name = get_global_context_value("namedata").get(id(boundary), "zero")
define_boundary_function(boundary, name)
return name
@preamble
def interpolate_vector(name):
define_vector(name)
element, (_, boundary) = get_constraints_metadata(_form.coefficients()[0].ufl_element())
b = name_boundary_function(boundary)
element = _form.coefficients()[0].ufl_element()
bf = name_boundary_function(element)
gfs = name_gfs(element)
return "Dune::PDELab::interpolate({}, {}, {});".format(b,
return "Dune::PDELab::interpolate({}, {}, {});".format(bf,
gfs,
name,
)
def maybe_interpolate_vector(name):
element, (_, boundary) = get_constraints_metadata(_form.coefficients()[0].ufl_element())
if boundary:
element = _form.coefficients()[0].ufl_element()
if has_constraints(element):
interpolate_vector(name)
else:
define_vector(name)
......
dune_symlink_to_source_files(FILES hagenpoiseuille_ref.vtu)
add_generated_executable(UFLFILE stokes.ufl
TARGET stokes_numdiff
FORM_COMPILER_ARGS --numerical-jacobian
......@@ -5,4 +7,5 @@ add_generated_executable(UFLFILE stokes.ufl
dune_add_system_test(TARGET stokes_numdiff
INIFILE stokes_numdiff.mini
SCRIPT dune_vtkcompare.py
)
File added
v_bctype = Expression("if (x[0] < 1e-8) return 1; else return 0;", on_intersection=True)
v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
v_dirichlet = Expression("Dune::FieldVector<double, 2> y(0.0); y[0] = 4*x[1]*(1.-x[1]); return y;")
cell = triangle
......
......@@ -2,6 +2,11 @@ __name = stokes_numdiff
lowerleft = 0.0 0.0
upperright = 1.0 1.0
elements = 4 4
elements = 16 16
elementType = simplical
printmatrix = true
printmatrix = false
[wrapper.vtkcompare]
name = {__name}
reference = hagenpoiseuille_ref
extension = vtu
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