Commit d4805639 authored by Dominic Kempf's avatar Dominic Kempf
Browse files

Merge branch 'bugfix/flipping' into 'master'

Fix Symmetry of Fiber Operator

See merge request !2
parents e0c4c4e2 7e30239f
......@@ -297,8 +297,8 @@ class EulerBernoulli2DLocalOperator
using namespace Dune::Indices;
auto child_0 = child(lfsu_s, _0);
auto child_1 = child(lfsu_s, _1);
auto cellgeo_s = (flipped ? ig.outside() : ig.inside()).geometry();
auto cellgeo_n = (flipped ? ig.inside() : ig.outside()).geometry();
auto cellgeo_s = ig.inside().geometry();
auto cellgeo_n = ig.outside().geometry();
using Evaluator = BasisEvaluator<typename LFSU::template Child<0>::Type::Traits::FiniteElementType::Traits::LocalBasisType>;
Evaluator basis_s(child_0.finiteElement().localBasis());
......@@ -320,8 +320,8 @@ class EulerBernoulli2DLocalOperator
for (std::size_t c=0; c<2; ++c)
for (std::size_t i=0; i<child_0.size(); i++)
{
u_s[c] += (flipped ? x_n : x_s)(child(lfsu_s, c), i) * basis_s.function(i);
u_n[c] += (flipped ? x_s : x_n)(child(lfsu_n, c), i) * basis_n.function(i);
u_s[c] += x_s(child(lfsu_s, c), i) * basis_s.function(i);
u_n[c] += x_n(child(lfsu_n, c), i) * basis_n.function(i);
}
// And the jacobian of displacement
......@@ -330,8 +330,8 @@ class EulerBernoulli2DLocalOperator
for (std::size_t d=0; d<2; ++d)
for (std::size_t i=0; i<child_0.size(); ++i)
{
d1u_s[c][d] += (flipped ? x_n : x_s)(child(lfsu_s, c), i) * basis_s.jacobian(i, d);
d1u_n[c][d] += (flipped ? x_s : x_n)(child(lfsu_n, c), i) * basis_n.jacobian(i, d);
d1u_s[c][d] += x_s(child(lfsu_s, c), i) * basis_s.jacobian(i, d);
d1u_n[c][d] += x_n(child(lfsu_n, c), i) * basis_n.jacobian(i, d);
}
// And finally the hessian of the displacement
......@@ -341,8 +341,8 @@ class EulerBernoulli2DLocalOperator
for (std::size_t d1=0; d1<2; ++d1)
for (std::size_t i=0; i<child_0.size(); ++i)
{
d2u_s[c][d0][d1] += (flipped ? x_n : x_s)(child(lfsu_s, c), i) * basis_s.hessian(i, d0, d1);
d2u_n[c][d0][d1] += (flipped ? x_s : x_n)(child(lfsu_n, c), i) * basis_n.hessian(i, d0, d1);
d2u_s[c][d0][d1] += x_s(child(lfsu_s, c), i) * basis_s.hessian(i, d0, d1);
d2u_n[c][d0][d1] += x_n(child(lfsu_n, c), i) * basis_n.hessian(i, d0, d1);
}
// Evaluate the jacobian inverse transposed in inside/outside cell
......@@ -359,7 +359,11 @@ class EulerBernoulli2DLocalOperator
auto I = (d*d*d) / 12.0;
// Compute the penalty factor
auto penalty = beta / ig.geometry().volume();
auto fiber_segment_s = fibintersection.element_fibre_intersections.find(is.index(ig.inside()));
auto fiber_segment_n = fibintersection.element_fibre_intersections.find(is.index(ig.outside()));
auto len_s = (fibre->eval(fiber_segment_s->second.second) - fibre->eval(fiber_segment_s->second.first)).two_norm();
auto len_n = (fibre->eval(fiber_segment_n->second.second) - fibre->eval(fiber_segment_n->second.first)).two_norm();
auto penalty = beta / std::min(len_s, len_n);
// The needed tangential derivative quantities. These expressions are generated
// using the generate_tangential_derivatives Python script.
......@@ -380,12 +384,12 @@ class EulerBernoulli2DLocalOperator
auto dt2vn_s_0 = (t[1] * (jit_s[1][1] * (jit_s[1][1] * basis_s.hessian(i, 1, 1) + jit_s[1][0] * basis_s.hessian(i, 1, 0)) + jit_s[1][0] * (jit_s[1][1] * basis_s.hessian(i, 0, 1) + jit_s[1][0] * basis_s.hessian(i, 0, 0))) * (-1) * t[1] + t[0] * (jit_s[0][1] * (jit_s[1][1] * basis_s.hessian(i, 1, 1) + jit_s[1][0] * basis_s.hessian(i, 1, 0)) + jit_s[0][0] * (jit_s[1][1] * basis_s.hessian(i, 0, 1) + jit_s[1][0] * basis_s.hessian(i, 0, 0))) * (-1) * t[1]) * t[1] + (t[1] * (jit_s[1][1] * (jit_s[0][1] * basis_s.hessian(i, 1, 1) + jit_s[0][0] * basis_s.hessian(i, 1, 0)) + jit_s[1][0] * (jit_s[0][1] * basis_s.hessian(i, 0, 1) + jit_s[0][0] * basis_s.hessian(i, 0, 0))) * (-1) * t[1] + t[0] * (jit_s[0][1] * (jit_s[0][1] * basis_s.hessian(i, 1, 1) + jit_s[0][0] * basis_s.hessian(i, 1, 0)) + jit_s[0][0] * (jit_s[0][1] * basis_s.hessian(i, 0, 1) + jit_s[0][0] * basis_s.hessian(i, 0, 0))) * (-1) * t[1]) * t[0];
auto dt2vn_s_1 = (t[1] * t[0] * (jit_s[1][1] * (jit_s[1][1] * basis_s.hessian(i, 1, 1) + jit_s[1][0] * basis_s.hessian(i, 1, 0)) + jit_s[1][0] * (jit_s[1][1] * basis_s.hessian(i, 0, 1) + jit_s[1][0] * basis_s.hessian(i, 0, 0))) + t[0] * t[0] * (jit_s[0][1] * (jit_s[1][1] * basis_s.hessian(i, 1, 1) + jit_s[1][0] * basis_s.hessian(i, 1, 0)) + jit_s[0][0] * (jit_s[1][1] * basis_s.hessian(i, 0, 1) + jit_s[1][0] * basis_s.hessian(i, 0, 0)))) * t[1] + (t[1] * t[0] * (jit_s[1][1] * (jit_s[0][1] * basis_s.hessian(i, 1, 1) + jit_s[0][0] * basis_s.hessian(i, 1, 0)) + jit_s[1][0] * (jit_s[0][1] * basis_s.hessian(i, 0, 1) + jit_s[0][0] * basis_s.hessian(i, 0, 0))) + t[0] * t[0] * (jit_s[0][1] * (jit_s[0][1] * basis_s.hessian(i, 1, 1) + jit_s[0][0] * basis_s.hessian(i, 1, 0)) + jit_s[0][0] * (jit_s[0][1] * basis_s.hessian(i, 0, 1) + jit_s[0][0] * basis_s.hessian(i, 0, 0)))) * t[0];
double flip_factor = flipped ? 1.0 : -1.0;
(flipped ? r_s : r_n).accumulate(lfsu_n.child(0), i, E*I*(-sk_dt2un*dtvn_n_0 - flip_factor * dt2vn_n_0 * sk_dtun + flip_factor * penalty * sk_dtun * dtvn_n_0));
(flipped ? r_s : r_n).accumulate(lfsu_n.child(1), i, E*I*(-sk_dt2un*dtvn_n_1 - flip_factor * dt2vn_n_1 * sk_dtun + flip_factor * penalty * sk_dtun * dtvn_n_1));
double flip_factor = flipped ? -1.0 : 1.0;
r_n.accumulate(lfsu_n.child(0), i, -E*I*(- flip_factor * sk_dt2un*dtvn_n_0 + flip_factor * 0.5 * dt2vn_n_0 * sk_dtun + penalty * sk_dtun * dtvn_n_0));
r_n.accumulate(lfsu_n.child(1), i, -E*I*(- flip_factor * sk_dt2un*dtvn_n_1 + flip_factor * 0.5 * dt2vn_n_1 * sk_dtun + penalty * sk_dtun * dtvn_n_1));
(flipped ? r_n : r_s).accumulate(lfsu_s.child(0), i, -E*I*(-sk_dt2un*dtvn_s_0 - flip_factor * dt2vn_s_0 * sk_dtun + flip_factor * penalty * sk_dtun * dtvn_s_0));
(flipped ? r_n : r_s).accumulate(lfsu_s.child(1), i, -E*I*(-sk_dt2un*dtvn_s_1 - flip_factor * dt2vn_s_1 * sk_dtun + flip_factor * penalty * sk_dtun * dtvn_s_1));
r_s.accumulate(lfsu_s.child(0), i, -E*I*(flip_factor * sk_dt2un*dtvn_s_0 + flip_factor * 0.5 * dt2vn_s_0 * sk_dtun - penalty * sk_dtun * dtvn_s_0));
r_s.accumulate(lfsu_s.child(1), i, -E*I*(flip_factor * sk_dt2un*dtvn_s_1 + flip_factor * 0.5 * dt2vn_s_1 * sk_dtun - penalty * sk_dtun * dtvn_s_1));
}
}
}
......
......@@ -69,8 +69,8 @@ class FibreReinforcedBulkOperator
virtual void alpha_skeleton(const IG& ig, const LFSU& lfsu_s, const X& x_s, const LFSV& lfsv_s, const LFSU& lfsu_n, const X& x_n, const LFSV& lfsv_n, R& r_s, R& r_n) const
{
bulkoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
//fibreoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
//bulkoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
fibreoperator.alpha_skeleton(ig, lfsu_s, x_s, lfsv_s, lfsu_n, x_n, lfsv_n, r_s, r_n);
}
private:
......
""" Code generator for some differential geometry quantities.
This is not embedded into the dune-codegen code generation pipeline.
Instead, it is meant to be run by hand and the resulting expressions
are copied and pasted into the 1D fibre operator.
Currently, this script needs two adjustments in codegen to be run:
* In python/ufl/transformations/__init__.py the transformation printing needs to be disabled
* In python/ufl/modified_terminals.py the assertion in the grad handler needs to be disabled
I am too lazy right now to find proper upstream fixes for these minor issues.
"""
from dune.codegen.ufl.execution import *
from dune.codegen.pdelab.restriction import restricted_name
from dune.codegen.pdelab.localoperator import GenericAccumulationMixin
from ufl.algorithms import expand_derivatives, extract_arguments
from ufl.algorithms.apply_restrictions import apply_restrictions
from ufl.algorithms.apply_algebra_lowering import apply_algebra_lowering
from ufl.algorithms.remove_complex_nodes import remove_complex_nodes
from ufl.algorithms.apply_function_pullbacks import apply_function_pullbacks
from dune.codegen.ufl.transformations.indexpushdown import pushdown_indexed
from dune.codegen.ufl.visitor import UFL2LoopyVisitor
from pymbolic.mapper import IdentityMapper
import pymbolic.primitives as prim
def generate_tangential_derivatives():
# The 2 here is the polynomial order of the finite element
FE = VectorElement("CG", triangle, 2)
u = Coefficient(FE, cargo={"name": "u"})
v = TestFunction(FE)
# Tangential and normal vector
t = Coefficient(FE, cargo={"name": "t", "diff": 0, "restriction": 0})
n = perp(t)
cell_quantities = {
"dtut" : inner(t, grad(inner(u, t))),
"dtvt" : inner(t, grad(inner(v, t))),
"dt2un" : inner(t, grad(inner(t, grad(inner(u, n))))),
"dt2vn" : inner(t, grad(inner(t, grad(inner(v, n))))),
}
facet_quantities = {
"sk_dt2un" : avg(inner(t, grad(inner(t, grad(inner(u, n)))))),
"dt2vn_n" : inner(t, grad(inner(t, grad(inner(v('-'), n))))),
"dt2vn_s" : inner(t, grad(inner(t, grad(inner(v('+'), n))))),
"sk_dtun" : jump(inner(t, grad(inner(u, n)))),
"dtvn_n" : inner(t, grad(inner(v('-'), n))),
"dtvn_s" : inner(t, grad(inner(v('+'), n))),
}
def print_code(name, expr, measure):
# Do the preprocessing
expr = apply_function_pullbacks(expr)
expr = expand_derivatives(expr)
expr = apply_algebra_lowering(expr)
expr = remove_complex_nodes(expr)
# This weird check would go away if we reproduced more of the UFL preprocessing here...
if name.startswith('sk'):
expr = apply_restrictions(expr)
expr = pushdown_indexed(expr)
print("\n\n\nQuantity: {}".format(name))
for snippet in ufl_to_code(expr, measure):
print(snippet + '\n\n')
for name, expr in cell_quantities.items():
print_code(name, expr, "cell")
for name, expr in facet_quantities.items():
print_code(name, expr, "interior_facet")
class AdHocVisitor(GenericAccumulationMixin, UFL2LoopyVisitor):
def __init__(self, measure):
UFL2LoopyVisitor.__init__(self, "cell", {})
self.grad_count = 0
self.measure = measure
def __call__(self, o, do_split=False):
if len(extract_arguments(o)) == 0:
return [self._call(o, False)]
else:
rets = []
for info in self.list_accumulation_infos(o):
self.current_info = info
expr = self._call(o, False)
if expr != 0:
rets.append(expr)
return rets
def lfs_inames(self, *args):
return []
def argument(self, o):
info = self.get_accumulation_info(o)
if info != self.current_info[0]:
self.indices = None
return 0
name = restricted_name("basis", self.restriction)
if self.grad_count == 0:
name = "{}.function".format(name)
elif self.grad_count == 1:
name = "{}.jacobian".format(name)
elif self.grad_count == 2:
name = "{}.hessian".format(name)
else:
raise NotImplementedError
indices = self.indices[1:]
self.indices = None
return prim.Call(prim.Variable(name), (prim.Variable("i"),) + indices)
def coefficient(self, o):
name = restricted_name(o.cargo["name"], o.cargo.get("restriction", self.restriction))
if self.grad_count > 0:
if "diff" in o.cargo:
self.indices = None
return o.cargo["diff"]
name = "d{}{}".format(self.grad_count, name)
return prim.Variable(name)
def reference_grad(self, o):
self.grad_count = self.grad_count + 1
ret = self.call(o.ufl_operands[0])
self.grad_count = self.grad_count - 1
return ret
def jacobian_inverse(self, o):
i, j = self.indices
self.indices = None
return prim.Subscript(prim.Variable(restricted_name("jit", self.restriction)), (j, i))
class IndexNester(IdentityMapper):
def map_subscript(self, expr):
if len(expr.index) > 1:
return self.rec(prim.Subscript(prim.Subscript(expr.aggregate, expr.index[:1]), expr.index[1:]))
else:
return expr
def ufl_to_code(expr, measure):
from dune.codegen.generation import global_context
with global_context(integral_type="cell", form_identifier="foo"):
visitor = AdHocVisitor(measure)
nester = IndexNester()
from pymbolic.mapper.c_code import CCodeMapper
ccm = CCodeMapper()
exprs = visitor(expr)
exprs = [nester(e) for e in exprs]
return [ccm(e) for e in exprs]
......@@ -6,8 +6,8 @@ grid:
- 4.0
- 1.0
N:
- 20
- 20
- 40
- 10
refinement: 0
material:
-
......@@ -35,7 +35,7 @@ reinforced_operator:
- -0.5
- 0.749
youngs_modulus: 10000
stabilization_parameter: 1000
stabilization_parameter: 1.0
geneo:
subdomains:
- 4
......
......@@ -182,6 +182,7 @@ int main(int argc, char** argv)
}
bool sym = is_symmetric(istlA,"A",true,1e-12);
std::cout << "matrix is symmetric: " << sym << std::endl;
std::cout << "matrix norm: " << istlA.infinity_norm() << std::endl;
// set up ISTL solver
int itmax = 100000;
......@@ -196,7 +197,7 @@ int main(int argc, char** argv)
ISTLV& istlv = Dune::PDELab::Backend::native(v);
v = 1.0;
set_constrained_dofs(cc,0.0,v); // dirichlet dofs are zero, others are 1, so we can use it as a mask!
solver.apply(istlv,istlr,stat);
solver.apply(istlv,istlr,stat);
istldisplacement -= istlv;
//==================
......
Markdown is supported
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