Merge branch 'feature_quadrature' into 'master'

Feature quadrature

See merge request !2

code/include/option_structure.h

+/*!
+ * \file option_structure.h
+ * \brief All global defined (physical) constants, enums etc
+ * \author <blank>
+ * \version 0.0
+ *
+ */
+
+#ifndef OPTION_STRUCTURE_H
+#define OPTION_STRUCTURE_H
+
+#include <map>
+
+// Definition for global constants goes here
+
+
+// Definition of enums goes here
+/*! \brief Enum for all currently available quadratures in rtsn.
+ *         Option enums are written in capital letters with underscores as spaces (e.g option "time integration" has option enum "TIME_INTEGRATION")
+ */
+enum QUAD_NAME{
+    QUAD_MonteCarlo,
+    QUAD_GaussLegendreTensorized,
+    QUAD_LevelSymmetric,
+    QUAD_Lebedev,
+    QUAD_LDFESA
+};
+
+/*! \brief Conversion Map String to enum
+ */
+inline std::map <std::string,QUAD_NAME> Quadrature_Map {
+  { "MONTE_CARLO", QUAD_MonteCarlo },
+  { "GAUSS_LEGENDRE_TENSORIZED", QUAD_GaussLegendreTensorized },
+  { "LEVEL_SYMMETRIC", QUAD_LevelSymmetric },
+  { "LEBEDEV", QUAD_Lebedev },
+  { "LDFESA",QUAD_LDFESA }
+};
+
+#endif // OPTION_STRUCTURE_H

code/include/qdummy.h

@@ -20,11 +20,10 @@ class QDummy : public Quadrature
     QDummy( unsigned order );
     virtual ~QDummy() {}
 
-    std::string ComputeName();
-    unsigned ComputeNq();
-    VectorVector ComputePoints();
-    Vector ComputeWeights();
-    VectorVectorU ComputeConnectivity();
+    void SetName() override;
+    void SetNq() override;
+    void SetPointsAndWeights() override;
+    void SetConnectivity()override;
 };
 
 #endif    // QDUMMY_H

code/include/qgausslegendretensorized.h

@@ -9,11 +9,10 @@ class QGaussLegendreTensorized : public Quadrature
     QGaussLegendreTensorized( unsigned order );
     virtual ~QGaussLegendreTensorized() {}
 
-    std::string ComputeName() override;
-    unsigned ComputeNq() override;
-    VectorVector ComputePoints() override;
-    Vector ComputeWeights() override;
-    VectorVectorU ComputeConnectivity() override;
+    inline void SetName() override { _name =  "Tensorized Gauss-Legendre quadrature."; }
+    inline void SetNq() override   { _nq = pow( GetOrder(), 2 ); }
+    void SetPointsAndWeights() override;
+    void SetConnectivity() override;
 };
 
 #endif    // QGAUSSLEGENDRETENSORIZED_H

code/include/qldfesa.h

+#ifndef QLDFESA_H
+#define QLDFESA_H
+
+#include "qlookupquadrature.h"
+
+class QLDFESA : public QLookupQuadrature
+{
+  public:
+    QLDFESA( unsigned order );
+    virtual ~QLDFESA() {}
+
+    inline void SetName() override { _name =  "LDFESA quadrature"; }
+    void SetAvailOrders() override;
+    void SetDataInfo() override;
+    void SetConnectivity() override;
+};
+
+#endif // QLDFESA_H

code/include/qlebedev.h

+#ifndef QLEBEDEV_H
+#define QLEBEDEV_H
+
+#include "qlookupquadrature.h"
+
+class QLebedev : public QLookupQuadrature
+{
+  public:
+    QLebedev( unsigned order );
+    virtual ~QLebedev() {}
+
+    inline void SetName() override { _name =  "Lebedev quadrature"; }
+    void SetAvailOrders() override;
+    void SetDataInfo() override;
+    void SetConnectivity() override;
+
+};
+
+#endif // QLEBEDEV_H

code/include/qlevelsymmetric.h

+#ifndef QLEVELSYMMETRIC_H
+#define QLEVELSYMMETRIC_H
+
+#include "qlookupquadrature.h"
+
+class QLevelSymmetric : public QLookupQuadrature
+{
+  public:
+    QLevelSymmetric( unsigned order );
+    virtual ~QLevelSymmetric() {}
+
+    inline void SetName() override { _name =  "Level Symmetric quadrature"; }
+    void SetAvailOrders() override;
+    void SetDataInfo() override;
+    void SetConnectivity() override;
+};
+
+#endif    // QLEVELSYMMETRIC_H

code/include/qlookupquadrature.h

+#ifndef QLOOKUPQUADRATURE_H
+#define QLOOKUPQUADRATURE_H
+
+#include "quadrature.h"
+
+class QLookupQuadrature : public Quadrature
+{
+  public:
+    QLookupQuadrature( unsigned order );
+    virtual ~QLookupQuadrature() {}
+
+    //helper
+    inline std::vector<unsigned> getAvailOrders() const {return _availableOrders;}    /*! @returns: Vector with avalable orders for level symmetric quadrature */
+
+  protected:
+    bool CheckOrder();        /*! @brief checks if given order is available for this quadrature rule. */
+    void SetNq() override;    /*! @brief: Assumes, that _order is available in lookup table */
+    void SetPointsAndWeights() override; /*! @brief reads in n_points gridpoints and -weights from given filename. */
+
+    virtual void SetAvailOrders() = 0; /*! @brief Sets vector with avaialbe Orders and corresponding vector with Nq. */
+    virtual void SetDataInfo() = 0;    /*! @brief Sets data path for lookup table */
+
+    std::vector<unsigned> _availableOrders; /*! @brief: Vector with available orders for lookup table */
+    std::vector<unsigned> _nqByOrder; /*! @brief: Vector with number of quadrature points of each listed order */
+    std::string _dataFiles; /*! path to the lookup table */
+    std::string _dataFileSuffix; /*! suffix to data files of this quadrature */
+
+};
+#endif // QLOOKUPQUADRATURE_H

code/include/qmontecarlo.h

@@ -9,11 +9,10 @@ class QMonteCarlo : public Quadrature
     QMonteCarlo( unsigned order );
     virtual ~QMonteCarlo() {}
 
-    std::string ComputeName() override;
-    unsigned ComputeNq() override;
-    VectorVector ComputePoints() override;
-    Vector ComputeWeights() override;
-    VectorVectorU ComputeConnectivity() override;
+    inline void SetName() override { _name =  "Monte Carlo Quadrature."; }
+    void SetNq() override;
+    void SetPointsAndWeights() override;
+    void SetConnectivity()override;
 };
 
 #endif    // QMONTECARLO_H

code/include/quadrature.h

 #ifndef QUADRATURE_H
 #define QUADRATURE_H
 
-#include "typedef.h"
 #include <iostream>
 #include <string>
-using namespace std;
+
+#include "settings.h"
+#include "typedef.h"
+
 class Quadrature
 {
   public:
     Quadrature( unsigned order );
     virtual ~Quadrature(){}
 
-    /*! @brief Gives the name of the quadrature rule
-     *  @returns string name : Name of quadrature rule     */
-    virtual std::string ComputeName()           = 0;
-
-    /*! @brief Computes the number of gridpoints of the quadrature rule
-     *  @returns unsigned nq : number of gridpoints of the quadrature rule     */
-    virtual unsigned ComputeNq()                = 0;
-
-    /*! @brief Computes the a vector (length: nq) of (coordinates of) gridpoints used for the quadrature rule
-     *  @returns VectorVector coordinates : A Vector of coordinates the gridpoints.     */
-    virtual VectorVector ComputePoints()        = 0;
-
-    /*! @brief Computes the a vector (length: nq) of weights for the gridpoints. The indices match the gridpoints VectorVector.
-     *  @returns Vector weights : A Vector of weights of the gridpoints.     */
-    virtual Vector ComputeWeights()             = 0;
-
-    /*! @brief TODO: How is connectivity defined?.
-     *  @returns VectorVectorU connectivity : TODO */
-    virtual VectorVectorU ComputeConnectivity() = 0;
-
     // Aux functions
     void PrintWeights();            /*! @brief prints: Weight vector */
     void PrintPoints();             /*! @brief prints: Point vectorVector */
@@ -40,7 +22,7 @@ class Quadrature
       *  @returns sum of all weights */
     double SumUpWeights();
 
-     /*! @brief computes the weighted sum of a given density function f over all quadrature points.
+     /*! @brief Integrates f(x,y,z) with the quadrature.
       *  @param double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
       *  @returns double result: result of the quadrature rule */
     double Integrate( double( f )( double x0, double x1, double x2 ) );
@@ -50,15 +32,7 @@ class Quadrature
      *  @param: std::string name: Name of the quadrature rule
      *  @param: unsigned order: Order of the quadrature rule
      *  @returns Quadrature* quadrature: returns pointer to instance of the given derived quadrature class */
-    static Quadrature* CreateQuadrature( std::string name, unsigned order );
-
-    // Setter
-    inline void SetName( std::string name ) { _name = name; }           /*! @brief sets: name of the quadrature */
-    inline void SetOrder( unsigned order ) { _order = order; }          /*! @brief sets: order of the quadrature */
-    inline void SetNq( unsigned nq ) { _nq = nq; }                      /*! @brief sets: number of gridpoints of the quadrature */
-    inline void SetPoints( VectorVector points ) { _points = points; }  /*! @brief sets: coordinates of gridpoints of the quadrature */
-    inline void SetWeights( Vector weights ) { _weights = weights; }    /*! @brief sets: weights of gridpoints of the quadrature */
-    inline void SetConnectivity( VectorVectorU connectivity ) { _connectivity = connectivity; } /*! @brief sets: connectivity vector */
+    static Quadrature* CreateQuadrature( QUAD_NAME name, unsigned order );
 
     // Getter
     inline std::string GetName() const { return _name; }      /*! @returns std::string _name:  name of the quadrature */
@@ -69,6 +43,19 @@ class Quadrature
     inline VectorVectorU GetConnectivity() const { return _connectivity; } /*! @returns VectorVectorU _connectivity:  connectivity of gridpoints of the quadrature */
 
   protected:
+
+    // Setter
+    inline void SetOrder( unsigned order ) { _order = order; } /*! @brief sets: order of the quadrature */
+    virtual void SetName()         = 0; /*! @brief Sets: name of the quadrature */
+    virtual void SetNq()           = 0; /*! @brief sets: number of gridpoints of the quadrature */
+    virtual void SetConnectivity() = 0; /*! @brief sets: Connectivity Adjacency Matrix as VektorVektor*/
+
+    /*! @brief Computes the a vector (length: nq) of (coordinates of) gridpoints used for the quadrature rule.
+     *         Computes the a vector (length: nq) of weights for the gridpoints. The indices match the gridpoints VectorVector.
+     *         Sets computed values for _points and _weights. */
+    virtual void SetPointsAndWeights()        = 0;
+
+    // Member variables
     std::string _name;          /*! @brief name of the quadrature */
     unsigned _order;            /*! @brief order of the quadrature */
     unsigned _nq;               /*! @brief number of gridpoints of the quadrature */

code/include/settings.h

@@ -4,6 +4,8 @@
 #include <filesystem>
 #include <string>
 
+#include "option_structure.h"
+
 class Settings
 {
   private:
@@ -20,7 +22,7 @@ class Settings
 
     // Solver
     unsigned _quadOrder;
-    std::string _quadName;
+    QUAD_NAME _quadName;
     double _CFL;
     double _tEnd;
 
@@ -39,7 +41,7 @@ class Settings
 
     // solver Getters
     unsigned GetQuadOrder() const;
-    std::string GetQuadName() const;
+    QUAD_NAME GetQuadName() const;
     double GetCFL() const;
     double GetTEnd() const;
 

code/input/example.toml

@@ -9,5 +9,5 @@ meshFile = "foo"
 [solver]
 CFL = 0.9
 tEnd = 0.3
-quadType = "bar"
+quadType = "MONTE_CARLO"
 quadOrder = 42 

code/src/io.cpp

 #include "io.h"
+#include <iostream>
 
 Settings* ReadInputFile( std::string inputFile ) {
     bool validConfig = true;
@@ -80,7 +81,13 @@ Settings* ReadInputFile( std::string inputFile ) {
         }
         auto quadType = solver->get_as<std::string>( "quadType" );
         if( quadType ) {
-            settings->_quadName = *quadType;
+            std::string quadTypeString = *quadType;
+            try{
+                settings->_quadName  = Quadrature_Map.at(quadTypeString);
+            }catch (const std::exception& e) {
+                std::cout << "Error: <" << quadTypeString << "> is not a feasible input. Please check the config template!\n";
+                exit(EXIT_FAILURE); // Quit RTSN
+            }
         }
         else {
             spdlog::error( "[inputfile] [solver] 'quadType' not set!" );

code/src/main.cpp

@@ -15,8 +15,9 @@ int main( int argc, char** argv ) {
     PrintLogHeader( settings->GetInputFile() );
 
     // playground and demo
-    Quadrature* Q = Quadrature::CreateQuadrature( "montecarlo", 10 );
+    Quadrature* Q = Quadrature::CreateQuadrature( QUAD_LDFESA, 2 );
     Q->PrintWeights();
+    cout << "====================\n";
     std::cout << Q->Integrate( f ) << std::endl;
 
     auto log = spdlog::get( "event" );

code/src/qdummy.cpp

@@ -12,41 +12,28 @@
 #include "qdummy.h"
 
 QDummy::QDummy( unsigned order ) : Quadrature( order ) {
-    SetName( ComputeName() );
-    SetNq( ComputeNq() );
-    SetPoints( ComputePoints() );
-    SetWeights( ComputeWeights() );
-    SetConnectivity( ComputeConnectivity() );
+    SetName();
+    SetNq();
+    SetPointsAndWeights();
+    SetConnectivity();
 }
 
-std::string QDummy::ComputeName() {
-    std::string name = "";
+void QDummy::SetName() {
+    _name = "";
     // @TODO: DUMMYCOMPUTATION
-    return name;
 }
 
-unsigned QDummy::ComputeNq() {
-    unsigned nq = 0;
-    // @TODO: DUMMYCOMPUTATION
-    return nq;
-};
-
-VectorVector QDummy::ComputePoints() {
-    unsigned nq = GetNq();
-    VectorVector points( nq );
-    // @TODO: DUMMYCOMPUTATION
-    return points;
-};
+void QDummy::SetNq() {
+    _nq = 0;
+}
 
-Vector QDummy::ComputeWeights() {
-    unsigned nq = GetNq();
-    Vector weights( nq, 4.0 * M_PI / nq );
+void QDummy::SetPointsAndWeights() {
+    _points = VectorVector( GetNq() );
+    _weights = Vector ( GetNq(), 4.0 * M_PI / GetNq() );
     // @TODO: DUMMYCOMPUTATION
-    return weights;
-};
+}
 
-VectorVectorU QDummy::ComputeConnectivity() {
+void QDummy::SetConnectivity() { //TODO
     VectorVectorU connectivity;
-    // @TODO: DUMMYCOMPUTATION
-    return connectivity;
-};
+    _connectivity = connectivity;
+}

code/src/qgausslegendretensorized.cpp

 #include "qgausslegendretensorized.h"
 
 QGaussLegendreTensorized::QGaussLegendreTensorized( unsigned order ) : Quadrature( order ) {
-    SetName( ComputeName() );
-    SetNq( ComputeNq() );
-    SetPoints( ComputePoints() );
-    SetWeights( ComputeWeights() );
-    SetConnectivity( ComputeConnectivity() );
+    SetName();
+    SetNq();
+    SetPointsAndWeights();
+    SetConnectivity();
 }
 
-std::string QGaussLegendreTensorized::ComputeName() { return "Tensorized Gauss-Legendre quadrature."; }
-
-unsigned QGaussLegendreTensorized::ComputeNq() { return pow( GetOrder(), 2 ); }
-
-VectorVector QGaussLegendreTensorized::ComputePoints() {
+void QGaussLegendreTensorized::SetPointsAndWeights(){ //TODO
+    //Compute Points
     // Nq random points on the sphere.
-    unsigned nq = GetNq();
-    VectorVector points( nq );
-    return points;
-}
+    _points = VectorVector(GetNq());
 
-Vector QGaussLegendreTensorized::ComputeWeights() {
-    // Equal weights
-    unsigned nq = GetNq();
-    Vector weights( nq, 4.0 * M_PI / nq );
-    return weights;
+    //Compute Weights
+    _weights = Vector( GetNq(), 4.0 * M_PI / GetNq());
 }
 
-VectorVectorU QGaussLegendreTensorized::ComputeConnectivity() {
+void QGaussLegendreTensorized::SetConnectivity() { //TODO
     // Not initialized for this quadrature.
     VectorVectorU connectivity;
-    return connectivity;
-};
+    _connectivity = connectivity;
+}

code/src/qldfesa.cpp

+#include "qldfesa.h"
+
+QLDFESA::QLDFESA( unsigned order ) : QLookupQuadrature( order ){
+
+    SetAvailOrders();
+    SetDataInfo();
+
+    SetName();
+    CheckOrder(); //Check if order is available
+    SetNq(); //Set number of quadrature points
+    SetPointsAndWeights();
+    SetConnectivity();
+}
+
+void QLDFESA::SetAvailOrders() {
+    _availableOrders = {1, 2, 3};
+    _nqByOrder = {32, 128, 512};
+}
+
+void QLDFESA::SetDataInfo() {
+    _dataFiles = "../ext/sphericalquadpy/sphericalquadpy/ldfesa/data/";
+    _dataFileSuffix = "_ldfesa.txt";
+}
+
+void QLDFESA::SetConnectivity() { //TODO
+    VectorVectorU connectivity;
+    _connectivity = connectivity;
+}

code/src/qlebedev.cpp

+#include "qlebedev.h"
+
+QLebedev::QLebedev( unsigned order ) : QLookupQuadrature( order ){
+
+    SetAvailOrders();
+    SetDataInfo();
+
+    SetName();
+    CheckOrder(); //Check if order is available
+    SetNq(); //Set number of quadrature points
+    SetPointsAndWeights();
+    SetConnectivity();
+}
+
+void QLebedev::SetAvailOrders()  {
+    _availableOrders = {3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35, 41, 47, 53,
+                        59, 65, 71, 77, 83, 89, 95, 101, 107, 113, 119, 125, 131}; //Available orders in lookuptable
+    _nqByOrder = {6, 14, 26, 38, 50, 74, 86, 110, 146, 170, 194, 230, 266, 302, 350, 434, 590,
+                  770, 974, 1202, 1454, 1730, 2030, 2354, 2702, 3074, 3470, 3890, 4334, 4802, 5294, 5810};
+}
+
+void QLebedev::SetDataInfo() {
+    _dataFiles = "../ext/sphericalquadpy/sphericalquadpy/lebedev/data/";
+   _dataFileSuffix = "_lebedev.txt";
+}
+
+void QLebedev::SetConnectivity() { //TODO
+    VectorVectorU connectivity;
+    _connectivity = connectivity;
+}
+

code/src/qlevelsymmetric.cpp

+#include "qlevelsymmetric.h"
+
+QLevelSymmetric::QLevelSymmetric( unsigned order ) : QLookupQuadrature( order ){
+
+    SetAvailOrders();
+    SetDataInfo();
+
+    SetName();
+    CheckOrder(); //Check if order is available
+    SetNq(); //Set number of quadrature points
+    SetPointsAndWeights();
+    SetConnectivity();
+}
+
+void QLevelSymmetric::SetAvailOrders() {
+    _availableOrders = {2, 4, 6, 8, 10, 12, 14, 16, 18, 20}; //Available orders in lookuptable
+    _nqByOrder = {8, 24, 48, 80, 120, 168, 224, 288, 360, 432};
+}
+
+void QLevelSymmetric::SetDataInfo() {
+    _dataFiles = "../ext/sphericalquadpy/sphericalquadpy/levelsymmetric/data/";
+    _dataFileSuffix = "_levelsym.txt";
+}
+
+void QLevelSymmetric::SetConnectivity() { //TODO
+    VectorVectorU connectivity;
+    _connectivity = connectivity;
+}

code/src/qlookupquadrature.cpp

+#include <fstream>
+#include <sstream>
+#include "qlookupquadrature.h"
+
+QLookupQuadrature::QLookupQuadrature( unsigned order ) : Quadrature( order ){
+}
+
+bool QLookupQuadrature::CheckOrder(){
+   std::vector<unsigned>::iterator it = std::find(_availableOrders.begin(), _availableOrders.end(), _order);
+
+   if (it == _availableOrders.end()){
+        std::cout << "ERROR! Order "<< _order << " not available. (Replace this error message by a proper exeption handler!)" << std::endl; //TODO: throw proper error!
+        exit(1);
+        return false;
+    }
+    std::cout << _name << " with order " << _order << " chosen.\n";
+    return true;
+}
+
+void QLookupQuadrature::SetNq() {
+    //find iterator of current order
+    std::vector<unsigned>::iterator it = std::find(_availableOrders.begin(), _availableOrders.end(), _order);
+    // Get index of element from iterator
+    unsigned index = std::distance(_availableOrders.begin(), it);
+    // Get _nq
+    _nq = _nqByOrder[index];
+}
+
+void QLookupQuadrature::SetPointsAndWeights() {
+
+    _weights.resize(_nq);
+    _points.resize(_nq);
+
+    std::string filename =  _dataFiles + std::to_string(_order) + _dataFileSuffix;
+    std::string line;
+    unsigned count;
+
+    std::ifstream in(filename); //give file to filestream
+    for(unsigned idx_point = 0; idx_point < _nq ; idx_point++){
+
+        count = 0;
+        _points[idx_point].resize(3);
+
+        in >> line; //Get line of CSV file
+        std::stringstream ss(line); //give line to stringstream
+
+        for (double double_in; ss >> double_in;) { //parse line
+            if(count <3) _points[idx_point][count] = double_in;
+            else _weights[idx_point] = double_in;
+
+            count ++;
+            if (ss.peek() == ',') ss.ignore();
+        }
+    }
+    in.close();
+}

code/src/qmontecarlo.cpp

 #include "qmontecarlo.h"
 
 QMonteCarlo::QMonteCarlo( unsigned order ) : Quadrature( order ) {
-    SetName( ComputeName() );
-    SetNq( ComputeNq() );
-    SetPoints( ComputePoints() );
-    SetWeights( ComputeWeights() );
-    SetConnectivity( ComputeConnectivity() );
+    SetName();
+    SetNq();
+    SetPointsAndWeights();
+    SetConnectivity();
 }
 
-std::string QMonteCarlo::ComputeName() {
-    // Human readable description.
-    return "Monte Carlo Quadrature.";
-}
 
-unsigned QMonteCarlo::ComputeNq() {
+void QMonteCarlo::SetNq() {
     // For Monte Carlo Quadrature, nq = order.
-    return GetOrder();