Commit 1f80bf0c authored by Dominic Kempf's avatar Dominic Kempf
Browse files

WIP on correct the flipping for symmetry

parent e0c4c4e2
......@@ -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
......@@ -352,6 +352,11 @@ class EulerBernoulli2DLocalOperator
// The tangential vector for the curve
auto t = fibre->tangent(it->second);
auto n = ig.intersection().centerUnitOuterNormal();
if (( t*n < 0) != flipped) {
std::cout << "FLIPPED IS WRONG" << std::endl;
}
// Extract physical parameter of the fibre
auto E = fibre_modulus[fibindex];
auto d = fibre_radii[fibindex];
......@@ -380,12 +385,25 @@ 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*( penalty * sk_dtun * dtvn_n_0));
// r_n.accumulate(lfsu_n.child(1), i, E*I*( penalty * sk_dtun * dtvn_n_1));
// r_s.accumulate(lfsu_s.child(0), i, -E*I*( penalty * sk_dtun * dtvn_s_0));
// r_s.accumulate(lfsu_s.child(1), i, -E*I*( penalty * sk_dtun * dtvn_s_1));
r_n.accumulate(lfsu_n.child(0), i, E*I*(- flip_factor * sk_dt2un*dtvn_n_0 - flip_factor * dt2vn_n_0 * sk_dtun));
r_n.accumulate(lfsu_n.child(1), i, E*I*(- flip_factor * sk_dt2un*dtvn_n_1 - flip_factor * dt2vn_n_1 * sk_dtun));
r_s.accumulate(lfsu_s.child(0), i, -E*I*(- flip_factor * sk_dt2un*dtvn_s_0 - flip_factor * dt2vn_s_0 * sk_dtun));
r_s.accumulate(lfsu_s.child(1), i, -E*I*(- flip_factor * sk_dt2un*dtvn_s_1 - flip_factor * dt2vn_s_1 * sk_dtun));
// r_n.accumulate(lfsu_n.child(0), i, E*I*(- flip_factor * sk_dt2un*dtvn_n_0 - flip_factor * 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 * dt2vn_n_1 * sk_dtun + penalty * sk_dtun * dtvn_n_1));
// r_s.accumulate(lfsu_s.child(0), i, -E*I*(- flip_factor * sk_dt2un*dtvn_s_0 - flip_factor * 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 * dt2vn_s_1 * sk_dtun + penalty * sk_dtun * dtvn_s_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));
}
}
}
......
......@@ -58,19 +58,19 @@ class FibreReinforcedBulkOperator
virtual void alpha_volume(const EG& eg, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const override
{
bulkoperator.alpha_volume(eg, lfsu, x, lfsv, r);
fibreoperator.alpha_volume(eg, lfsu, x, lfsv, r);
//bulkoperator.alpha_volume(eg, lfsu, x, lfsv, r);
//fibreoperator.alpha_volume(eg, lfsu, x, lfsv, r);
}
virtual void alpha_boundary(const IG& ig, const LFSU& lfsu, const X& x, const LFSV& lfsv, R& r) const override
{
bulkoperator.alpha_boundary(ig, lfsu, x, lfsv, r);
//bulkoperator.alpha_boundary(ig, lfsu, x, lfsv, r);
}
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:
......
......@@ -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