root/src/secOrd_op_rhs_u.cpp

DEFINITIONS

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

   1 /*! \file secOrd_op_rhs_u.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) February 2016
  10 
  11      Implementation details for secOrd_op_rhs_u.h
  12      Created by Christian Power on 04.02.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_u.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_u
  34   : public Dune::Fem::Function<Esfem::Grid::Grid_and_time::Function_space, RHS_data_u>
  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_u() = delete;
  42   explicit RHS_data_u(const Dune::Fem::TimeProviderBase&,
  43                     const Esfem::Io::Parameter&);
  44   RHS_data_u(const RHS_data_u&) = delete;
  45   RHS_data_u& operator=(const RHS_data_u&) = 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 a;
  52   const double gamma;
  53 };
  54 
  55 void assemble_RHS(const RHS_data_u&, FE_function&);
  56 void matrixFree_assembly(const Dune::Fem::TimeProviderBase&, const Geometry&,
  57                          const Quadrature&, FE_function::LocalFunctionType&);
  58 void massMatrix_for_entity(const Geometry&, const Quadrature&, const RHS_data_u&,
  59                            FE_function::LocalFunctionType&);
  60 
  61 // ----------------------------------------------------------------------
  62 // Implementation esfem.h
  63 
  64 struct Esfem::SecOrd_op::Rhs_u::Data{
  65   RHS_data_u rhs;
  66   const Dune::Fem::TimeProviderBase& tp;
  67   const Grid::Grid_and_time::FE_space& fe_space;
  68   Data(const Io::Parameter& p, const Grid::Grid_and_time& gt)
  69     : rhs {gt.time_provider(), p}, tp {gt.time_provider()},
  70       fe_space {gt.fe_space()}
  71   {}
  72 };
  73 
  74 Esfem::SecOrd_op::Rhs_u::Rhs_u(const Io::Parameter& p, const Grid::Grid_and_time& gt){
  75   try{
  76     d_ptr = new Data {p, gt};
  77   }
  78   catch(std::exception&){
  79     std::throw_with_nested(std::logic_error
  80                            {"Error in constructor of Grid_and_time."});
  81   }
  82   catch(...){
  83     throw std::logic_error {"Unknown error in constructor of Grid_and_time."};
  84   }
  85 }
  86 Esfem::SecOrd_op::Rhs_u::~Rhs_u(){
  87   delete d_ptr;
  88   d_ptr = nullptr;
  89 #ifdef DEBUG
  90   std::cerr << "~Rhs_u(): delete d_ptr.\n";
  91 #endif
  92 }
  93 void Esfem::SecOrd_op::Rhs_u::assemble_and_addScaled_to(Grid::Scal_FEfun& fef) const{
  94   static FE_function tmp {"local_tmp", d_ptr -> fe_space};
  95   assemble_RHS(d_ptr -> rhs, tmp);
  96   FE_function& dune_fef = fef;
  97   dune_fef.axpy(d_ptr -> tp.deltaT(), tmp); 
  98 }
  99 void Esfem::SecOrd_op::Rhs_u::assemble(Grid::Scal_FEfun& fef) const{
 100   FE_function& dune_fef = fef;
 101   assemble_RHS(d_ptr -> rhs, dune_fef);
 102 }
 103 
 104 // ----------------------------------------------------------------------
 105 // Internal implementation
 106 
 107 RHS_data_u::RHS_data_u(const Dune::Fem::TimeProviderBase& tpb,
 108                    const Esfem::Io::Parameter& p)
 109   : tp {tpb}, a {p.tg_a()}, gamma {p.tg_gamma()}
 110 {}
 111 void RHS_data_u::evaluate(const Domain& d, Range& r) const{
 112   static_assert(Domain::dimension == 3, "Bad domain dimension.");
 113   static_assert(Range::dimension == 1, "Bad range dimension.");
 114   const double x = d[0];
 115   const double y = d[1];
 116   const double z = d[2];
 117   const double t = tp.time();
 118   const double g = gamma;
 119   const double u = std::exp(-6. * t) * x * y;
 120   // const double w = std::exp(-6. * t) * y * z;
 121   r = u
 122   //r =
 123     // - g * a + g * u  // rhs for w == 0.
 124     // tp.deltaT() * g * a 
 125     // u
 126     // tp.deltaT() * g * u * u * w
 127     // #include "/Users/christianpower/cpp/syntax/data/brusselator_u_cpp.txt"
 128     // #include "/Users/christianpower/cpp/syntax/data/brusselator_u_v2_cpp.txt"
 129     ;
 130 }
 131 RHS_data_u::Range RHS_data_u::operator()(const Domain& d) const{
 132   Range r {0};
 133   evaluate(d,r);
 134   return r;
 135 }
 136 
 137 void assemble_RHS(const RHS_data_u& rhs_fun, FE_function& fef){
 138   fef.clear();
 139   const auto& df_space = fef.space();
 140   for(const auto& entity : df_space){
 141     const auto& geometry = entity.geometry();
 142     const Quadrature quad {entity, 2 * df_space.order() + 1};
 143     auto fef_local = fef.localFunction(entity);
 144     massMatrix_for_entity(geometry, quad, rhs_fun, fef_local);
 145   }
 146   fef.communicate();    
 147 }
 148 
 149 void massMatrix_for_entity(const Geometry& g, const Quadrature& q,
 150                            const RHS_data_u& rhs_fun,
 151                            FE_function::LocalFunctionType& f_loc){
 152   for(std::size_t pt = 0; pt < q.nop(); ++pt){
 153     const auto& x = q.point(pt);
 154     RHS_data_u::Range fx {rhs_fun(g.global(x))};
 155     fx *= q.weight(pt) * g.integrationElement(x);
 156     f_loc.axpy(q[pt], fx);
 157   }  
 158 }
 159 
 160 
 161 /*! Log:
 162  */