tensorization removed and problem class used

code/include/problems/electronrt.h

 #ifndef ELECTRONRT_H
 #define ELECTRONRT_H
 
+#include "physics.h"
 #include "problembase.h"
 
 class ElectronRT : public ProblemBase
@@ -21,6 +22,8 @@ class ElectronRT : public ProblemBase
 
     virtual VectorVector GetScatteringXS( const std::vector<double>& energies );
     virtual VectorVector GetTotalXS( const std::vector<double>& energies );
+    virtual VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles );
+    virtual Vector GetTotalXSE( const Vector& energies );
     virtual std::vector<VectorVector> GetExternalSource( const std::vector<double>& energies );
     virtual std::vector<double> GetStoppingPower( const std::vector<double>& energies );
     virtual VectorVector SetupIC();

code/include/problems/problembase.h

@@ -37,6 +37,19 @@ class ProblemBase
      * @param density is vector with patient densities (at different spatial cells)
      */
     virtual VectorVector GetTotalXS( const std::vector<double>& energies ) = 0;
+    /**
+     * @brief GetTotalXSE gives back vector of total cross sections for
+     *        energies in vector energy
+     * @param energy is the energy the cross section is queried for
+     */
+    virtual Vector GetTotalXSE( const Vector& energies ) { return Vector( 1 ); }
+
+    /**
+     * @brief GetScatteringXSE gives back vector (each energy) of scattering cross sections for energies
+     *        in vector energy
+     * @param energy is the energy the cross section is queried for
+     */
+    virtual VectorVector GetScatteringXSE( const Vector& energies, const Vector& angles ) { return VectorVector( 1, Vector( 1 ) ); }
 
     /**
      * @brief GetExternalSource gives back vector of vectors of source terms for each

code/include/solvers/csdsnsolver.h

@@ -16,6 +16,9 @@ class CSDSNSolver : public SNSolver
     Vector _density;   /*! @brief: patient density for each grid cell */
     Physics* _physics; /*! @brief: class that hold physics information */
 
+    VectorVector _sigmaSE; /*!  @brief scattering cross section for all energies*/
+    Vector _sigmaTE;       /*!  @brief total cross section for all energies*/
+
   public:
     /**
      * @brief CSDSNSolver constructor

code/include/solvers/solverbase.h

@@ -29,8 +29,8 @@ class Solver
     std::vector<double> _s;           // stopping power, dim(_s) = _nTimeSteps
     std::vector<VectorVector> _Q;     /*!  @brief  external source term */
 
-    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies */
-    VectorVector _sigmaT; /*!  @brief total cross section for all energies */
+    VectorVector _sigmaS; /*!  @brief scattering cross section for all energies and spatial cells*/
+    VectorVector _sigmaT; /*!  @brief total cross section for all energies and spatial cells*/
 
     // quadrature related numbers
     QuadratureBase* _quadrature; /*! @brief quadrature to create members below */

code/src/problems/electronrt.cpp

@@ -4,7 +4,7 @@
 
 ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     // @TODO
-    _physics = nullptr;
+    _physics = new Physics();
 }
 
 ElectronRT::~ElectronRT() {}
@@ -19,6 +19,17 @@ VectorVector ElectronRT::GetTotalXS( const std::vector<double>& energies ) {
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
+VectorVector ElectronRT::GetScatteringXSE( const Vector& energies, const Vector& angles ) {
+    // @TODO
+    _physics->GetScatteringXS( energies, angles );
+    return VectorVector( energies.size(), Vector( angles.size(), 0.0 ) );
+}
+
+Vector ElectronRT::GetTotalXSE( const Vector& energies ) {
+    // @TODO
+    return Vector( energies.size(), 0.0 );
+}
+
 std::vector<VectorVector> ElectronRT::GetExternalSource( const std::vector<double>& energies ) {
     // @TODO
     return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );

code/src/solvers/csdsnsolver.cpp

@@ -3,7 +3,7 @@
 #include "common/io.h"
 #include "fluxes/numericalflux.h"
 #include "kernels/scatteringkernelbase.h"
-#include "physics.h"
+#include "problems/problembase.h"
 
 CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
     _dose = std::vector<double>( _settings->GetNCells(), 0.0 );
@@ -13,7 +13,8 @@ CSDSNSolver::CSDSNSolver( Config* settings ) : SNSolver( settings ) {
     _energies = Vector( _nEnergies, 0.0 );
     // TODO: write meaningfull values for them!
 
-    _sigmaS = _physics->GetScatteringXS( _energies, _angle );
+    _sigmaSE = _problem->GetScatteringXSE( _energies, _angle );
+    _sigmaTE = _problem->GetTotalXSE( _energies );
 
     // Get patient density
     _density = Vector( _nCells, 0.0 );
@@ -79,10 +80,12 @@ void CSDSNSolver::Solve() {
                 }
                 // time update angular flux with numerical flux and total scattering cross section
                 psiNew[idx_cell][idx_ord] = _sol[idx_cell][idx_ord] - ( _dE / _areas[idx_cell] ) * psiNew[idx_cell][idx_ord] -
-                                            _dE * _sigmaT[idx_energy][idx_cell] * _sol[idx_cell][idx_ord];
+                                            _dE * _density[idx_cell] * _sigmaTE[idx_energy] * _sol[idx_cell][idx_ord];
             }
             // compute scattering effects
-            psiNew[idx_cell] += _dE * _sigmaS[idx_energy][idx_cell] * _scatteringKernel * _sol[idx_cell];    // multiply scattering matrix with psi
+            psiNew[idx_cell] +=
+                _dE * _density[idx_cell] * ( blaze::trans( _sigmaSE[idx_energy] ) * _sol[idx_cell] );    // multiply scattering matrix with psi
+            // TODO: Check if _sigmaS^T*psi is correct
 
             // TODO: figure out a more elegant way
             // add external source contribution