From 1313757e16be734fd388212bc00035bd344f13f6 Mon Sep 17 00:00:00 2001
From: Marcel Koch <marcel.koch@uni-muenster.de>
Date: Mon, 31 Jul 2017 14:53:04 +0200
Subject: [PATCH] Adds some comments
---
dune/perftool/localbasiscache.hh | 2 +-
.../dune/perftool/blockstructured/__init__.py | 2 +-
.../dune/perftool/blockstructured/argument.py | 1 +
python/dune/perftool/blockstructured/basis.py | 1 +
.../dune/perftool/blockstructured/geometry.py | 6 ++++++
python/dune/perftool/blockstructured/tools.py | 17 +++++++++++++++++
6 files changed, 27 insertions(+), 2 deletions(-)
diff --git a/dune/perftool/localbasiscache.hh b/dune/perftool/localbasiscache.hh
index f9370436..6a9312e7 100644
--- a/dune/perftool/localbasiscache.hh
+++ b/dune/perftool/localbasiscache.hh
@@ -31,7 +31,7 @@ class CacheReturnProxy
const std::vector<T>* v;
};
-//TODO FS hier entfernen und nur LocalBasisType benutzen
+
template<typename LocalBasisType >
class LocalBasisCacheWithoutReferences
{
diff --git a/python/dune/perftool/blockstructured/__init__.py b/python/dune/perftool/blockstructured/__init__.py
index a3ed2718..17d70256 100644
--- a/python/dune/perftool/blockstructured/__init__.py
+++ b/python/dune/perftool/blockstructured/__init__.py
@@ -24,7 +24,7 @@ class BlockStructuredInterface(PDELabInterface):
# Local function space related generator functions
#
- # TODO current way to squeeze subelem iname in, not really understandable
+ # TODO current way to squeeze subelem iname in, not really ideal
def lfs_inames(self, element, restriction, number=None, context=''):
return lfs_inames(element, restriction, number, context) + sub_element_inames()
diff --git a/python/dune/perftool/blockstructured/argument.py b/python/dune/perftool/blockstructured/argument.py
index 0f1ea01f..7ca3a94c 100644
--- a/python/dune/perftool/blockstructured/argument.py
+++ b/python/dune/perftool/blockstructured/argument.py
@@ -19,4 +19,5 @@ def pymbolic_coefficient(container, lfs, element, index):
if not isinstance(lfs, prim.Expression):
lfs = prim.Variable(lfs)
+ # use higher order FEM index instead of Q1 index
return prim.Call(CoefficientAccess(container), (lfs, micro_index_to_macro_index(element, index),))
diff --git a/python/dune/perftool/blockstructured/basis.py b/python/dune/perftool/blockstructured/basis.py
index 4533179e..4751b89c 100644
--- a/python/dune/perftool/blockstructured/basis.py
+++ b/python/dune/perftool/blockstructured/basis.py
@@ -20,6 +20,7 @@ from dune.perftool.pdelab.quadrature import pymbolic_quadrature_position_in_cell
from dune.perftool.pdelab.spaces import type_gfs
+# define FE basis explicitly in localoperator
@backend(interface="typedef_localbasis", name="blockstructured")
@class_member(classtag="operator")
def typedef_localbasis(element, name):
diff --git a/python/dune/perftool/blockstructured/geometry.py b/python/dune/perftool/blockstructured/geometry.py
index 8ebf8d2f..58512145 100644
--- a/python/dune/perftool/blockstructured/geometry.py
+++ b/python/dune/perftool/blockstructured/geometry.py
@@ -12,25 +12,31 @@ from dune.perftool.blockstructured.tools import sub_element_inames
import pymbolic.primitives as prim
+# scale determinant according to the order of the blockstructure
def pymbolic_jacobian_determinant():
return prim.Quotient(prim.Variable(name_jacobian_determinant()),
prim.Power(get_option("number_of_blocks"), local_dimension()))
+# scale Jacobian according to the order of the blockstructure
def pymbolic_jacobian_inverse_transposed(i, j, restriction):
return prim.Product((get_option("number_of_blocks"),
prim.Subscript(prim.Variable(name_jacobian_inverse_transposed(restriction)), (j, i))))
+# scale determinant according to the order of the blockstructure
def pymbolic_facet_jacobian_determinant():
return prim.Quotient(prim.Variable(name_facet_jacobian_determinant()),
prim.Power(get_option("number_of_blocks"), local_dimension()))
+# translate a point in the micro element into macro coordinates
def define_point_in_macro(name, point_in_micro):
dim = local_dimension()
temporary_variable(name, shape_impl=('fv',), shape=(dim,))
+ # point_macro = (point_micro + index_micro) / number_of_blocks
+ # where index_micro = tensor index of the micro element
subelem_inames = sub_element_inames()
for i in range(dim):
if isinstance(point_in_micro, prim.Subscript):
diff --git a/python/dune/perftool/blockstructured/tools.py b/python/dune/perftool/blockstructured/tools.py
index cf247599..b7d83f6e 100644
--- a/python/dune/perftool/blockstructured/tools.py
+++ b/python/dune/perftool/blockstructured/tools.py
@@ -15,6 +15,8 @@ from dune.perftool.options import get_option
import pymbolic.primitives as prim
+# add inames over the micro elements in tensor representation,
+# i.e. each element has (i_1,i_2,...,i_d) indices
@iname
def sub_element_inames():
name = "subel"
@@ -27,20 +29,30 @@ def sub_element_inames():
return inames
+# define inames for boundary integration
+# In the case of boundary integrationsince we have only d-1 loops,
+# but we need always d inames to compute the macro index.
+# Therefore we must find out which iname is constant.
+# Since loo.py cannot handle if-else blocks very well, this whole
+# computation seems a bit cumbersome.
def sub_facet_inames():
subelem_inames = sub_element_inames()
center = pymbolic_in_cell_coordinates(prim.Variable(name_localcenter()), Restriction.NEGATIVE)
+ # check if iname[index] must be constant or not
def predicate(index):
return prim.Comparison(prim.Subscript(center, index), "==", 0.5)
def conditional_instruction(index):
+ # instruction for not constant iname
+ # special case for the third iname
instruction(assignee=prim.Variable(inames[index]),
expression=prim.Variable(subelem_inames[1 if index == 2 else 0]),
within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
predicates=frozenset([predicate(index)])
)
+ # instruction for constant iname
instruction(assignee=prim.Variable(inames[index]),
expression=prim.Product(((k - 1), prim.Subscript(center, (index,)))),
within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
@@ -59,6 +71,7 @@ def sub_facet_inames():
else:
inames = inames + ("y",)
temporary_variable(inames[1])
+ # one additional condition is needed in 3d for the second iname
instruction(assignee=prim.Variable(inames[1]),
expression=prim.Variable(subelem_inames[0]),
within_inames=frozenset(subelem_inames).union(frozenset(quadrature_inames())),
@@ -82,18 +95,22 @@ def sub_facet_inames():
return inames
+# compute sequential index for given tensor index, the same as index in base-k to base-10
def tensor_index_to_sequential_index(indices, k):
return prim.Sum(tuple(index * k ** i for i, index in enumerate(indices)))
+# compute tensor index for given sequential index, only for Q1
def to_tensor_index(iname, k):
return tuple(prim.RightShift(prim.BitwiseAnd((prim.Variable(iname), 2**n)), n) for n in range(world_dimension()))
+# compute tensor index for given sequential index, the same as index in base-10 to base-k
def to_tensor_index_v2(iname, k):
return tuple(prim.Remainder(prim.Variable(iname) / (k**i), k) for i in range(world_dimension()))
+# compute index for higher order FEM for a given Q1 index
def micro_index_to_macro_index(element, iname):
it = get_global_context_value("integral_type")
if it == "cell":
--
GitLab