meancurvaturemodel.hh

#ifndef __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__
#define __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__
 
#include "../../tensors/tensor.hh"
#include "../modelbase.hh"
#include "../expressiontraits.hh"
 
namespace Dune {
 
  namespace ACFem::PDEModel {
 
    using namespace Literals;
 
    template<class FunctionSpace, class Regularization = Tensor::FieldVectorTensor<typename FunctionSpace::RangeFieldType> >
    class MeanCurvatureModel
      : public ModelBase<FunctionSpace>
      , public Expressions::SelfExpression<MeanCurvatureModel<FunctionSpace, Regularization> >
      , public MPL::UniqueTags<ConditionalType<ExpressionTraits<Regularization>::isVolatile, VolatileExpression, void>,
                               ConditionalType<IsConstantExprArg<Regularization>::value, ConstantExpression, void>,
                               ConditionalType<IsTypedValue<Regularization>::value, TypedValueExpression, void>,
                               SemiPositiveExpression, SymmetricModel>
    {
      using ThisType = MeanCurvatureModel;
      using BaseType = ModelBase<FunctionSpace>;
     public:
      using RegularizationType = Regularization;
 
      using typename BaseType::RangeFieldType;
      using typename BaseType::RangeType;
      using typename BaseType::JacobianRangeType;
      using typename BaseType::HessianRangeType;
      using BaseType::dimDomain;
 
      template<class RegularizationArg,
               std::enable_if_t<(std::is_constructible<RegularizationType, RegularizationArg>::value
                                 && std::is_constructible<RegularizationArg, decltype(1_f)>::value
        ), int> = 0>
      MeanCurvatureModel(RegularizationArg&& eps = 1_f, const std::string& name = "")
        : regularization_(std::forward<RegularizationArg>(eps)),
          name_(name == "" ? "MC" : name)
      {}
 
      template<class RegularizationArg,
               std::enable_if_t<(std::is_constructible<RegularizationType, RegularizationArg>::value
                                 && !std::is_constructible<RegularizationArg, decltype(1_f)>::value
        ), int> = 0>
      MeanCurvatureModel(RegularizationArg&& eps, const std::string& name = "")
        : regularization_(std::forward<RegularizationArg>(eps)),
          name_(name == "" ? "MC" : name)
      {}
 
      std::string name() const
      {
        return name_;
      }
 
      JacobianRangeType flux(const JacobianRangeType& jacobian) const
      {
        auto eps = regularization_;
        return jacobian / sqrt(eps*eps + jacobian.frobenius_norm2());
      }
 
      JacobianRangeType linearizedFlux(const JacobianRangeType& DuBar,
                                       const JacobianRangeType& jacobian) const
      {
        RangeFieldType eps = regularization_;
        auto factor = sqrt(eps*eps + DuBar.frobenius_norm2());
        return jacobian / factor - DuBar * contractInner<2>(DuBar, jacobian) / factor / factor / factor;
      }
 
      RangeType fluxDivergence(const JacobianRangeType& jacobian,
                               const HessianRangeType& hessian) const
      {
        auto eps = regularization_;
        auto factor = std::sqrt(eps*eps + jacobian.frobenius_norm2());
 
        // This is then the trace of the Hessian
        RangeType result = -trace(hessian[0]);
 
        // + the correction from the non-linearity
        RangeFieldType correction = 0;
        for (int i = 0; i < dimDomain; ++i) {
          correction += jacobian[0][i] * hessian[0][i][i] * jacobian[0][i];
          for (int j = i+1; j < dimDomain; ++j) {
            correction += 2.0*jacobian[0][i] * hessian[0][i][j] * jacobian[0][j];
          }
        }
        result += correction/factor/factor;
        return result /= factor;
      }
 
     protected:
      RegularizationType regularization_;
      std::string name_;
    };
 
 
    template<class Object, class Regularization = Tensor::FieldVectorTensor<typename Object::FunctionSpaceType::RangeFieldType> >
    auto
    meanCurvatureModel(Regularization&& regularization,
                       const Object& object,
                       const std::string& name = "")
    {
      return expressionClosure(MeanCurvatureModel<typename Object::FunctionSpaceType, Regularization>(std::forward<Regularization>(regularization), name));
    }
 
 
 
 
  } // namespace ACFem::PDEModel
 
  namespace ACFem {
 
    using PDEModel::meanCurvatureModel;
 
  }
 
}  //Namespace Dune
 
 
#endif // __DUNE_ACFEM_MODELS_MODULES_MEANCURVATUREMODEL_HH__