Merge branch 'mini_sprint' of https://git.scc.kit.edu/rtsn/rtsn into entropy_closure

README.md

@@ -77,7 +77,7 @@ tar xzf VTK-8.2.0.tar.gz
 mkdir VTK-build
 cd VTK-build
 cmake -DCMAKE_BUILD_TYPE=Release -DBUILD_DOCUMENTATION=OFF -DBUILD_TESTING=OFF -DCMAKE_INSTALL_PREFIX=~/VTK-install ../VTK-8.2.0
-make -j10
+make -j
 make install
 cd -
 rm -r VTK-8.2.0 VTK-build
@@ -95,10 +95,6 @@ Append `HINTS VTK_INSTALL_DIR` to the `find_package( VTK ... )` line in the CMak
 ```bash
 find_package( VTK REQUIRED COMPONENTS vtkIOGeometry vtkFiltersCore HINTS ~/VTK-install )
 ```
-As a temporary fix, also set the `BLAZE_BLAS_MODE` preprocessor definition to `0`
-```bash
-add_compile_definitions( BLAZE_BLAS_MODE=0 )
-```
 
 Compile it
 ```bash
@@ -107,7 +103,7 @@ module load compiler/gnu/9.2
 module load mpi/openmpi/4.0
 cd code/build/release/
 cmake -DCMAKE_BUILD_TYPE=Release ../../
-make -j10
+make -j
 ```
 
 ---

code/CMakeLists.txt

@@ -3,9 +3,9 @@ project( RTSN VERSION 0.1.0 LANGUAGES CXX )
 set( CMAKE_CXX_STANDARD 17 )
 set( CMAKE_CXX_STANDARD_REQUIRED ON )
 
-set( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -march=native -Wno-dev" CACHE STRING "Flags used by the compiler during release builds." FORCE )
-set( CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -march=native -pg -no-pie" CACHE STRING "Flags used by the compiler during release builds with debug information." FORCE )
-set( CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -Wall -Werror" CACHE STRING "Flags used by the compiler during debug builds." FORCE )
+set( CMAKE_CXX_FLAGS_RELEASE "${CMAKE_CXX_FLAGS_RELEASE} -march=native -Wno-dev" )
+set( CMAKE_CXX_FLAGS_RELWITHDEBINFO "${CMAKE_CXX_FLAGS_RELWITHDEBINFO} -march=native -pg -no-pie" )
+set( CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -O0 -Wall -Werror" )
 
 set( CMAKE_RUNTIME_OUTPUT_DIRECTORY "${CMAKE_SOURCE_DIR}/bin" )
 
@@ -25,8 +25,18 @@ include_directories( ${MPI_INCLUDE_PATH} )
 find_package( LAPACK REQUIRED )
 include_directories( ${LAPACK_INCLUDE_DIR} )
 
+find_package( BLAS )
+if( BLAS_FOUND )
+    message( STATUS "Blaze: BLAS mode enabled" )
+    add_compile_definitions( BLAZE_BLAS_MODE=1 )
+    include_directories( ${BLAS_INCLUDE_DIR} )
+else()
+    message( STATUS "Blaze: BLAS mode disabled" )
+    add_compile_definitions( BLAZE_BLAS_MODE=0 )
+endif()
+
 add_compile_definitions( BLAZE_USE_SHARED_MEMORY_PARALLELIZATION=0 )
-add_compile_definitions( BLAZE_BLAS_MODE=1 )
+
 set(flag 1)
 if (WIN32)
     execute_process(COMMAND wmic cpu get L3CacheSize
@@ -97,7 +107,7 @@ endif (flag)
 string( REGEX MATCH "([0-9][0-9]+)" tmp ${tmp} )
 math( EXPR BLAZE_CACHE_SIZE ${tmp}*1024 )
 add_compile_definitions( BLAZE_CACHE_SIZE=${BLAZE_CACHE_SIZE}UL )
-message( "-- Blaze: Automatic Cache Size Configuration = ${BLAZE_CACHE_SIZE} KiB" )
+message( STATUS "Blaze: Automatic Cache Size Configuration = ${BLAZE_CACHE_SIZE} KiB" )
 include_directories( ${CMAKE_SOURCE_DIR}/ext/blaze )
 
 add_compile_definitions( METIS_EXPORT= )
@@ -120,7 +130,7 @@ execute_process(
 )
 add_compile_definitions( GIT_HASH="${GIT_HASH}" )
 
-set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${MPI_LIBRARIES} ${VTK_LIBRARIES} parmetis -lstdc++fs )
+set( CORE_LIBRARIES ${LAPACK_LIBRARIES} ${BLAS_LIBRARIES} ${MPI_LIBRARIES} ${VTK_LIBRARIES} parmetis -lstdc++fs )
 
 target_link_libraries( ${CMAKE_PROJECT_NAME} ${CORE_LIBRARIES} )
 

parmetis @ 54affdf3

-Subproject commit e68c8e3947f9383e28cb8623d116ca996b3651a9
+Subproject commit 54affdf3f9e1d55c0047c0401401da308aea8878

code/include/mesh.h

@@ -6,14 +6,14 @@
 #include <omp.h>
 #include <vector>
 
-#include "blaze/Blaze.h"
+#include "blaze/math/CompressedMatrix.h"
 #include "metis.h"
 #include "parmetis.h"
 #include "spdlog/spdlog.h"
 
-#include "typedef.h"
-
 #include "settings/globalconstants.h"
+#include "toolboxes/errormessages.h"
+#include "typedef.h"
 
 class Mesh
 {
@@ -52,18 +52,57 @@ class Mesh
           std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries );
     ~Mesh();
 
-    unsigned GetDim() const;
-    unsigned GetNumCells() const;
-    unsigned GetNumNodes() const;
-    unsigned GetNumNodesPerCell() const;
+    inline unsigned GetDim() const { return _dim; }
+    inline unsigned GetNumCells() const { return _numCells; }
+    inline unsigned GetNumNodes() const { return _numNodes; }
+    inline unsigned GetNumNodesPerCell() const { return _numNodesPerCell; }
+
+    /**
+     * @brief Returns all node coordinates
+     * @return dimension: numNodes x dim
+     */
     const std::vector<Vector>& GetNodes() const;
+
+    /**
+     * @brief  Returns the mid point coordinates of each cell
+     * @return dimension: numCells x dim
+     */
     const std::vector<Vector>& GetCellMidPoints() const;
+
+    /**
+     * @brief Returns all node IDs that construct up each cell
+     * @return dimension: numCells x numNodes
+     */
     const std::vector<std::vector<unsigned>>& GetCells() const;
+
+    /**
+     * @brief Returns the cell area of each cell
+     * @return dimension: numCells
+     */
     const std::vector<double>& GetCellAreas() const;
+
+    /**
+     * @brief Return the color/ID of the mesh partition
+     * @return dimension: numCells
+     */
     const std::vector<unsigned>& GetPartitionIDs() const;
+
+    /**
+     * @brief Returns the neighbor cell IDs for every cell
+     * @return dimension: numCells x numNodes
+     */
     const std::vector<std::vector<unsigned>>& GetNeighbours() const;
+
+    /**
+     * @brief Returns the edge length scaled normal vectors of each cell
+     * @return dimension: numCells x numNodes x dim
+     */
     const std::vector<std::vector<Vector>>& GetNormals() const;
 
+    /**
+     * @brief Returns the boundary enum for each cell. BOUNDARY_TYPE::NONE is the default.
+     * @return dimension: numCells
+     */
     const std::vector<BOUNDARY_TYPE>& GetBoundaryTypes() const;
 
     double GetDistanceToOrigin( unsigned idx_cell ) const;

code/include/quadratures/qgausslegendretensorized.h

@@ -5,6 +5,10 @@
 
 class QGaussLegendreTensorized : public QuadratureBase
 {
+  private:
+    double Pythag( const double a, const double b );
+    std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
+
   public:
     QGaussLegendreTensorized( unsigned order );
     virtual ~QGaussLegendreTensorized() {}

code/include/quadratures/quadraturebase.h

@@ -2,6 +2,7 @@
 #define QUADRATURE_H
 
 #include "settings/globalconstants.h"
+#include "toolboxes/errormessages.h"
 #include "typedef.h"
 #include <iostream>
 #include <string>

code/include/settings/config.h

@@ -21,93 +21,93 @@
 #include "optionstructure.h"
 
 /*!
- * \class Config
- * \brief Main class for defining the problem; basically this class reads the configuration file, and
+ * @class Config
+ * @brief Main class for defining the problem; basically this class reads the configuration file, and
  *        stores all the information.
  */
 
 class Config
 {
   private:
-    std::string _fileName; /*!< \brief Name of the current file without extension */
+    std::string _fileName; /*!< @rief Name of the current file without extension */
     bool _baseConfig;
 
-    int _commRank, _commSize; /*!< \brief MPI rank and size.*/    // Not yet used!!
+    int _commRank, _commSize; /*!< @brief MPI rank and size.*/    // Not yet used!!
 
     // --- Options ---
     // File Structure
-    std::string _inputDir;   /*!< \brief Directory for input files*/
-    std::string _outputDir;  /*!< \brief Directory for output files*/
-    std::string _outputFile; /*!< \brief Name of output file*/
-    std::string _logDir;     /*!< \brief Directory of log file*/
-    std::string _meshFile;   /*!< \brief Name of mesh file*/
+    std::string _inputDir;   /*!< @brief Directory for input files*/
+    std::string _outputDir;  /*!< @brief Directory for output files*/
+    std::string _outputFile; /*!< @brief Name of output file*/
+    std::string _logDir;     /*!< @brief Directory of log file*/
+    std::string _meshFile;   /*!< @brief Name of mesh file*/
 
     // Quadrature
-    QUAD_NAME _quadName;       /*!< \brief Quadrature Name*/
-    unsigned short _quadOrder; /*!< \brief Quadrature Order*/
+    QUAD_NAME _quadName;       /*!< @brief Quadrature Name*/
+    unsigned short _quadOrder; /*!< @brief Quadrature Order*/
     unsigned _nQuadPoints;
 
     // Mesh
     unsigned _nCells;
 
     // Solver
-    double _CFL;  /*!< \brief CFL Number for Solver*/
-    double _tEnd; /*!< \brief Final Time for Simulation */
+    double _CFL;  /*!< @brief CFL Number for Solver*/
+    double _tEnd; /*!< @brief Final Time for Simulation */
     PROBLEM_NAME _problemName;
     SOLVER_NAME _solverName;
 
     bool _cleanFluxMat;
 
     // Boundary Conditions
-    /*!< \brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
+    /*!< @brief List of all Pairs (marker, BOUNDARY_TYPE), e.g. (farfield,DIRICHLET).
          Each Boundary Conditions must have an entry in enum BOUNDARY_TYPE*/
     std::vector<std::pair<std::string, BOUNDARY_TYPE>> _boundaries;
-    unsigned short _nMarkerDirichlet;          /*!< \brief Number of Dirichlet BC markers. Enum entry: DIRICHLET */
-    unsigned short _nMarkerNeumann;            /*!< \brief Number of Neumann BC markers. Enum entry: Neumann */
-    std::vector<std::string> _MarkerDirichlet; /*!< \brief Dirichlet BC markers. */
-    std::vector<std::string> _MarkerNeumann;   /*!< \brief Neumann BC markers. */
+    unsigned short _nMarkerDirichlet;          /*!< @brief Number of Dirichlet BC markers. Enum entry: DIRICHLET */
+    unsigned short _nMarkerNeumann;            /*!< @brief Number of Neumann BC markers. Enum entry: Neumann */
+    std::vector<std::string> _MarkerDirichlet; /*!< @brief Dirichlet BC markers. */
+    std::vector<std::string> _MarkerNeumann;   /*!< @brief Neumann BC markers. */
 
     // Scattering Kernel
-    KERNEL_NAME _kernelName; /*!< \brief Scattering Kernel Name*/
+    KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
 
     // --- Parsing Functionality and Initializing of Options ---
     /*!
-     * \brief Set default values for all options not yet set.
+     * @brief Set default values for all options not yet set.
      */
     void SetDefault( void );
 
     /*!
-     * \brief Set the config options.
+     * @brief Set the config options.
      *        ==> Set new config options here.
      */
     void SetConfigOptions( void );
 
     /*!
-     * \brief Set the config file parsing.
+     * @brief Set the config file parsing.
      */
     void SetConfigParsing( char case_filename[MAX_STRING_SIZE] );
 
     /*!
-     * \brief Config file screen output.
+     * @brief Config file screen output.
      */
     void SetOutput( void );
 
     /*!
-     * \brief Initializes pointers to null
+     * @brief Initializes pointers to null
      */
     void SetPointersNull( void );
 
     /*!
-     * \brief Config file postprocessing.
+     * @brief Config file postprocessing.
      */
     void SetPostprocessing( void );
 
     /*!
-     * \brief breaks an input line from the config file into a set of tokens
-     * \param[in] str - the input line string
-     * \param[out] option_name - the name of the option found at the beginning of the line
-     * \param[out] option_value - the tokens found after the "=" sign on the line
-     * \return false if the line is empty or a commment, true otherwise
+     * @brief breaks an input line from the config file into a set of tokens
+     * @param[in] str - the input line string
+     * @param[out] option_name - the name of the option found at the beginning of the line
+     * @param[out] option_value - the tokens found after the "=" sign on the line
+     * @return false if the line is empty or a commment, true otherwise
      */
     bool TokenizeString( std::string& str, std::string& option_name, std::vector<std::string>& option_value );
 
@@ -137,7 +137,7 @@ class Config
     // List and Array options should also be able to be specified with the string "NONE" indicating that there
     // are no elements. This allows the option to be present in a config file but left blank.
 
-    /*!< \brief addDoubleOption creates a config file parser for an option with the given name whose
+    /*!< @brief addDoubleOption creates a config file parser for an option with the given name whose
      value can be represented by a su2double.*/
 
     // Simple Options
@@ -168,19 +168,19 @@ class Config
 
   public:
     /*!
-     * \brief Constructor of the class which reads the input file.
+     * @brief Constructor of the class which reads the input file.
      */
     Config( char case_filename[MAX_STRING_SIZE] );
 
     /*!
-     * \brief Destructor of the class.
+     * @brief Destructor of the class.
      */
     ~Config( void );
 
     // ---- Getters for option values ----
 
     /*!
-     * \brief Get Value of this option.
+     * @brief Get Value of this option.
      *        Please keep alphabetical order within each subcategory
      */
     // File structure
@@ -206,7 +206,7 @@ class Config
     SOLVER_NAME inline GetSolverName() const { return _solverName; }
 
     // Boundary Conditions
-    BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! \brief Get Boundary Type of given marker */
+    BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
 
     // Scattering Kernel
     KERNEL_NAME inline GetKernelName() const { return _kernelName; }

code/include/settings/globalconstants.h

@@ -15,8 +15,8 @@
 
 // --- Definition for global constants goes here ---
 
-const double PI_NUMBER             = M_PI; /*!< \brief Pi number. */
-const unsigned int MAX_STRING_SIZE = 200;  /*!< \brief Maximum size for strings. */
+const double PI_NUMBER             = M_PI; /*!< @brief Pi number. */
+const unsigned int MAX_STRING_SIZE = 200;  /*!< @brief Maximum size for strings. */
 
 // --- Definition of enums goes here ---
 
@@ -24,12 +24,12 @@ enum BOUNDARY_TYPE { DIRICHLET, NEUMANN, NONE, INVALID };
 
 // --- Definition of enums for EnumOptions goes here ---
 
-/*! \brief Enum for all currently available quadratures in rtsn.
+/*! @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
+/*! @brief Conversion Map String to enum
  */
 inline std::map<std::string, QUAD_NAME> Quadrature_Map{ { "MONTE_CARLO", QUAD_MonteCarlo },
                                                         { "GAUSS_LEGENDRE_TENSORIZED", QUAD_GaussLegendreTensorized },

code/include/typedef.h

 #ifndef TYPEDEFS_H
 #define TYPEDEFS_H
 
+#include <blaze/math/CompressedMatrix.h>
 #include <blaze/math/DynamicMatrix.h>
 #include <blaze/math/DynamicVector.h>
 

code/input/checkerboard.cfg

@@ -11,5 +11,5 @@ PROBLEM = CHECKERBOARD
 
 % ---- Boundary Conditions ----
 BC_NEUMANN = ( void )
-QUAD_TYPE = LEBEDEV
-QUAD_ORDER = 35
+QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
+QUAD_ORDER = 10

code/src/io.cpp

@@ -34,7 +34,7 @@ void ExportVTK( const std::string fileName,
                 pts->SetPoint( nodeID++, node[0], node[1], node[2] );
             }
             else {
-                exit( EXIT_FAILURE );
+                ErrorMessages::Error( "Unsupported dimension (d=" + std::to_string( dim ) + ")!", CURRENT_FUNCTION );
             }
         }
         vtkCellArraySP cellArray = vtkCellArraySP::New();
@@ -72,8 +72,8 @@ void ExportVTK( const std::string fileName,
                     break;
                 default:
                     auto log = spdlog::get( "event" );
-                    log->error( "[ERROR][IO::ExportVTK] Invalid dimension" );
-                    exit( EXIT_FAILURE );
+                    ErrorMessages::Error( "Please implement output for results of size " + std::to_string( results[i].size() ) + "!",
+                                          CURRENT_FUNCTION );
             }
             grid->GetCellData()->AddArray( cellData );
         }
@@ -106,7 +106,8 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
     std::vector<std::vector<unsigned>> cells;
     std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries;
 
-    if( !std::filesystem::exists( settings->GetMeshFile() ) ) exit( EXIT_FAILURE );
+    if( !std::filesystem::exists( settings->GetMeshFile() ) )
+        ErrorMessages::Error( "Cannot find mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
     std::ifstream ifs( settings->GetMeshFile(), std::ios::in );
     std::string line;
 
@@ -181,7 +182,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
                             }
                         }
                         else {
-                            exit( EXIT_FAILURE );
+                            ErrorMessages::Error( "Invalid mesh file detected! Make sure boundaries are provided.'", CURRENT_FUNCTION );
                         }
                     }
                     std::sort( bnodes.begin(), bnodes.end() );
@@ -215,8 +216,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
                             break;
                         }
                         default: {
-                            log->error( "[LoadSU2MeshFromFile] Unsupported mesh type!" );
-                            exit( EXIT_FAILURE );
+                            ErrorMessages::Error( "Unsupported mesh type!'", CURRENT_FUNCTION );
                         }
                     }
                 }
@@ -225,8 +225,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
         }
         bool mixedElementMesh = !std::equal( numNodesPerCell.begin() + 1, numNodesPerCell.end(), numNodesPerCell.begin() );
         if( mixedElementMesh ) {
-            log->error( "[LoadSU2MeshFromFile] Mixed element meshes are currently not supported!" );
-            exit( EXIT_FAILURE );
+            ErrorMessages::Error( "Mixed element meshes are currently not supported!'", CURRENT_FUNCTION );
         }
         ifs.clear();
         ifs.seekg( 0, std::ios::beg );
@@ -254,8 +253,7 @@ Mesh* LoadSU2MeshFromFile( const Config* settings ) {
         }
     }
     else {
-        log->error( "[LoadSU2MeshFromFile] File not found" );
-        exit( EXIT_FAILURE );
+        ErrorMessages::Error( "Cannot open mesh file '" + settings->GetMeshFile() + "!", CURRENT_FUNCTION );
     }
     ifs.close();
     return new Mesh( nodes, cells, boundaries );
@@ -281,8 +279,7 @@ std::string ParseArguments( int argc, char* argv[] ) {
             inputFile = std::string( argv[i] );
             std::ifstream f( inputFile );
             if( !f.is_open() ) {
-                std::cerr << "[ERROR] Unable to open inputfile '" << inputFile << "'!" << std::endl;
-                exit( EXIT_FAILURE );
+                ErrorMessages::OptionNotSetError( "Unable to open inputfile '" + inputFile + "' !", CURRENT_FUNCTION );
             }
         }
     }

code/src/main.cpp

@@ -6,7 +6,6 @@
 #include "settings/config.h"
 
 int main( int argc, char** argv ) {
-
     MPI_Init( &argc, &argv );
 
     std::string filename = "default.cfg";

code/src/mesh.cpp

@@ -82,10 +82,9 @@ void Mesh::ComputeConnectivity() {
 
     if( std::any_of( neighborsFlat.begin(), neighborsFlat.end(), []( int i ) { return i == -1; } ) ) {
         for( unsigned idx = 0; idx < neighborsFlat.size(); ++idx ) {
-            if( neighborsFlat[idx] == -1 ) _log->error( "[mesh] Detected unassigned faces at index {0}", idx );
-            break;
+            if( neighborsFlat[idx] == -1 )
+                ErrorMessages::Error( "Detected unassigned faces at index " + std::to_string( idx ) + " !", CURRENT_FUNCTION );
         }
-        exit( EXIT_FAILURE );
     }
 
     _cellNeighbors.resize( _numCells );
@@ -153,7 +152,8 @@ void Mesh::ComputeCellAreas() {
                 break;
             }
             default: {
-                exit( EXIT_FAILURE );
+                ErrorMessages::Error( "Area computation for cells with " + std::to_string( _numNodesPerCell ) + " nodes is not implemented yet!",
+                                      CURRENT_FUNCTION );
             }
         }
     }
@@ -320,11 +320,6 @@ void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiD
     }
 }
 
-unsigned Mesh::GetDim() const { return _dim; }
-unsigned Mesh::GetNumCells() const { return _numCells; }
-unsigned Mesh::GetNumNodes() const { return _numNodes; }
-unsigned Mesh::GetNumNodesPerCell() const { return _numNodesPerCell; }
-
 const std::vector<Vector>& Mesh::GetNodes() const { return _nodes; }
 const std::vector<Vector>& Mesh::GetCellMidPoints() const { return _cellMidPoints; }
 const std::vector<std::vector<unsigned>>& Mesh::GetCells() const { return _cells; }

code/src/quadratures/qgausslegendretensorized.cpp

@@ -7,13 +7,60 @@ QGaussLegendreTensorized::QGaussLegendreTensorized( unsigned order ) : Quadratur
     SetConnectivity();
 }
 
-void QGaussLegendreTensorized::SetPointsAndWeights() {    // TODO
-    // Compute Points
-    // Nq random points on the sphere.
-    _points = VectorVector( GetNq() );
+void QGaussLegendreTensorized::SetPointsAndWeights() {
+    Vector nodes( _order ), weights( _order );
 
-    // Compute Weights
-    _weights = Vector( GetNq(), 4.0 * M_PI / GetNq() );
+    // construct companion matrix
+    Matrix CM( _order, _order, 0.0 );
+
+    for( unsigned i = 0; i < _order - 1; ++i ) {
+        CM( i + 1, i ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
+        CM( i, i + 1 ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
+    }
+
+    auto evSys = ComputeEigenValTriDiagMatrix( CM );
+
+    for( unsigned i = 0; i < _order; ++i ) {
+        if( std::fabs( evSys.first[i] ) < 1e-15 )
+            nodes[i] = 0;
+        else
+            nodes[i] = evSys.first[i];
+        weights[i] = 2 * std::pow( evSys.second( 0, i ), 2 );
+    }
+    for( unsigned i = 0; i < _order; ++i ) {
+        nodes[i] = ( nodes[i] + 1.0 ) * 0.5;
+    }
+
+    std::vector<unsigned> p( nodes.size() );
+    std::iota( p.begin(), p.end(), 0 );
+    std::sort( p.begin(), p.end(), [&]( unsigned i, unsigned j ) { return nodes[i] < nodes[j]; } );
+    Vector sorted_nodes( static_cast<unsigned>( p.size() ) ), sorted_weights( static_cast<unsigned>( p.size() ) );
+    std::transform( p.begin(), p.end(), sorted_nodes.begin(), [&]( unsigned i ) { return nodes[i]; } );
+    std::transform( p.begin(), p.end(), sorted_weights.begin(), [&]( unsigned i ) { return weights[i]; } );
+    nodes   = sorted_nodes;
+    weights = sorted_weights;
+
+    Vector phi( 2 * _order );
+    for( unsigned i = 0; i < 2 * _order; ++i ) {
+        phi[i] = ( i + 0.5 ) * M_PI / _order;
+    }
+
+    unsigned range = std::floor( _order / 2.0 );
+
+    _points.resize( _nq );
+    for( auto& p : _points ) {
+        p.resize( 3 );
+    }
+    _weights.resize( _nq );
+
+    for( unsigned j = 0; j < range; ++j ) {
+        for( unsigned i = 0; i < 2 * _order; ++i ) {
+            _points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes[j] * nodes[j] ) * std::cos( phi[i] );
+            _points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes[j] * nodes[j] ) * std::sin( phi[i] );
+            _points[j * ( 2 * _order ) + i][2] = nodes[j];
+            _weights[j * ( 2 * _order ) + i]   = 2.0 * M_PI / _order * weights[j];
+        }
+    }
 }
 
 void QGaussLegendreTensorized::SetConnectivity() {    // TODO
@@ -21,3 +68,74 @@ void QGaussLegendreTensorized::SetConnectivity() {    // TODO
     VectorVectorU connectivity;
     _connectivity = connectivity;
 }
+
+std::pair<Vector, Matrix> QGaussLegendreTensorized::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
+    unsigned n = mat.rows();
+
+    Vector d( n, 0.0 ), e( n, 0.0 );
+    Matrix z( n, n, 0.0 );
+    for( unsigned i = 0; i < n; ++i ) {
+        d[i]          = mat( i, i );
+        z( i, i )     = 1.0;
+        i == 0 ? e[i] = 0.0 : e[i] = mat( i, i - 1 );
+    }
+
+    int m, l, iter, i, k;
+    m = l = iter = i = k = 0;
+    double s, r, p, g, f, dd, c, b;
+    s = r = p = g = f = dd = c = b = 0.0;
+    const double eps               = std::numeric_limits<double>::epsilon();
+    for( i = 1; i < static_cast<int>( n ); i++ ) e[i - 1] = e[i];
+    e[n - 1] = 0.0;
+    for( l = 0; l < static_cast<int>( n ); l++ ) {
+        iter = 0;
+        do {
+            for( m = l; m < static_cast<int>( n ) - 1; m++ ) {
+                dd = std::fabs( d[m] ) + std::fabs( d[m + 1] );