root/src/secOrd_op_solutionDriven_impl.cpp

DEFINITIONS

1. brusselator_rhs
2. mcf_lhs_matrixFree_assembly
3. mcf_rhs_matrixFree_assembly
4. surface_normal
5. oriented_basis

   1 /*! \file secOrd_op_solutionDriven_impl.cpp
   2     \brief Implementation of `secOrd_op_solutionDriven_impl.h`
   3 
   4     Revision history
   5     --------------------------------------------------
   6 
   7          Revised by Christian Power April 2016
   8          Revised by Christian Power March 2016
   9          Originally written by Christian Power
  10               (power22c@gmail.com) March 2016
  11 
  12     Idea
  13     --------------------------------------------------
  14 
  15     Implementation of helper classes of `secOrd_op_solutionDriven.h`.
  16 
  17 
  18     \author Christian Power
  19     \date 26. April 2016
  20     \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  21 */
  22 
  23 #include <dassert.h>
  24 #include "esfem_error.h"
  25 #include "secOrd_op_solutionDriven_impl.h"
  26 
  27 using std::size_t;
  28 using Esfem::Impl::MCF_op; 
  29 using Vector_fef = MCF_op::Vector_fef;
  30 using Scalar_fef = MCF_op::Scalar_fef;
  31 //! Vector valued finite element space
  32 using Vec_fe_space = Vector_fef::DiscreteFunctionSpaceType;
  33 template<typename T>
  34 using Local_function = typename T::LocalFunctionType;
  35 template<typename T>
  36 using Domain = typename Local_function<T>::DomainType;
  37 template<typename T>
  38 using Range = typename Local_function<T>::RangeType;
  39 using Geometry = MCF_op::Geometry;
  40 using Grid_part = MCF_op::Grid_part;
  41 using Quadrature = MCF_op::Quadrature;
  42 template<typename T>
  43 using Jacobian_range = typename Local_function<T>::JacobianRangeType;
  44 // Forward typedef
  45 // using Linear_operator = MCF_op::Linear_operator;
  46 
  47 
  48 
  49 // static void prepare_linearOperator_matrix(const Vec_fe_space&, Linear_operator&);
  50 // static std::size_t calculate_matrix_row_size(const Vector_fef&);
  51 
  52 // ----------------------------------------------------------------------
  53 // Implementation of Esfem::Impl::MCF_op
  54 
  55 // ------------------------------------------------------------
  56 // Public interface
  57 
  58 MCF_op::MCF_op(const Io::Parameter& p,
  59                const Grid::Grid_and_time& g,
  60                const Scalar_fef& u_input)
  61   :alpha {p.velocity_regularization()},
  62    delta {p.surface_growthFactor()},
  63    epsilon {p.mcf_regularization()},
  64    eps {p.eps()},
  65    tp {g.time_provider()},
  66    u {u_input}
  67 {}
  68 
  69 void MCF_op::operator()(const Vector_fef& rhs, Vector_fef& lhs) const{
  70   // (M + (alpha + epsilon * tau) A) X
  71   lhs.clear();
  72   const auto& df_space = lhs.space();
  73   for(const auto& entity : df_space){
  74     const auto& geometry = entity.geometry();
  75     const auto rhs_loc = rhs.localFunction(entity);
  76     auto lhs_loc = lhs.localFunction(entity);
  77     Quadrature quad {entity, rhs_loc.order() + lhs_loc.order()};
  78     mcf_lhs_matrixFree_assembly(geometry, quad, rhs_loc, lhs_loc);
  79   }
  80   lhs.communicate();
  81 }
  82 
  83 void MCF_op::brusselator_rhs(const Vector_fef& rhs, Vector_fef& lhs) const{
  84   // (M + alpha * A) X + tau * delta * M(u^n, surfaceNormal)
  85   lhs.clear();
  86   const auto& df_space = lhs.space();
  87   for(const auto& entity : df_space){
  88     const auto& geometry = entity.geometry();
  89     const auto rhs_loc = rhs.localFunction(entity);
  90     const auto u_loc = u.localFunction(entity); // Scalar valued
  91     auto lhs_loc = lhs.localFunction(entity);
  92     Quadrature quad {entity, rhs_loc.order() + lhs_loc.order()};
  93     mcf_rhs_matrixFree_assembly(geometry, quad, rhs_loc, u_loc, lhs_loc);
  94   }
  95   lhs.communicate();
  96 }
  97 
  98 // void MCF_op::jacobian_matrix(const Vector_fef& fef, Linear_operator& matrix) const{
  99 //   const auto& df_space = fef.space();
 100 //   prepare_linearOperator_matrix(df_space, matrix);
 101 //   for(const auto& entity : df_space){
 102 //     const auto fef_loc = fef.localFunction(entity);
 103 //     auto matrix_loc = matrix.localMatrix(entity, entity);
 104 //     // jacobian_matrix_heat(entity, fef_loc, matrix_loc);
 105 //     // jacobian_matrix_quadMass(entity, fef_loc, matrix_loc);
 106 //   }
 107 //   matrix.communicate();
 108 // }
 109 
 110 // ------------------------------------------------------------
 111 // Private interface
 112 
 113 void MCF_op::mcf_lhs_matrixFree_assembly(const Geometry& g,
 114                                          const Quadrature& q,
 115                                          const Local_function<Vector_fef>& cf,
 116                                          Local_function<Vector_fef>& f) const{
 117   for(size_t pt = 0; pt < q.nop(); ++pt){
 118     // (M + (alpha + epsilon * tau) A) X
 119     const auto& x = q.point(pt);
 120     const auto integral_factor = q.weight(pt) * g.integrationElement(x);
 121 
 122     auto X_p = evaluate(pt, q, cf);
 123     X_p *= integral_factor;
 124     
 125     auto dX_p = jacobian(pt, q, cf);
 126     dX_p *= integral_factor * (alpha + epsilon * tp.deltaT());
 127 
 128     f.axpy(q[pt], X_p, dX_p);
 129   }
 130 }
 131 void MCF_op::mcf_rhs_matrixFree_assembly(const Geometry& g,
 132                                          const Quadrature& q,
 133                                          const Local_function<Vector_fef>& cf,
 134                                          const Local_function<Scalar_fef>& u_loc,
 135                                          Local_function<Vector_fef>& f) const{
 136   for(size_t pt = 0; pt < q.nop(); ++pt){
 137     // (M + alpha * A) X + tau * delta * M(u^n, surfaceNormal)
 138     const auto& x = q.point(pt);
 139     const auto integral_factor = q.weight(pt) * g.integrationElement(x);
 140 
 141     auto n_p = surface_normal(g);
 142       
 143     const auto u_p = evaluate(pt, q, u_loc);
 144     n_p *= u_p * tp.deltaT() * delta * integral_factor;
 145 
 146     auto X_p = evaluate(pt, q, cf);
 147     X_p *= integral_factor;
 148     X_p += n_p;
 149     
 150     auto dX_p = jacobian(pt, q, cf);
 151     dX_p *= alpha * integral_factor;
 152 
 153     f.axpy(q[pt], X_p, dX_p);   
 154   }
 155 }
 156 
 157 MCF_op::Range<Vector_fef>
 158 MCF_op::surface_normal(const Geometry& g) const{
 159   static_assert(dim_vec_domain == 3, "Bad dimension");
 160   const auto basis = oriented_basis(g);
 161   auto normal = nonUnit_normal(basis);
 162   const auto norm = euclidean_norm(normal);
 163   Assert::dynamic<Assert::level(7), Esfem::SolutionDriven_error>
 164     (norm > eps,
 165      Assert::compose(__FILE__, __LINE__, "Norm of normal vector almost vanishes."));
 166   normal /= norm;
 167   return normal;
 168 }
 169 
 170 // ----------------------------------------------------------------------
 171 // Implementation helper function
 172 
 173 std::vector<MCF_op::Domain<Vector_fef> >
 174 Esfem::Impl::oriented_basis(const MCF_op::Geometry& g){
 175   const Domain<Vector_fef> p0 = g.corner(0);
 176   const Domain<Vector_fef> p1 = g.corner(1);
 177   const Domain<Vector_fef> p2 = g.corner(2);
 178 
 179   // Oriented basis of tangent space
 180   // Draw a picture to understand this.
 181   const Domain<Vector_fef> v = p2 - p0;
 182   const Domain<Vector_fef> w = p1 - p0;
 183   return {v, w};
 184 }
 185 
 186 MCF_op::Range<Vector_fef> Esfem::Impl::nonUnit_normal
 187 (const std::vector<MCF_op::Domain<Vector_fef> >& basis){
 188   Assert::dynamic<Assert::level(7), Esfem::SolutionDriven_error>
 189     (basis.size() == 2, "Dimension of basis bad.");
 190   const auto& v = basis[0];
 191   const auto& w = basis[1];
 192 
 193   // Cross product formula
 194   Range<Vector_fef> normal;
 195   normal[0] = v[1] * w[2] - v[2] * w[1];
 196   normal[1] = - v[0] * w[2] + v[2] * w[0];
 197   normal[2] = v[0] * w[1] - v[1] * w[0];
 198 
 199   return normal;
 200 }
 201 
 202 // ----------------------------------------------------------------------
 203 // Internal Implementation
 204 
 205 // std::size_t calculate_matrix_row_size(const Vector_fef& fef){
 206 //   const auto& df_space = fef.space();
 207 //   return df_space.blockMapper().maxNumDofs() * df_space.localBlockSize;
 208 //   // localBlockSize is equal to 1 for scalar functions
 209 // }
 210 
 211 // void prepare_linearOperator_matrix(const Vec_fe_space& df_space,
 212 //                                 Linear_operator& matrix){
 213 //   Dune::Fem::DiagonalStencil<Vec_fe_space, Vec_fe_space>
 214 //     stencil {df_space, df_space};
 215 //   matrix.reserve(stencil);
 216 //   matrix.clear();
 217 // }