From 5b0942d4036a2abae836c1068d23fee60521e9e6 Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Tue, 26 Apr 2016 15:06:23 +0200
Subject: [PATCH] Work poisson_neumann_numdiff example!
---
python/dune/perftool/__init__.py | 4 +-
python/dune/perftool/loopy/transformer.py | 18 +-
python/dune/perftool/pdelab/argument.py | 56 +++++-
python/dune/perftool/pdelab/basis.py | 9 +-
python/dune/perftool/pdelab/driver.py | 10 +-
python/dune/perftool/pdelab/geometry.py | 173 +++++++++++++++++--
python/dune/perftool/pdelab/localoperator.py | 59 +++++--
python/dune/perftool/pdelab/parameter.py | 67 +++++--
python/dune/perftool/pdelab/quadrature.py | 7 +-
python/dune/perftool/ufl/execution.py | 3 +-
test/poisson/CMakeLists.txt | 1 +
test/poisson/poisson_neumann.ufl | 4 +-
12 files changed, 352 insertions(+), 59 deletions(-)
diff --git a/python/dune/perftool/__init__.py b/python/dune/perftool/__init__.py
index 175c0adc..eb279cf1 100644
--- a/python/dune/perftool/__init__.py
+++ b/python/dune/perftool/__init__.py
@@ -1,4 +1,4 @@
class Restriction:
NONE = 0
- POSITIVE = 1
- NEGATIVE = 2
+ INSIDE = 1
+ OUTSIDE = 2
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 6776df96..21c2265d 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -48,8 +48,13 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper):
name = get_global_context_value("namedata")[id(o)]
# Trigger the generation of code for this thing in the parameter class
- from dune.perftool.pdelab.parameter import parameter_function
- parameter_function(name, o)
+ from dune.perftool.pdelab.parameter import (cell_parameter_function,
+ intersection_parameter_function,
+ )
+ if o.on_intersection:
+ intersection_parameter_function(name, o)
+ else:
+ cell_parameter_function(name, o)
# And return a symbol
from pymbolic.primitives import Variable
@@ -151,8 +156,9 @@ def transform_accumulation_term(term, measure, subdomain_id):
# Generate the code for the modified arguments:
for arg in test_ma:
- from dune.perftool.pdelab.argument import name_argumentspace, pymbolic_argument
- accumargs.append(name_argumentspace(arg))
+ from dune.perftool.pdelab.argument import pymbolic_argument
+ from dune.perftool.pdelab.basis import name_lfs
+ accumargs.append(name_lfs(arg.argexpr.element()))
accumargs.append(lfs_iname(arg.argexpr.element(), argcount=arg.argexpr.count()))
# Determine the shape
@@ -215,8 +221,8 @@ def transform_accumulation_term(term, measure, subdomain_id):
name = get_global_context_value("namedata")[id(subdomain_data)]
# Trigger the generation of code for this thing in the parameter class
- from dune.perftool.pdelab.parameter import parameter_function
- parameter_function(name, subdomain_data, t='int')
+ from dune.perftool.pdelab.parameter import intersection_parameter_function
+ intersection_parameter_function(name, subdomain_data, t='int')
code = "if ({} == {})\n {}".format(name, subdomain_id, code)
diff --git a/python/dune/perftool/pdelab/argument.py b/python/dune/perftool/pdelab/argument.py
index 86e1ccba..3ea5c1d8 100644
--- a/python/dune/perftool/pdelab/argument.py
+++ b/python/dune/perftool/pdelab/argument.py
@@ -68,17 +68,35 @@ def pymbolic_trialfunction(ma):
@symbol
-def name_testfunctionspace(*a):
- # TODO
+def name_testfunctionspace():
return "lfsv"
@symbol
-def name_trialfunctionspace(*a):
- # TODO
+def name_trialfunctionspace():
return "lfsu"
+@symbol
+def type_testfunctionspace():
+ return "LFSV"
+
+
+@symbol
+def type_trialfunctionspace():
+ return "LFSU"
+
+
+@symbol
+def name_coefficientcontainer():
+ return "x"
+
+
+@symbol
+def type_coefficientcontainer():
+ return "X"
+
+
def name_argumentspace(ma):
if ma.argexpr.number() == 0:
return name_testfunctionspace(ma)
@@ -110,6 +128,36 @@ def name_jacobian():
return "jac"
+@symbol
+def type_jacobian():
+ return "J"
+
+
@symbol
def name_residual():
return "r"
+
+
+@symbol
+def type_residual():
+ return "R"
+
+
+def name_accumulation_variable():
+ from dune.perftool.generation import get_global_context_value
+ ft = get_global_context_value("form_type")
+ if ft == 'residual':
+ return name_residual()
+ if ft == 'jacobian':
+ return name_jacobian()
+ assert False
+
+
+def type_accumulation_variable():
+ from dune.perftool.generation import get_global_context_value
+ ft = get_global_context_value("form_type")
+ if ft == 'residual':
+ return type_residual()
+ if ft == 'jacobian':
+ return type_jacobian()
+ assert False
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index f8f7484e..2015ac73 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -11,11 +11,12 @@ from dune.perftool.generation import (cached,
symbol,
temporary_variable,
)
-from dune.perftool.pdelab.quadrature import (name_quadrature_position,
+from dune.perftool.pdelab.quadrature import (name_quadrature_position_in_cell,
quadrature_iname,
)
from dune.perftool.pdelab.geometry import (name_dimension,
name_jacobian_inverse_transposed,
+ to_cell_coordinates,
)
from dune.perftool.pdelab.localoperator import (lop_template_ansatz_gfs,
lop_template_test_gfs,
@@ -192,7 +193,7 @@ def evaluate_basis(element, name):
declare_cache_temporary(element, name, which='Function')
cache = name_localbasis_cache(element)
lfs = name_lfs(element)
- qp = name_quadrature_position()
+ qp = name_quadrature_position_in_cell()
instruction(inames=(quadrature_iname(),
),
code='{} = {}.evaluateFunction({}, {}.finiteElement().localBasis());'.format(name,
@@ -201,6 +202,7 @@ def evaluate_basis(element, name):
lfs,
),
assignees=frozenset({name}),
+ read_variables=frozenset({qp}),
)
@@ -216,7 +218,7 @@ def evaluate_reference_gradient(element, name):
declare_cache_temporary(element, name, which='Jacobian')
cache = name_localbasis_cache(element)
lfs = name_lfs(element)
- qp = name_quadrature_position()
+ qp = name_quadrature_position_in_cell()
instruction(inames=(quadrature_iname(),
),
code='{} = {}.evaluateJacobian({}, {}.finiteElement().localBasis());'.format(name,
@@ -225,6 +227,7 @@ def evaluate_reference_gradient(element, name):
lfs,
),
assignees=frozenset({name}),
+ read_variables=frozenset({qp}),
)
diff --git a/python/dune/perftool/pdelab/driver.py b/python/dune/perftool/pdelab/driver.py
index e5c25af6..556fbfb8 100644
--- a/python/dune/perftool/pdelab/driver.py
+++ b/python/dune/perftool/pdelab/driver.py
@@ -339,7 +339,10 @@ def name_constraintscontainer(element, dirichlet):
@preamble
def define_intersection_lambda(expression, name):
- return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr)
+ if expression.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, expression.c_expr)
+ else:
+ return "auto {} = [&](const auto& is, const auto& x){{ {} }};".format(name, expression.c_expr)
@symbol
@@ -601,7 +604,10 @@ def define_vector(name):
@preamble
def define_boundary_lambda(boundary, name):
- return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr)
+ if boundary.is_global:
+ return "auto {} = [&](const auto& x){{ {} }};".format(name, boundary.c_expr)
+ else:
+ return "auto {} = [&](const auto& e, const auto& x){{ {} }};".format(name, boundary.c_expr)
@symbol
diff --git a/python/dune/perftool/pdelab/geometry.py b/python/dune/perftool/pdelab/geometry.py
index bc42d0af..4e0f5ab6 100644
--- a/python/dune/perftool/pdelab/geometry.py
+++ b/python/dune/perftool/pdelab/geometry.py
@@ -1,3 +1,4 @@
+from dune.perftool import Restriction
from dune.perftool.generation import (preamble,
symbol,
temporary_variable,
@@ -8,38 +9,184 @@ from dune.perftool.pdelab.quadrature import (name_quadrature_position,
@symbol
-def name_entitygeometry():
+def name_element_geometry_wrapper():
return 'eg'
@symbol
-def name_entity():
- eg = name_entitygeometry()
- return "{}.entity()".format(eg)
+def type_element_geometry_wrapper():
+ return 'EG'
+
+
+@symbol
+def name_intersection_geometry_wrapper():
+ return 'ig'
+
+
+@symbol
+def type_intersection_geometry_wrapper():
+ return 'IG'
+
+
+def name_geometry_wrapper():
+ """ Selects the 'native' geometry wrapper of the kernel """
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+
+ if it == 'cell':
+ return name_element_geometry_wrapper()
+ if it == 'exterior_facet' or it == 'interior_facet':
+ return name_intersection_geometry_wrapper()
+ assert False
+
+
+def type_geometry_wrapper():
+ """ Selects the 'native' geometry wrapper of the kernel """
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+
+ if it == 'cell':
+ return type_element_geometry_wrapper()
+ if it == 'exterior_facet' or it == 'interior_facet':
+ return type_intersection_geometry_wrapper()
+ assert False
@preamble
-def define_geometry(name):
- eg = name_entitygeometry()
+def define_restricted_cell(name, restriction):
+ ig = name_intersection_geometry_wrapper()
+ which = "inside" if restriction == Restriction.INSIDE else "outside"
+ return "const auto& {} = {}.{}();".format(name,
+ ig,
+ which,
+ )
+
+
+@symbol
+def _name_cell(restriction):
+ if restriction == Restriction.NONE:
+ eg = name_element_geometry_wrapper()
+ return "{}.entity()".format(eg)
+
+ which = "inside" if restriction == Restriction.INSIDE else "outside"
+ name = "cell_{}".format(which)
+ define_restricted_cell(name, restriction)
+ return name
+
+
+def name_cell():
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+
+ if it == 'cell':
+ r = Restriction.NONE
+ if it == 'exterior_facet':
+ r = Restriction.INSIDE
+ if it == 'interior_facet':
+ raise NotImplementedError
+
+ return _name_cell(r)
+
+
+@preamble
+def define_cell_geometry(name):
+ eg = name_element_geometry_wrapper()
return "auto {} = {}.geometry();".format(name,
eg
)
@symbol
+def name_cell_geometry():
+ define_cell_geometry("cell_geo")
+ return "cell_geo"
+
+
+@preamble
+def define_intersection_geometry(name):
+ ig = name_intersection_geometry_wrapper()
+ return "auto {} = {}.geometry();".format(name,
+ ig,
+ )
+
+
+@symbol
+def name_intersection_geometry():
+ define_intersection_geometry("is_geo")
+ return "is_geo"
+
+
+@symbol
+def name_intersection():
+ ig = name_intersection_geometry_wrapper()
+ return "{}.intersection()".format(ig)
+
+
def name_geometry():
- define_geometry("geo")
- return "geo"
+ """ Selects the 'native' geometry of the kernel """
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+
+ if it == 'cell':
+ return name_cell_geometry()
+ if it == 'exterior_facet' or it == 'interior_facet':
+ return name_intersection_geometry()
+ assert False
+
+
+@preamble
+def define_in_cell_geometry(restriction, name):
+ cell = _name_cell(restriction)
+ ig = name_intersection_geometry_wrapper()
+ which = "In" if restriction == Restriction.INSIDE else "Out"
+ return "auto {} = {}.geometryIn{}side();".format(name,
+ ig,
+ which
+ )
@symbol
-def type_geometry():
- return 'EG'
+def name_in_cell_geometry(restriction):
+ assert restriction is not Restriction.NONE
+
+ name = "geo_in_{}side".format("in" if restriction is Restriction.INSIDE else "out")
+ define_in_cell_geometry(restriction, name)
+ return name
+
+
+def apply_in_cell_transformation(name, local, restriction):
+ geo = name_in_cell_geometry(restriction)
+ return quadrature_preamble("{} = {}.global({});".format(name,
+ geo,
+ local,
+ ),
+ assignees=name,
+ )
+
+
+@symbol
+def name_in_cell_coordinates(local, basename, restriction):
+ name = "{}_in_inside".format(basename)
+ temporary_variable(name, shape=(name_dimension(),), shape_impl=("fv",))
+ apply_in_cell_transformation(name, local, restriction=Restriction.INSIDE)
+ return name
+
+
+def to_cell_coordinates(local, basename):
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+
+ if it == 'cell':
+ return local
+ if it == 'exterior_facet':
+ return name_in_cell_coordinates(local, basename, Restriction.INSIDE)
+ if it == 'interior_facet':
+ raise NotImplementedError
@preamble
def define_dimension(name):
- geo = type_geometry()
+ geo = type_geometry_wrapper()
return 'const int {} = {}::Entity::dimension;'.format(name, geo)
@@ -51,7 +198,7 @@ def name_dimension():
@symbol
def type_jacobian_inverse_transposed():
- geo = type_geometry()
+ geo = type_element_geometry_wrapper()
return "typename {}::Geometry::JacobianInverseTransposed".format(geo)
@@ -66,7 +213,7 @@ def define_jacobian_inverse_transposed_temporary(name, shape, shape_impl):
def define_jacobian_inverse_transposed(name):
dim = name_dimension()
temporary_variable(name, decl_method=define_jacobian_inverse_transposed_temporary, shape=(dim, dim))
- geo = name_geometry()
+ geo = name_cell_geometry()
pos = name_quadrature_position()
return quadrature_preamble("{} = {}.jacobianInverseTransposed({});".format(name,
geo,
diff --git a/python/dune/perftool/pdelab/localoperator.py b/python/dune/perftool/pdelab/localoperator.py
index 3301357b..370ef617 100644
--- a/python/dune/perftool/pdelab/localoperator.py
+++ b/python/dune/perftool/pdelab/localoperator.py
@@ -117,20 +117,55 @@ def assembly_routine_signature():
form_type = get_global_context_value("form_type")
aj = "alpha" if form_type == 'residual' else "jacobian"
-
- if integral_type == 'cell':
- return ['template<typename EG, typename LFSU, typename X, typename LFSV, typename R>',
- 'void {}_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const'.format(aj)]
-
- if integral_type == 'exterior_facet':
- return ['template<typename IG, typename LFSV0, typename X, typename LFSV1, typename R>',
- 'void {}_boundary(const IG& ig, const LFSV0& lfsv0, const X& x, const LFSV1& lfsv1, R& r) const'.format(aj)]
+ which = ufl_measure_to_pdelab_measure(integral_type).lower()
if integral_type == 'interior_facet':
- return ['template<typename IG, typename LFSV0_S, typename X, typename LFSV1_S, typename LFSV0_N, typename R, typename LFSV1_N>',
- 'void {}_skeleton(const IG& ig, const LFSV0_S& lfsv0_s, const X& x_s, const LFSV1_S& lfsv1_s, const LFSV0_N& lfsv0_n, const X& x_n, const LFSV1_N& lfsv1_n, R& r_s, R& r_n) const'.format(aj)]
-
- assert False
+ raise NotImplementedError
+
+ from dune.perftool.pdelab.geometry import (name_geometry_wrapper,
+ type_geometry_wrapper,
+ )
+ geot = type_geometry_wrapper()
+ geo = name_geometry_wrapper()
+
+ from dune.perftool.pdelab.argument import (name_accumulation_variable,
+ type_accumulation_variable,
+ name_coefficientcontainer,
+ type_coefficientcontainer,
+ name_testfunctionspace,
+ type_testfunctionspace,
+ name_trialfunctionspace,
+ type_trialfunctionspace,
+ )
+ lfsut = type_trialfunctionspace()
+ lfsu = name_trialfunctionspace()
+ lfsvt = type_testfunctionspace()
+ lfsv = name_testfunctionspace()
+ cct = type_coefficientcontainer()
+ cc = name_coefficientcontainer()
+ avt = type_accumulation_variable()
+ av = name_accumulation_variable()
+
+ return ['template<typename {}, typename {}, typename {}, typename {}, typename {}>'.format(geot,
+ lfsut,
+ cct,
+ lfsvt,
+ avt,
+ ),
+ 'void {}_{}(const {}& {}, const {}& {}, const {}& {}, const {}& {}, {}& {}) const'.format(aj,
+ which,
+ geot,
+ geo,
+ lfsut,
+ lfsu,
+ cct,
+ cc,
+ lfsvt,
+ lfsv,
+ avt,
+ av,
+ )
+ ]
def generate_kernel(integral):
diff --git a/python/dune/perftool/pdelab/parameter.py b/python/dune/perftool/pdelab/parameter.py
index 9b4885d3..dc717063 100644
--- a/python/dune/perftool/pdelab/parameter.py
+++ b/python/dune/perftool/pdelab/parameter.py
@@ -5,11 +5,15 @@ from dune.perftool.generation import (class_basename,
constructor_parameter,
generator_factory,
initializer_list,
+ preamble,
symbol,
temporary_variable
)
-from dune.perftool.pdelab.geometry import name_entity
+from dune.perftool.pdelab.geometry import (name_cell,
+ name_intersection,
+ )
from dune.perftool.pdelab.quadrature import (name_quadrature_position,
+ name_quadrature_position_in_cell,
quadrature_preamble,
)
from dune.perftool.cgen.clazz import AccessModifier
@@ -39,14 +43,16 @@ def name_paramclass():
@class_member("parameterclass", access=AccessModifier.PUBLIC)
-def define_parameter_function_class_member(name, expr, t):
- result = ["template<typename E, typename X>",
- "{} {}(const E& e, const X& local) const".format(t, name),
+def define_parameter_function_class_member(name, expr, t, cell):
+ geot = "E" if cell else "I"
+ geo = geot.lower()
+ result = ["template<typename {}, typename X>".format(geot),
+ "{} {}(const {}& {}, const X& local) const".format(t, name, geot, geo),
"{",
]
if expr.is_global:
- result.append(" auto x = e.geometry().global(local);")
+ result.append(" auto x = {}.geometry().global(local);".format(geo))
else:
result.append(" auto x = local;")
@@ -56,21 +62,56 @@ def define_parameter_function_class_member(name, expr, t):
return result
-def evaluate_parameter_function(name):
+def evaluate_cell_parameter_function(name):
param = name_paramclass()
- entity = name_entity()
- pos = name_quadrature_position()
+ entity = name_cell()
+ pos = name_quadrature_position_in_cell()
return quadrature_preamble('{} = {}.{}({}, {});'.format(name,
name_paramclass(),
name,
entity,
pos,
),
- assignees=frozenset({name})
+ assignees=frozenset({name}),
+ read_variables=frozenset({pos}),
+ )
+
+
+def evaluate_intersection_parameter_function(name):
+ # Check that this is not a volume term, as that would not be well-defined
+ from dune.perftool.generation import get_global_context_value
+ it = get_global_context_value("integral_type")
+ assert it is not 'cell'
+
+ param = name_paramclass()
+ intersection = name_intersection()
+ pos = name_quadrature_position()
+ return quadrature_preamble('{} = {}.{}({}, {});'.format(name,
+ name_paramclass(),
+ name,
+ intersection,
+ pos,
+ ),
+ assignees=frozenset({name}),
+ read_variables=frozenset({pos}),
)
-def parameter_function(name, expr, t='double'):
- temporary_variable(name, shape=())
- define_parameter_function_class_member(name, expr, t)
- evaluate_parameter_function(name)
+def define_parameter_temporary(t):
+ @preamble
+ def _define(name, shape, shape_impl):
+ assert len(shape) == 0
+ return "{} {};".format(t, name)
+ return _define
+
+
+def cell_parameter_function(name, expr, t='double'):
+ temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
+ define_parameter_function_class_member(name, expr, t, True)
+ evaluate_cell_parameter_function(name)
+
+
+def intersection_parameter_function(name, expr, t='double'):
+ temporary_variable(name, shape=(), decl_method=define_parameter_temporary(t))
+ define_parameter_function_class_member(name, expr, t, False)
+ evaluate_intersection_parameter_function(name)
diff --git a/python/dune/perftool/pdelab/quadrature.py b/python/dune/perftool/pdelab/quadrature.py
index b4555eef..2e717f76 100644
--- a/python/dune/perftool/pdelab/quadrature.py
+++ b/python/dune/perftool/pdelab/quadrature.py
@@ -29,6 +29,12 @@ def name_quadrature_position():
return "{}.position()".format(qp)
+@symbol
+def name_quadrature_position_in_cell():
+ from dune.perftool.pdelab.geometry import to_cell_coordinates
+ return to_cell_coordinates(name_quadrature_position(), name_quadrature_point())
+
+
@symbol
def name_quadrature_weight():
qp = name_quadrature_point()
@@ -63,7 +69,6 @@ def name_order():
return "2"
-@cached
def quadrature_loop_statement():
include_file('dune/pdelab/common/quadraturerules.hh', filetag='operatorfile')
qp = name_quadrature_point()
diff --git a/python/dune/perftool/ufl/execution.py b/python/dune/perftool/ufl/execution.py
index a3397c0a..28cacedc 100644
--- a/python/dune/perftool/ufl/execution.py
+++ b/python/dune/perftool/ufl/execution.py
@@ -48,10 +48,11 @@ _dg0 = FiniteElement("DG", "triangle", 0)
class Expression(Coefficient):
- def __init__(self, expr, is_global=True):
+ def __init__(self, expr, is_global=True, on_intersection=False):
assert isinstance(expr, str)
self.c_expr = expr
self.is_global = is_global
+ self.on_intersection = on_intersection
# Initialize a coefficient with a dummy finite element map.
Coefficient.__init__(self, _dg0)
diff --git a/test/poisson/CMakeLists.txt b/test/poisson/CMakeLists.txt
index 2ccd6cce..c4e3a69d 100644
--- a/test/poisson/CMakeLists.txt
+++ b/test/poisson/CMakeLists.txt
@@ -31,4 +31,5 @@ add_generated_executable(UFLFILE poisson_neumann.ufl
dune_add_system_test(TARGET poisson_neumann_numdiff
INIFILE poisson.mini
+ SCRIPT dune_vtkcompare.py
)
diff --git a/test/poisson/poisson_neumann.ufl b/test/poisson/poisson_neumann.ufl
index 385a1bb8..e86bb53b 100644
--- a/test/poisson/poisson_neumann.ufl
+++ b/test/poisson/poisson_neumann.ufl
@@ -1,7 +1,7 @@
f = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return 4.*(1.-c.two_norm2())*std::exp(-1.*c.two_norm2());")
g = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; return std::exp(-1.*c.two_norm2());")
-j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return s*4*x[1]*c.two_norm2()*std::exp(-1.*c.two_norm2());")
-bctype = Expression("if ((x[1]<1e-8) || (x[1]>1-1e-8)) return 0; else return 1;")
+j = Expression("Dune::FieldVector<double,2> c(0.5); c-= x; double s; if (x[1]>0.5) s=1.; else s=-1.; return -2.*s*(x[1]-0.5)*std::exp(-1.*c.two_norm2());", on_intersection=True)
+bctype = Expression("if ((x[1]<1e-8) || (x[1]>1.-1e-8)) return 0; else return 1;", on_intersection=True)
V = FiniteElement("CG", "triangle", 1, dirichlet_expression=g, dirichlet_constraints=bctype)
u = TrialFunction(V)
--
GitLab