root/src/secOrd_op_rhs_w.cpp

DEFINITIONS

1. assemble_and_addScaled_to
2. evaluate
3. assemble_RHS
4. massMatrix_for_entity

   1 /*! \file secOrd_op_rhs_w.cpp
   2 
   3     \brief <Program Name>
   4 
   5      Revision history:
   6 
   7           Revised by Christian Power dd.mm.yyyy
   8           Originally written by Christian Power
   9                (power22c@gmail.com) 30. Januar 2016
  10 
  11      Implementation details for secOrd_op_rhs_w.h
  12      Created by Christian Power on 30.01.2016
  13      Copyright (c) 2016 Christian Power. All rights reserved.
  14  */
  15 
  16 #include <cmath>
  17 #include <config.h>
  18 #include <dune/fem/quadrature/cachingquadrature.hh>
  19 #include "secOrd_op_rhs_w.h"
  20 #include "grid.h"
  21 #include "io_parameter.h"
  22 
  23 #ifdef DEBUG
  24 #include <iostream>
  25 #endif
  26 
  27 using FE_function = Esfem::Grid::Scal_FEfun::Dune_FEfun;
  28 using Geometry
  29 = FE_function::DiscreteFunctionSpaceType::IteratorType::Entity::Geometry;
  30 using Grid_part = FE_function::GridPartType;
  31 using Quadrature = Dune::Fem::CachingQuadrature<Grid_part, 0>;
  32 
  33 class RHS_data_w
  34   : public Dune::Fem::Function<Esfem::Grid::Grid_and_time::Function_space, RHS_data_w>
  35 {
  36 public:
  37   using Base = Esfem::Grid::Grid_and_time::Function_space;
  38   using Domain = Base::DomainType;
  39   using Range = Base::RangeType;
  40   
  41   RHS_data_w() = delete;
  42   explicit RHS_data_w(const Dune::Fem::TimeProviderBase&,
  43                     const Esfem::Io::Parameter&);
  44   RHS_data_w(const RHS_data_w&) = delete;
  45   RHS_data_w& operator=(const RHS_data_w&) = delete;
  46 
  47   void evaluate(const Domain&, Range&) const;
  48   Range operator()(const Domain&) const;
  49 private:
  50   const Dune::Fem::TimeProviderBase& tp;
  51   const double d;
  52   const double b;
  53   const double gamma;
  54 };
  55 
  56 void assemble_RHS(const RHS_data_w&, FE_function&);
  57 void matrixFree_assembly(const Dune::Fem::TimeProviderBase&, const Geometry&,
  58                          const Quadrature&, FE_function::LocalFunctionType&);
  59 void massMatrix_for_entity(const Geometry&, const Quadrature&, const RHS_data_w&,
  60                            FE_function::LocalFunctionType&);
  61 
  62 // ----------------------------------------------------------------------
  63 // Implementation esfem.h
  64 
  65 struct Esfem::SecOrd_op::Rhs_w::Data{
  66   RHS_data_w rhs;
  67   const Dune::Fem::TimeProviderBase& tp;
  68   const Grid::Grid_and_time::FE_space& fe_space;
  69   Data(const Io::Parameter& p, const Grid::Grid_and_time& gt)
  70     : rhs {gt.time_provider(), p}, tp {gt.time_provider()},
  71       fe_space {gt.fe_space()}
  72   {}
  73 };
  74 
  75 Esfem::SecOrd_op::Rhs_w::Rhs_w(const Io::Parameter& p,
  76                                const Grid::Grid_and_time& gt){
  77   try{
  78     d_ptr = new Data {p, gt};
  79   }
  80   catch(std::exception&){
  81     std::throw_with_nested(std::logic_error
  82                            {"Error in constructor of Grid_and_time."});
  83   }
  84   catch(...){
  85     throw std::logic_error {"Unknown error in constructor of Grid_and_time."};
  86   }
  87 }
  88 Esfem::SecOrd_op::Rhs_w::~Rhs_w(){
  89   delete d_ptr;
  90   d_ptr = nullptr;
  91 #ifdef DEBUG
  92   std::cerr << "~Rhs_w(): delete d_ptr.\n";
  93 #endif
  94 }
  95 void Esfem::SecOrd_op::Rhs_w::assemble_and_addScaled_to(Grid::Scal_FEfun& fef) const{
  96   static FE_function tmp {"local_tmp", d_ptr -> fe_space};
  97   assemble_RHS(d_ptr -> rhs, tmp);
  98   FE_function& dune_fef = fef;
  99   dune_fef.axpy(d_ptr -> tp.deltaT(), tmp); 
 100   // dune_fef.assign(tmp);      // for testing
 101 }
 102 
 103 // ----------------------------------------------------------------------
 104 // Internal implementation
 105 
 106 RHS_data_w::RHS_data_w(const Dune::Fem::TimeProviderBase& tpb,
 107                    const Esfem::Io::Parameter& p)
 108   : tp {tpb}, d {p.tg_Dc()}, b {p.tg_b()}, gamma {p.tg_gamma()}
 109 {}
 110 void RHS_data_w::evaluate(const Domain& dom, Range& r) const{
 111   static_assert(Domain::dimension == 3, "Bad domain dimension.");
 112   static_assert(Range::dimension == 1, "Bad range dimension.");
 113   const double x = dom[0];
 114   const double y = dom[1];
 115   const double z = dom[2];
 116   const double t = tp.time();
 117   const double g = gamma;
 118   r =
 119     - g * b     // for w == 0
 120     // #include "/Users/christianpower/cpp/syntax/data/brusselator_w_cpp.txt"
 121     // #include "/Users/christianpower/cpp/syntax/data/brusselator_w_v2_cpp.txt"
 122     ;
 123 }
 124 RHS_data_w::Range RHS_data_w::operator()(const Domain& d) const{
 125   Range r {0};
 126   evaluate(d,r);
 127   return r;
 128 }
 129 
 130 void assemble_RHS(const RHS_data_w& rhs_fun, FE_function& fef){
 131   fef.clear();
 132   const auto& df_space = fef.space();
 133   for(const auto& entity : df_space){
 134     const auto& geometry = entity.geometry();
 135     const Quadrature quad {entity, 2 * df_space.order() + 1};
 136     auto fef_local = fef.localFunction(entity);
 137     massMatrix_for_entity(geometry, quad, rhs_fun, fef_local);
 138   }
 139   fef.communicate();    
 140 }
 141 
 142 void massMatrix_for_entity(const Geometry& g, const Quadrature& q,
 143                            const RHS_data_w& rhs_fun,
 144                            FE_function::LocalFunctionType& f_loc){
 145   for(std::size_t pt = 0; pt < q.nop(); ++pt){
 146     const auto& x = q.point(pt);
 147     RHS_data_w::Range fx {rhs_fun(g.global(x))};
 148     // rhs_fun.evaluate(g.global(x), fx);
 149     fx *= q.weight(pt) * g.integrationElement(x);
 150     f_loc.axpy(q[pt], fx);
 151   }  
 152 }
 153 
 154 
 155 /*! Log:
 156  */