root/src/secOrd_op_solutionDriven_impl.h

/* [<][>][^][v][top][bottom][index][help] */

INCLUDED FROM


DEFINITIONS

This source file includes following definitions.
  1. euclidean_norm
  2. evaluate
  3. evaluate
  4. jacobian

   1 /*! \file secOrd_op_solutionDriven_impl.h
   2     \author Christian Power
   3     \date 18. March 2016
   4 
   5     \brief Helper classes for `Esfem::SecOrd_op::Solution_driven`
   6 
   7      Revision history
   8      --------------------------------------------------
   9 
  10           Revised by Christian Power March 2016
  11           Originally written by Christian Power
  12                (power22c@gmail.com) March 2016
  13 
  14      Idea
  15      --------------------------------------------------
  16 
  17      Implementation of `secOrd_op_solutionDriven.h`
  18 
  19          Created by Christian Power on 18.03.2016
  20          Copyright (c) 2016 Christian Power.  All rights reserved.
  21 
  22 */
  23 
  24 #ifndef SECORD_OP_SOLUTIONDRIVEN_IMPL_H
  25 #define SECORD_OP_SOLUTIONDRIVEN_IMPL_H
  26 
  27 #include <cmath>
  28 #include <vector>
  29 #include <config.h>
  30 #include <dune/fem/operator/common/operator.hh>
  31 #include <dune/fem/operator/linear/istloperator.hh>
  32 #include <dune/fem/solver/istlsolver.hh>
  33 #include <dune/fem/operator/common/differentiableoperator.hh>
  34 #include <dune/fem/quadrature/cachingquadrature.hh>
  35 #include "io_parameter.h"
  36 #include "grid.h"
  37 
  38 namespace Esfem{
  39   namespace Impl{
  40     class MCF_op 
  41       : public Dune::Fem::Operator<Esfem::Grid::Vec_FEfun::Dune_FEfun>
  42     {
  43     public:
  44       using Scalar_fef = Esfem::Grid::Scal_FEfun::Dune_FEfun;
  45       /*!< \brief \f$ f\colon \R^3 \to \R \f$ */
  46       using Vector_fef = Esfem::Grid::Vec_FEfun::Dune_FEfun;
  47       /*!< \brief \f$ f\colon \R^3 \to \R^3 \f$ */
  48       // Big Matrix
  49       // using Linear_operator
  50       // = Dune::Fem::ISTLLinearOperator<Vector_fef, Vector_fef>;
  51       // Sub matrix for an element 
  52       // using Local_matrix = Linear_operator::LocalMatrixType;
  53       template<typename T>
  54       using Local_function = typename T::LocalFunctionType;
  55       /*!< \brief Finite element function restricted to a triangle. */
  56       template<typename T>
  57       using Domain = typename Local_function<T>::DomainType;
  58       template<typename T>
  59       using Range = typename Local_function<T>::RangeType;
  60       using RangeField = typename Local_function<Vector_fef>::RangeFieldType;
  61       /*!< \brief Simply \f$ \R \f$ */
  62 
  63       MCF_op(const Io::Parameter&,
  64              const Grid::Grid_and_time&,
  65              const Scalar_fef& u_input);
  66       /*!< \sa Constructor of `Esfem::SecOrd_op::Solution_driven`
  67        */
  68       MCF_op(const MCF_op&) = delete;
  69       MCF_op& operator=(const MCF_op&) = delete;
  70       
  71       void operator()(const Vector_fef& rhs, Vector_fef& lhs) const override;
  72       /*!< \brief Applying the mean curvature operator.  
  73                   The cg solver uses this.
  74 
  75         The precise formulation reads as
  76         \f{equation*}{
  77           \parentheses[\big]{M_3^n + (\alpha + \varepsilon\tau) A_3^n} 
  78           \nodalValue{X}^{n+1} 
  79         \f}
  80        */
  81       void brusselator_rhs(const Vector_fef& rhs, Vector_fef& lhs) const;
  82       /*!< \brief Generates rhs for the linear system.
  83 
  84         The new value of the finite element function will be
  85         \f{equation*}{
  86           (M_3^n + \alpha A_3^n) \nodalValue{X}^n 
  87           + \tau \delta M_3^n(\nodalValue{u}^n, 
  88           \nodalValue{\surfaceNormal}).
  89         \f}
  90        */
  91 
  92       // Assemble matrix for linear system
  93       // void jacobian_matrix(const Vector_fef&, Linear_operator&) const;
  94       
  95       /*! \name Variable containing dimensions
  96        */
  97       //@{
  98       static constexpr int dim_vec_domain = Domain<Vector_fef>::dimension;
  99       static constexpr int dim_vec_range = Range<Vector_fef>::dimension;
 100       static constexpr int dim_scalar_domain = Domain<Scalar_fef>::dimension;
 101       static constexpr int dim_scalar_range = Range<Scalar_fef>::dimension;
 102       //@}
 103 
 104       // ------------------------------------------------------------
 105       // Typedef for helper functions
 106       using Geometry 
 107       = Vector_fef::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
 108       /*!< \brief Geometry means access to the finite element. */
 109       using Grid_part = Vector_fef::GridPartType;
 110       /*!< \brief Auxiliary template argument */
 111       using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
 112       template<typename T>
 113       using Jacobian_range = typename Local_function<T>::JacobianRangeType;
 114       /*!< \brief Jacobian of the local finite element function
 115         \warning `Jacobian_range` is used like this: `jac[0][1]`.
 116       */
 117       
 118       static_assert(dim_vec_domain == dim_vec_range, "Bad dimension");
 119       static_assert(dim_scalar_range == 1, "Bad dimension");
 120     private:
 121       void mcf_lhs_matrixFree_assembly(const Geometry&,
 122                                        const Quadrature&,
 123                                        const Local_function<Vector_fef>&,
 124                                        Local_function<Vector_fef>&) const;
 125       /*!< \brief Used by operator()
 126         
 127         The matrix vector formulation reads as
 128         \f{equation*}{
 129           \parentheses[\big]{M + (\alpha + \varepsilon \tau) A}
 130           \nodalValue{X}
 131         \f}
 132        */
 133       void mcf_rhs_matrixFree_assembly(const Geometry&,
 134                                        const Quadrature&,
 135                                        const Local_function<Vector_fef>&,
 136                                        const Local_function<Scalar_fef>& u_loc,
 137                                        Local_function<Vector_fef>&) const;
 138       /*!< \brief Used by rhs()
 139 
 140         The matrix vector formulation reads as
 141         \f{equation*}{
 142           (M + \alpha A) \nodalValue{X} + \tau \delta 
 143           M(\nodalValue{u}, \surfaceNormal)
 144         \f}
 145        */      
 146       Range<Vector_fef> surface_normal(const Geometry&) const;
 147       /*!< \brief Calculates the outwards pointing normal vector. */
 148       
 149       const double alpha; 
 150       /*!< \brief Velocity Laplace regularization parameter by Lubich */
 151       const double delta; /*!< \brief Tumor growth parameter */
 152       const double epsilon; 
 153       /*!< \brief Mean curvature regularization parameter by Elliott */
 154       const double eps; /*!< \brief Generic tolerance */
 155       const Dune::Fem::TimeProviderBase& tp; /*!< \brief Time step provider */
 156       const Scalar_fef& u; 
 157       /*!< \brief Concentration of the growth promoting chemical substance */
 158     };
 159     /*!< \sa `Esfem::SecOrd_op::Solution_driven` */
 160 
 161     std::vector<MCF_op::Domain<MCF_op::Vector_fef> >
 162     oriented_basis(const MCF_op::Geometry&);
 163     /*!< \brief Calculates an oriented basis for the tangent space.
 164                 Assumes `grid dimension == 2` and
 165                 `world dimension == 3`.
 166          \return Two (numerical) vectors which represent an
 167                  oriented basis of the tangent space.
 168     */
 169     MCF_op::Range<MCF_op::Vector_fef> nonUnit_normal
 170     (const std::vector<MCF_op::Domain<MCF_op::Vector_fef> >& basis);
 171     /*!< \brief Calculates an non-normalized normal vector via
 172                 the cross product formula.
 173          \warning Assumes `basis.size() == 2`
 174                    and that the basis is correctly oriented.
 175          \return Outward pointing non unit normal vector.
 176      */
 177     inline MCF_op::RangeField
 178     euclidean_norm(const MCF_op::Range<MCF_op::Vector_fef>& v){
 179       static_assert(MCF_op::dim_vec_range == 3, "Bad dimension");
 180       return std::sqrt(v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
 181     }
 182     /*!< \brief Euclidean norm for a 3-dimensional vector. */
 183 
 184     /*! \name Evaluation helper functions
 185         \param pt Quadrature point obtained by `MCF_op::Quadrature`.
 186         \param q The quadrature object itself
 187         \param cf The local finit element function
 188 
 189         Evaluates the local finite element function
 190         on the quadrature points.  I have to pass the point and
 191         the quadrature
 192         object itself since I do not know what its precise type is.
 193         I do not consider using template a good practise for this
 194         problem.
 195 
 196         \todo Figure out the type of `q[pt]`.
 197      */
 198     //@{
 199     inline MCF_op::Range<MCF_op::Vector_fef>
 200     evaluate(const std::size_t pt,
 201              const MCF_op::Quadrature& q,
 202              const MCF_op::Local_function<MCF_op::Vector_fef>& cf){
 203       MCF_op::Range<MCF_op::Vector_fef> X;
 204       cf.evaluate(q[pt], X);
 205       return X;
 206     }
 207     /*! \brief Returns `X(p)` */
 208     inline MCF_op::Range<MCF_op::Scalar_fef>
 209     evaluate(const std::size_t pt,
 210              const MCF_op::Quadrature& q,
 211              const MCF_op::Local_function<MCF_op::Scalar_fef>& cf){
 212       MCF_op::Range<MCF_op::Scalar_fef> u;
 213       cf.evaluate(q[pt], u);
 214       return u;
 215     }
 216     /*!< \brief Returns `u(p)` */
 217     inline MCF_op::Jacobian_range<MCF_op::Vector_fef>
 218     jacobian(const std::size_t pt,
 219              const MCF_op::Quadrature& q,
 220              const MCF_op::Local_function<MCF_op::Vector_fef>& cf){
 221       MCF_op::Jacobian_range<MCF_op::Vector_fef> nabla_X;
 222       cf.jacobian(q[pt], nabla_X);
 223       return nabla_X;
 224     }
 225     /*!< \brief Returns `dX(p)`*/
 226     //@}
 227 
 228   }     // namespace Impl
 229 }       // namespace Esfem
 230 
 231 #endif // SECORD_OP_SOLUTIONDRIVEN_IMPL_H

/* [<][>][^][v][top][bottom][index][help] */