newton_cotes_open_driver.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/cubature/scalar/driver_base.hpp>
 
namespace FEAT
{
  namespace Cubature
  {
    namespace Scalar
    {
      class NewtonCotesOpenDriver :
        public DriverBase
      {
      public:
        static constexpr bool variadic = true;
        static constexpr int min_points = 1;
        static constexpr int max_points = 7;
 
        static String name()
        {
          return "newton-cotes-open";
        }
 
        template<
          typename Weight_,
          typename Coord_>
        static void fill(Rule<Weight_, Coord_>& rule, int num_points)
        {
          // how many points do we have?
          switch(num_points)
          {
          case 1:
            rule.get_coord(0) = Coord_(0);
 
            rule.get_weight(0) = Weight_(2);
            break;
 
          case 2:
            rule.get_coord(0) = -Coord_(1) / Coord_(3);
            rule.get_coord(1) =  Coord_(1) / Coord_(3);
 
            rule.get_weight(0) = Weight_(1);
            rule.get_weight(1) = Weight_(1);
            break;
 
          case 3:
            rule.get_coord(0) = -Coord_(1) / Coord_(2);
            rule.get_coord(1) = Coord_(0);
            rule.get_coord(2) = Coord_(1) / Coord_(2);
 
            rule.get_weight(0) = Weight_(4) / Weight_(3);
            rule.get_weight(1) = - Weight_(2) / Weight_(3);
            rule.get_weight(2) = Weight_(4) / Weight_(3);
            break;
 
          case 4:
            rule.get_coord(0) = -Coord_(3) / Coord_(5);
            rule.get_coord(1) = -Coord_(1) / Coord_(5);
            rule.get_coord(2) =  Coord_(1) / Coord_(5);
            rule.get_coord(3) =  Coord_(3) / Coord_(5);
 
            rule.get_weight(0) = Weight_(11) / Weight_(12);
            rule.get_weight(1) = Weight_(1) / Weight_(12);
            rule.get_weight(2) = Weight_(1) / Weight_(12);
            rule.get_weight(3) = Weight_(11) / Weight_(12);
            break;
 
          case 5:
            rule.get_coord(0) = -Coord_(2) / Coord_(3);
            rule.get_coord(1) = -Coord_(1) / Coord_(3);
            rule.get_coord(2) =  Coord_(0);
            rule.get_coord(3) =  Coord_(1) / Coord_(3);
            rule.get_coord(4) =  Coord_(2) / Coord_(3);
 
            rule.get_weight(0) = Weight_(11) / Weight_(10);
            rule.get_weight(1) = - Weight_(14) / Weight_(10);
            rule.get_weight(2) = Weight_(26) / Weight_(10);
            rule.get_weight(3) = - Weight_(14) / Weight_(10);
            rule.get_weight(4) = Weight_(11) / Weight_(10);
            break;
 
          case 6:
            rule.get_coord(0) = -Coord_(5) / Coord_(7);
            rule.get_coord(1) = -Coord_(3) / Coord_(7);
            rule.get_coord(2) = -Coord_(1) / Coord_(7);
            rule.get_coord(3) =  Coord_(1) / Coord_(7);
            rule.get_coord(4) =  Coord_(3) / Coord_(7);
            rule.get_coord(5) =  Coord_(5) / Coord_(7);
 
            rule.get_weight(0) =  Weight_(611) / Weight_(720);
            rule.get_weight(1) = -Weight_(453) / Weight_(720);
            rule.get_weight(2) =  Weight_(562) / Weight_(720);
            rule.get_weight(3) =  Weight_(562) / Weight_(720);
            rule.get_weight(4) = -Weight_(453) / Weight_(720);
            rule.get_weight(5) =  Weight_(611) / Weight_(720);
            break;
 
          case 7:
            rule.get_coord(0) = -Coord_(3) / Coord_(4);
            rule.get_coord(1) = -Coord_(1) / Coord_(2);
            rule.get_coord(2) = -Coord_(1) / Coord_(4);
            rule.get_coord(3) =  Coord_(0);
            rule.get_coord(4) =  Coord_(1) / Coord_(4);
            rule.get_coord(5) =  Coord_(1) / Coord_(2);
            rule.get_coord(6) =  Coord_(3) / Coord_(4);
 
            rule.get_weight(0) =  Weight_(920)  / Weight_(945);
            rule.get_weight(1) = -Weight_(1908) / Weight_(945);
            rule.get_weight(2) =  Weight_(4392) / Weight_(945);
            rule.get_weight(3) = -Weight_(4918) / Weight_(945);
            rule.get_weight(4) =  Weight_(4392) / Weight_(945);
            rule.get_weight(5) = -Weight_(1908) / Weight_(945);
            rule.get_weight(6) =  Weight_(920)  / Weight_(945);
            break;
          }
        }
      }; // class NewtonCotesOpenDriver<...>
    } // namespace Scalar
  } // namespace Cubature
} // namespace FEAT