newton_cotes_closed_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 NewtonCotesClosedDriver :
        public DriverBase
      {
      public:
        static constexpr bool variadic = true;
        static constexpr int min_points = 2;
        static constexpr int max_points = 7;
 
        static String name()
        {
          return "newton-cotes-closed";
        }
 
        template<typename Functor_>
        static void alias(Functor_& functor)
        {
          functor.alias("simpson", 3);
          functor.alias("pulcherrima", 4);
          functor.alias("milne-boole", 5);
          functor.alias("6-point", 6);
          functor.alias("weddle", 7);
        }
 
        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 2:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = Coord_(1);
 
            rule.get_weight(0) = Weight_(1);
            rule.get_weight(1) = Weight_(1);
            break;
 
          case 3:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = Coord_(0);
            rule.get_coord(2) = Coord_(1);
 
            rule.get_weight(0) = Weight_(1) / Weight_(3);
            rule.get_weight(1) = Weight_(4) / Weight_(3);
            rule.get_weight(2) = Weight_(1) / Weight_(3);
            break;
 
          case 4:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = -Coord_(1) / Coord_(3);
            rule.get_coord(2) =  Coord_(1) / Coord_(3);
            rule.get_coord(3) =  Coord_(1);
 
            rule.get_weight(0) = Weight_(1) / Weight_(4);
            rule.get_weight(1) = Weight_(3) / Weight_(4);
            rule.get_weight(2) = Weight_(3) / Weight_(4);
            rule.get_weight(3) = Weight_(1) / Weight_(4);
            break;
 
          case 5:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = -Coord_(1) / Coord_(2);
            rule.get_coord(2) =  Coord_(0);
            rule.get_coord(3) =  Coord_(1) / Coord_(2);
            rule.get_coord(4) =  Coord_(1);
 
            rule.get_weight(0) = Weight_(7)  / Weight_(45);
            rule.get_weight(1) = Weight_(32) / Weight_(45);
            rule.get_weight(2) = Weight_(12) / Weight_(45);
            rule.get_weight(3) = Weight_(32) / Weight_(45);
            rule.get_weight(4) = Weight_(7)  / Weight_(45);
            break;
 
          case 6:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = -Coord_(3) / Coord_(5);
            rule.get_coord(2) = -Coord_(1) / Coord_(5);
            rule.get_coord(3) =  Coord_(1) / Coord_(5);
            rule.get_coord(4) =  Coord_(3) / Coord_(5);
            rule.get_coord(5) =  Coord_(1);
 
            rule.get_weight(0) = Weight_(19) / Weight_(144);
            rule.get_weight(1) = Weight_(75) / Weight_(144);
            rule.get_weight(2) = Weight_(50) / Weight_(144);
            rule.get_weight(3) = Weight_(50) / Weight_(144);
            rule.get_weight(4) = Weight_(75) / Weight_(144);
            rule.get_weight(5) = Weight_(19) / Weight_(144);
            break;
 
          case 7:
            rule.get_coord(0) = -Coord_(1);
            rule.get_coord(1) = -Coord_(2) / Coord_(3);
            rule.get_coord(2) = -Coord_(1) / Coord_(3);
            rule.get_coord(3) =  Coord_(0);
            rule.get_coord(4) =  Coord_(1) / Coord_(3);
            rule.get_coord(5) =  Coord_(2) / Coord_(3);
            rule.get_coord(6) =  Coord_(1);
 
            rule.get_weight(0) = Weight_(41)  / Weight_(420);
            rule.get_weight(1) = Weight_(216) / Weight_(420);
            rule.get_weight(2) = Weight_(27)  / Weight_(420);
            rule.get_weight(3) = Weight_(272) / Weight_(420);
            rule.get_weight(4) = Weight_(27)  / Weight_(420);
            rule.get_weight(5) = Weight_(216) / Weight_(420);
            rule.get_weight(6) = Weight_(41)  / Weight_(420);
            break;
          }
        }
      }; // class NewtonCotesClosedDriver<...>
    } // namespace Scalar
  } // namespace Cubature
} // namespace FEAT