From 589f840c8883833f8ed5e663872fae27846e090b Mon Sep 17 00:00:00 2001
From: Dominic Kempf <dominic.kempf@iwr.uni-heidelberg.de>
Date: Fri, 9 Sep 2016 16:59:55 +0200
Subject: [PATCH] Implement symmetry for TensorElements
and an example with stokes stress
---
python/dune/perftool/loopy/transformer.py | 2 +-
python/dune/perftool/pdelab/basis.py | 21 ++++++++++++++-----
test/stokes/CMakeLists.txt | 5 +++++
test/stokes/stokes_stress_sym.mini | 18 ++++++++++++++++
test/stokes/stokes_stress_sym.ufl | 25 +++++++++++++++++++++++
5 files changed, 65 insertions(+), 6 deletions(-)
create mode 100644 test/stokes/stokes_stress_sym.mini
create mode 100644 test/stokes/stokes_stress_sym.ufl
diff --git a/python/dune/perftool/loopy/transformer.py b/python/dune/perftool/loopy/transformer.py
index 0129f4d7..2e0f3e19 100644
--- a/python/dune/perftool/loopy/transformer.py
+++ b/python/dune/perftool/loopy/transformer.py
@@ -98,7 +98,7 @@ class UFL2LoopyVisitor(ModifiedTerminalTracker, UFL2PymbolicMapper, GeometryMapp
from dune.perftool.pdelab.geometry import dimension_iname
from dune.perftool.pdelab.basis import lfs_child
idims = tuple(dimension_iname(context='arg', count=i) for i in range(len(subel.value_shape())))
- lfs = lfs_child(lfs, idims, shape=subel.value_shape())
+ lfs = lfs_child(lfs, idims, shape=subel.value_shape(), symmetry=subel.symmetry)
residual_shape[icount] = name_lfs_bound(name_leaf_lfs(subel.sub_elements()[0], ma.restriction))
subel = subel.sub_elements()[0]
else:
diff --git a/python/dune/perftool/pdelab/basis.py b/python/dune/perftool/pdelab/basis.py
index edbf22a2..4c6c4011 100644
--- a/python/dune/perftool/pdelab/basis.py
+++ b/python/dune/perftool/pdelab/basis.py
@@ -75,14 +75,25 @@ def define_lfs(name, father, child):
return "auto {} = child({}, _{});".format(name, father, child)
-def lfs_child(lfs, children, shape=None):
+def lfs_child(lfs, children, shape=None, symmetry=False):
from pymbolic.primitives import Call, Product, Sum
# Old pre-TensorElement implementation kept for comaptibility
if shape is None:
indices = (Variable(children[0]),)
else:
- children = tuple(Variable(c) for c in children)
- indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
+ if symmetry and len(children) == 2:
+ # I do not want to think about tensors of rank > 2 right now
+ i, j = children
+ if i>j:
+ j, i = i, j
+ i = Variable(i)
+ j = Variable(j)
+ n = len(children)
+
+ indices = (Sum((Product((n-1, i)), Product((.5, i, 1-i)), j)),)
+ else:
+ children = tuple(Variable(c) for c in children)
+ indices = (Sum(tuple(Product((Product(tuple(shape[j] for j in range(i))), child)) for i, child in enumerate(children))),)
from dune.perftool.loopy.functions import LFSChild
return Call(LFSChild(lfs), indices)
@@ -367,7 +378,7 @@ def evaluate_coefficient(element, name, container, restriction, component):
if isinstance(sub_element, (VectorElement, TensorElement)):
# TODO we need a concept of symmetry here
- lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape))
+ lfs = lfs_child(lfs, idims, shape=shape_as_pymbolic(shape), symmetry=element.symmetry)
from dune.perftool.pdelab.argument import pymbolic_coefficient
coeff = pymbolic_coefficient(container, lfs, index)
@@ -406,7 +417,7 @@ def evaluate_coefficient_gradient(element, name, container, restriction, compone
if isinstance(sub_element, (VectorElement, TensorElement)):
# TODO we need a concept of symmetry here
- lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]))
+ lfs = lfs_child(lfs, idims[:-1], shape=shape_as_pymbolic(shape[:-1]), symmetry=element.symmetry)
from dune.perftool.pdelab.argument import pymbolic_coefficient
coeff = pymbolic_coefficient(container, lfs, index)
diff --git a/test/stokes/CMakeLists.txt b/test/stokes/CMakeLists.txt
index b8bc4b0a..0553ee4c 100644
--- a/test/stokes/CMakeLists.txt
+++ b/test/stokes/CMakeLists.txt
@@ -21,3 +21,8 @@ dune_add_formcompiler_system_test(UFLFILE stokes_stress.ufl
BASENAME stokes_stress
INIFILE stokes_stress.mini
)
+
+dune_add_formcompiler_system_test(UFLFILE stokes_stress_sym.ufl
+ BASENAME stokes_stress_sym
+ INIFILE stokes_stress_sym.mini
+ )
diff --git a/test/stokes/stokes_stress_sym.mini b/test/stokes/stokes_stress_sym.mini
new file mode 100644
index 00000000..9038ebfa
--- /dev/null
+++ b/test/stokes/stokes_stress_sym.mini
@@ -0,0 +1,18 @@
+__name = stokes_stress_sym_{__exec_suffix}
+__exec_suffix = symdiff, numdiff | expand num
+
+lowerleft = 0.0 0.0
+upperright = 1.0 1.0
+elements = 16 16
+elementType = simplical
+printmatrix = false
+
+[wrapper.vtkcompare]
+name = {__name}
+reference = hagenpoiseuille_ref
+extension = vtu
+
+[formcompiler]
+numerical_jacobian = 0, 1 | expand num
+exact_solution_expression = g
+compare_l2errorsquared = 1e-6
diff --git a/test/stokes/stokes_stress_sym.ufl b/test/stokes/stokes_stress_sym.ufl
new file mode 100644
index 00000000..47e41380
--- /dev/null
+++ b/test/stokes/stokes_stress_sym.ufl
@@ -0,0 +1,25 @@
+v_bctype = Expression("if (x[0] < 1. - 1e-8) return 1; else return 0;", on_intersection=True)
+g_v = Expression(("4*x[1]*(1.-x[1])", "0.0"))
+g_p = Expression("8*(1.-x[0])")
+g_S = Expression(("0.0", "-8*x[1] + 4.", "0.0"))
+g = g_v * g_p * g_S
+
+cell = triangle
+P2 = VectorElement("Lagrange", cell, 2, dirichlet_constraints=v_bctype, dirichlet_expression=g_v)
+P1 = FiniteElement("Lagrange", cell, 1)
+P2_stress = TensorElement("DG", cell, 1, symmetry={(0, 1): (1, 0)})
+
+TH = MixedElement(P2, P1, P2_stress)
+
+v, q, T = TestFunctions(TH)
+u, p, S = TrialFunctions(TH)
+
+n = FacetNormal(triangle)('+')
+ds = ds(subdomain_data=v_bctype, subdomain_id=0)
+
+r = (inner(grad(v), S) + inner(2*sym(grad(u)) - S, T) - div(v)*p - q*div(u))*dx \
+ - inner(grad(u).T*n, v)*ds
+# \
+# + inner(S.T*n, v)*ds
+
+forms = [r]
--
GitLab