root/src/io_l2h1Calculator.cpp

DEFINITIONS

1. assign_v1
2. assign_v2

   1 /*! \file io_l2h1Calculator.cpp
   2     \brief Implementation of io_l2h1Calculator.h
   3 
   4      Revision history
   5      --------------------------------------------------
   6 
   7           Revised by Christian Power April 2016
   8           Originally written by Christian Power
   9                (power22c@gmail.com) Januar 2016
  10 
  11      \author Christian Power
  12      \date 27. April 2016
  13      \copyright Copyright (c) 2016 Christian Power.  All rights reserved.
  14  */
  15 
  16 #include <config.h>
  17 #include <stdexcept>
  18 #include <dune/fem/misc/l2norm.hh>
  19 #include <dune/fem/misc/h1norm.hh>
  20 #include <dassert.h>
  21 #include "io_l2h1Calculator.h"
  22 #include "grid.h"
  23 
  24 
  25 using namespace std;
  26 using Esfem::Io::L2H1_calculator;
  27 
  28 //! \f$f\colon \R^3 \to \R\f$
  29 using FEfun = Esfem::Grid::Scal_FEfun::Dune_FEfun;
  30 //! \f$f\colon \R^3 \to \R^3\f$
  31 using Vec_FEfun = Esfem::Grid::Vec_FEfun::Dune_FEfun;
  32 //! \f$L^2\f$-norm
  33 using L2_norm = Dune::Fem::L2Norm<Esfem::Grid::Grid_and_time::Grid_part>;
  34 //! \f$H^1\f$-norm
  35 using H1_norm = Dune::Fem::H1Norm<Esfem::Grid::Grid_and_time::Grid_part>;
  36 
  37 // ----------------------------------------------------------------------
  38 // Actual implementations 
  39 
  40 //! %Data members of L2H1_calculator
  41 struct L2H1_calculator::Data{
  42   //! Dune \f$L^2\f$-norm functor
  43   L2_norm l2;
  44   //! Dune \f$H^1\f$-norm functor
  45   H1_norm h1;
  46   //! First vector valued finite element function
  47   /*! __Invariant:__ `size() == worlddim`*/
  48   vector<FEfun> vec_fefun1;
  49   //! Second vector valued finite element function
  50   /*! \copydetails vec_fefun1 */
  51   vector<FEfun> vec_fefun2;
  52   //! Get grid
  53   /*! \post Grid_and_time must outlive this object. 
  54    \note I could add try and catch block, to avoid cryptic 
  55    error messages. */
  56   Data(const Grid::Grid_and_time& gt);
  57   //! Change value of vec_fefun1
  58   /*! \copydetails assign() */
  59   void assign_v1(const Vec_FEfun& v){ assign(v, vec_fefun1); }
  60   //! Change_value of vec_fefun2
  61   /*! \copydetails assign() */
  62   void assign_v2(const Vec_FEfun& v){ assign(v, vec_fefun2); }
  63 private:
  64   //! Constructor helper
  65   void init_vec(vector<FEfun>& v, const Grid::Grid_and_time::FE_space& fes);
  66   //! Actual implementation of assign_v1() and assign_v2()
  67   /*! \post Dimension of Range should be `world_dim()`
  68     and the underlying grid should be the same. */
  69   void assign(const Vec_FEfun& v, vector<FEfun>& vec);
  70 };
  71 
  72 L2H1_calculator::Data::Data(const Grid::Grid_and_time& gt) 
  73 :l2 {gt.grid_part()}, h1{gt.grid_part()} {
  74   init_vec(vec_fefun1, gt.fe_space());
  75   init_vec(vec_fefun2, gt.fe_space());
  76 }
  77 void L2H1_calculator::Data::init_vec
  78 (vector<FEfun>& v, const Grid::Grid_and_time::FE_space& fes){
  79   constexpr auto d = Grid::world_dim();
  80   v.reserve(d);
  81   for(int it = 0; it < d; ++it){
  82     v.emplace_back("vec_fefun", fes);
  83     v.back().clear();
  84   }
  85 }
  86 void L2H1_calculator::Data::assign
  87 (const Vec_FEfun& v, vector<FEfun>& vec){
  88   constexpr auto d = Grid::world_dim();
  89   size_t no_nodes = v.size() / d; // assuming the post condition
  90   Assert::dynamic(v.size() == no_nodes * d, 
  91                   Assert::compose(__FILE__, __LINE__, "assign()"));
  92   for(int jt = 0; jt < d; ++jt){
  93     auto vfef_ptr = v.dbegin();
  94     auto fef_ptr = vec[jt].dbegin();
  95     vfef_ptr += jt;
  96     for(size_t it = 0; it < no_nodes; ++it, vfef_ptr += d, ++fef_ptr)
  97       *fef_ptr = *vfef_ptr;
  98   }
  99 }
 100 
 101 L2H1_calculator::L2H1_calculator(const Grid::Grid_and_time& gt)
 102   :d_ptr {make_unique<Data>(gt)} {}
 103 Esfem::Io::L2H1_calculator::~L2H1_calculator() = default;
 104 
 105 double L2H1_calculator::l2_err
 106 (const Grid::Scal_FEfun& u, const Grid::Scal_FEfun& uN) const{
 107   const FEfun& u1 = u;
 108   const FEfun& u2 = uN;
 109   return d_ptr->l2.distance(u1, u2);
 110 }
 111 double L2H1_calculator::l2_err
 112 (const Grid::Vec_FEfun& u, const Grid::Vec_FEfun& uN) const{
 113   d_ptr->assign_v1(u);
 114   d_ptr->assign_v2(uN);
 115   auto rv = 0.;
 116   for(size_t it = 0; it < d_ptr->vec_fefun1.size(); ++it){
 117     rv += d_ptr->
 118       l2.distance(d_ptr->vec_fefun1[it], d_ptr->vec_fefun2[it]);
 119   }
 120   return rv;
 121   /*
 122   const Vec_FEfun& u1 = u;
 123   const Vec_FEfun& u2 = uN;
 124   return d_ptr->l2.distance(u1,u2);
 125   */
 126 }
 127 double L2H1_calculator::h1_err
 128 (const Grid::Scal_FEfun& u, const Grid::Scal_FEfun& uN) const{
 129   const FEfun& u1 = u;
 130   const FEfun& u2 = uN;
 131   return d_ptr->h1.distance(u1, u2);
 132 }
 133 double L2H1_calculator::h1_err
 134 (const Grid::Vec_FEfun& u, const Grid::Vec_FEfun& uN) const{
 135   d_ptr->assign_v1(u);
 136   d_ptr->assign_v2(uN);
 137   auto rv = 0.;
 138   for(size_t it = 0; it < d_ptr->vec_fefun1.size(); ++it){
 139     rv += d_ptr->
 140       h1.distance(d_ptr->vec_fefun1[it], d_ptr->vec_fefun2[it]);
 141   }
 142   return rv;
 143   /*
 144   const Vec_FEfun& u1 = u;
 145   const Vec_FEfun& u2 = uN;
 146   return d_ptr->h1.distance(u1,u2);
 147   */
 148 }