referenceneumannfunction.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ELASTICITY_UTILITIES_REFERENCE_NEUMANN_FUNCTION_HH
#define DUNE_ELASTICITY_UTILITIES_REFERENCE_NEUMANN_FUNCTION_HH

#include <dune/common/fvector.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/hybridutilities.hh>

#include <dune/fufem/functions/basisgridfunction.hh>
#include <dune/matrix-vector/transpose.hh>
#include <dune/matrix-vector/addtodiagonal.hh>

/** A basis grid function that transforms Neumann boundary values, defined in deformed coordinates, back to reference coordinates.
 *
 *  The formular is proven in Ciarlet, Mathematical Elasticity Volume I: Three-Dimensional Elasticity, Thm 1.7-1 and Section 2.6
 */
template <class Basis, class CoefficientVector>
class ReferenceNeumannFunction :
    public BasisGridFunction<Basis, CoefficientVector>
{
  using Base = BasisGridFunction<Basis, CoefficientVector>;
  using Grid = typename Basis::GridView::Grid;
  using field_type = Dune::field_t<CoefficientVector>;

  static constexpr int dim = Grid::dimension;

public:
  using LocalDomainType = typename Base::LocalDomainType;
  using RangeType = typename Base::RangeType;
  using Element = typename Base::Element;
  using LocalFaceDomainType = typename Grid::template Codim<1>::LocalGeometry::LocalCoordinate;

  ReferenceNeumannFunction(const CoefficientVector& displacementCoefficient, const Basis& basis, const CoefficientVector& neumannValues) :
    Base(basis, neumannValues),
    displacementField_(basis, displacementCoefficient)
  {}

  void evaluateLocal(const Element& e, const LocalDomainType& x, RangeType& y) const override
  {

    typename Base::DerivativeType deformationGradient, invTransposedDefGradient;
    displacementField_.evaluateDerivativeLocal(e, x, deformationGradient);
    Dune::MatrixVector::addToDiagonal(deformationGradient, 1.0);

    deformationGradient.invert();
    Dune::MatrixVector::transpose(deformationGradient, invTransposedDefGradient);

    Dune::FieldVector<field_type, dim> normal;
    for (const auto& is : intersections(this->basis_.getGridView(), e)) {
      auto localCoordinates = is.geometryInInside().local(x);
      if (intersectionContainsVertex(localCoordinates)) {
        normal = is.unitOuterNormal(localCoordinates);
        break;
      }
    }

    LocalEvaluator<Basis, CoefficientVector, RangeType>::evaluateLocal(this->basis_, this->coefficients_, e, x, y);

    Dune::FieldVector<field_type, dim> transformedNormal(0);
    invTransposedDefGradient.mv(normal, transformedNormal);
    y *= deformationGradient.determinant()*transformedNormal.two_norm();
  }
private:

  // can't this be done nicer?
  bool intersectionContainsVertex(const LocalFaceDomainType& localPoint) const {
    field_type eps(1e-14);
    bool contains;
    Dune::Hybrid::ifElse((dim == 3),
      [&](auto&&) {contains = (-eps <= localPoint[0]) and (localPoint[0] <= 1 + eps) and
                        (-eps <= localPoint[1]) and (localPoint[1] <= 1 + eps) and
                        (localPoint[0] + localPoint[1] <= 1 + eps);},
      [&](auto&&) {contains = (-eps <= localPoint[0]) and (localPoint[0] <= 1 + eps);}
    );
    return contains;
  }

  Base displacementField_;
};


#endif