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

code/include/common/globalconstants.h

@@ -103,6 +103,7 @@ enum SOLVER_NAME {
     CSD_SN_FOKKERPLANCK_SOLVER,
     CSD_SN_FOKKERPLANCK_TRAFO_SOLVER,
     CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D,
+    CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D,
     PN_SOLVER,
     MN_SOLVER
 };
@@ -113,6 +114,7 @@ inline std::map<std::string, SOLVER_NAME> Solver_Map{ { "SN_SOLVER", SN_SOLVER }
                                                       { "CSD_SN_FOKKERPLANCK_SOLVER", CSD_SN_FOKKERPLANCK_SOLVER },
                                                       { "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER },
                                                       { "CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D", CSD_SN_FOKKERPLANCK_TRAFO_SOLVER_2D },
+                                                      { "CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D", CSD_SN_FOKKERPLANCK_TRAFO_SH_SOLVER_2D },
                                                       { "PN_SOLVER", PN_SOLVER },
                                                       { "MN_SOLVER", MN_SOLVER } };
 

code/include/common/io.h

 /*!
- * @file io.h
- * @brief Creates required input from given config file
+ * @file: io.h
+ * @brief: Creates required input from given config file
+ * @author: J. Wolters
  */
+
 #ifndef IO_H
 #define IO_H
 
@@ -14,7 +16,7 @@ class Config;
 class Mesh;
 
 /*!
- * @brief 
+ * @brief
  * @param[in] fileName -name of the file that results shall be exported to
  * @param[in] results - results of computation
  * @param[in] fieldNames - specify different output values
@@ -30,7 +32,6 @@ void ExportVTK( const std::string fileName,
  */
 Mesh* LoadSU2MeshFromFile( const Config* settings );
 
-
 /*!
  * @brief Parses arguments given when calling program from command line
  * @param[in] argc - number arguments

code/include/common/mesh.h

@@ -21,7 +21,8 @@ class Mesh
     const unsigned _numNodes;        /*! @brief: number of nodes in the mesh (for node centered view)*/
     const unsigned _numNodesPerCell; /*! @brief: number of nodes per cell */
     const unsigned _numBoundaries;   /*! @brief: number of boundary cells in the mesh */
-    const unsigned _ghostCellID;     // used for what?    /*! @brief: equal to _numCells and therefore has the ID of the last cell + 1 */
+    const unsigned _ghostCellID; /*! @brief: Id of the ghost cell. (we use only one ghost cell). equal to _numCells and therefore has the ID of the
+                                    last cell + 1 */
 
     std::vector<std::pair<double, double>> _bounds;    // ???
 
@@ -44,8 +45,8 @@ class Mesh
 
     void ComputeCellAreas();     /*! @brief: Computes only the areas of the mesh cells. Write to _cellAreas. */
     void ComputeCellMidpoints(); /*! @brief: Compute only the midpoints of the cells. Write to _cellMidPoints*/
-    void ComputeConnectivity();  /*! @brief: Computes ??? */
-    void ComputePartitioning();  /*! @brief: Computes ??? */
+    void ComputeConnectivity();  /*! @brief: Computes _cellNeighbors and _nodeNeighbors, i.e. neighborship relation in mesh*/
+    void ComputePartitioning();  /*! @brief: Computes local partitioning for openMP */
 
     /*! @brief: Computes outward facing normal of two neighboring nodes nodeA and nodeB with common cellCellcenter.
      *          Normals are scaled with their respective edge length
@@ -54,7 +55,7 @@ class Mesh
      *  @param: cellCenter: Center of the cell that has nodeA and nodeB as nodes.
      *  @return: outward facing normal */
     Vector ComputeOutwardFacingNormal( const Vector& nodeA, const Vector& nodeB, const Vector& cellCenter );
-    void ComputeBounds(); /*! @brief: Computes ??? */
+    void ComputeBounds(); /*! @brief: Computes the spatial bounds of a 2D domain. */
 
   public:
     Mesh() = delete;    //  no default constructor
@@ -111,7 +112,7 @@ class Mesh
      *  @return dimension: scalar */
     double GetDistanceToOrigin( unsigned idx_cell ) const;
 
-    /*! @brief ComputeSlopes calculates the slope in every cell into x and y direction
+    /*! @brief ComputeSlopes calculates the slope in every cell into x and y direction using the divergence theorem.
      *  @param nq is number of quadrature points
      *  @param psiDerX is slope in x direction (gets computed. Slope is stored here)
      *  @param psiDerY is slope in y direction (gets computed. Slope is stored here)
@@ -119,7 +120,17 @@ class Mesh
     // Not used
     void ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
 
+    /*! @brief:Structured mesh slope reconstruction with flux limiters.
+     *  @param nq is number of quadrature points
+     *  @param psiDerX is slope in x direction (gets computed. Slope is stored here)
+     *  @param psiDerY is slope in y direction (gets computed. Slope is stored here)
+     *  @param psi is solution for which slope is computed */
     void ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
+    /*! @brief: Use gauss theorem and limiters. For unstructured mesh *
+     *  @param nq is number of quadrature points
+     *  @param psiDerX is slope in x direction (gets computed. Slope is stored here)
+     *  @param psiDerY is slope in y direction (gets computed. Slope is stored here)
+     *  @param psi is solution for which slope is computed */
     void ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
 };
 

code/include/kernels/isotropic.h

+/*!
+ * @file isotropic1D.h
+ * @brief Class for computing an isotropic scattering kernel of kinetic equations
+ * @author ?
+ */
+
 #ifndef ISOTROPIC_H
 #define ISOTROPIC_H
 

code/include/kernels/isotropic1D.h

+/*!
+ * @file isotropic1D.h
+ * @brief Class for computing an isotropic scattering kernel of 1D kinetic equations
+ * @author ?
+ */
+
 #ifndef ISOTROPIC1D_H
 #define ISOTROPIC1D_H
 

code/include/kernels/scatteringkernelbase.h

+/*!
+ * @file scatteringkernelbase.h
+ * @brief Base class for computing the scattering kernel of kinetic equations
+ * @author ?
+ */
+
 #ifndef SCATTERINGKERNELBASE_CPP
 #define SCATTERINGKERNELBASE_CPP
 
@@ -13,14 +19,19 @@ class ScatteringKernel
     ScatteringKernel() = delete;
 
   protected:
-    QuadratureBase* _quad;
+    QuadratureBase* _quad; /*! @brief: Pointer to the quadrature used to compute the scattering integral */
 
   public:
+    /*! @brief: Copies the pointer of quad. Does not create an own quad */
     ScatteringKernel( QuadratureBase* quad );
+
     virtual ~ScatteringKernel();
 
+    /*! @brief: Computes the scattering kernel and for the whole SN system and stores it in a Matrix
+        @return: Matrix with discretized scattering kernel */
     virtual Matrix GetScatteringKernel() = 0;
 
+    /*! @brief: Creates an object of the child class of ScatteringKernelBase corresponding to the enum KERNEL_NAME */
     static ScatteringKernel* CreateScatteringKernel( KERNEL_NAME name, QuadratureBase* quad );
 };
 

code/include/optimizers/newtonoptimizer.h

+/*!
+ * @file newtonoptimizer.h
+ * @brief class for solving the minimal entropy optimization problem using a newton optimizer with line search.
+ * @author S. Schotthöfer
+ */
+
 #ifndef NEWTONOPTIMIZER_H
 #define NEWTONOPTIMIZER_H
 
@@ -10,7 +16,7 @@ class NewtonOptimizer : public OptimizerBase
   public:
     NewtonOptimizer( Config* settings );
 
-    inline ~NewtonOptimizer() {}
+    ~NewtonOptimizer();
 
     void Solve( Vector& lambda, Vector& u, VectorVector& moments ) override;
 

code/include/optimizers/optimizerbase.h

+/*!
+ * @file optimizerbase.h
+ * @brief Base class for solving the minimal entropy optimization problem
+ * @author S. Schotthöfer
+ */
+
 #ifndef OPTIMIZERBASE_H
 #define OPTIMIZERBASE_H
 
 #include "common/typedef.h"
-#include "entropies/entropybase.h"
 
 // Foward declaration
 class Config;
+class EntropyBase;
 
 class OptimizerBase
 {
   public:
     OptimizerBase( Config* settings );
 
-    virtual inline ~OptimizerBase() { delete _entropy; }
+    ~OptimizerBase();
 
+    /*! @brief: Optimizer creator: Depending on the chosen option, this function creates an object of the chosen child class of OptimizerBase */
     static OptimizerBase* Create( Config* settings );
 
     /*! @brief  : Computes the optimal Lagrange multilpiers for the dual entropy minimization problem
@@ -23,7 +30,7 @@ class OptimizerBase
 
   protected:
     EntropyBase* _entropy; /*! @brief: Class to handle entropy functional evaluations */
-    Config* _settings;
+    Config* _settings;     /*! @biref: Pointer to settings class of the solver */
 };
 
 #endif    // OPTIMIZERBASE_H

code/include/problems/checkerboard.h

@@ -6,13 +6,13 @@
 class Checkerboard_SN : public ProblemBase
 {
   private:
-    Vector _scatteringXS;
-    Vector _totalXS;
+    Vector _scatteringXS;     /*! @brief Vector of scattering crosssections */
+    Vector _totalXS;      /*! @brief Vector of total crosssections */
 
     Checkerboard_SN() = delete;
 
-    bool isAbsorption( const Vector& pos ) const;
-    bool isSource( const Vector& pos ) const;
+    bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
+    bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
 
   public:
     Checkerboard_SN( Config* settings, Mesh* mesh );
@@ -27,13 +27,22 @@ class Checkerboard_SN : public ProblemBase
 class Checkerboard_PN : public ProblemBase
 {
   private:
-    Vector _scatteringXS;
-    Vector _totalXS;
+    Vector _scatteringXS;     /*! @brief Vector of scattering crosssections */
+    Vector _totalXS;      /*! @brief Vector of total crosssections */
 
     Checkerboard_PN() = delete;
 
-    bool isAbsorption( const Vector& pos ) const;
-    bool isSource( const Vector& pos ) const;
+    bool isAbsorption( const Vector& pos ) const; /*! @return True if pos is in absorption region, False otherwise */
+    bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
+
+    /**
+     * @brief Gets the global index for given order l of Legendre polynomials and given
+     *        order k of Legendre functions.
+     *        Note: This is code doubling from PNSolver::GlobalIndex
+     * @param l : order of Legendre polynomial
+     * @param k : order of Legendre function
+     * @returns global index
+     */
     int GlobalIndex( int l, int k ) const;
 
   public:

code/include/problems/problembase.h

@@ -16,8 +16,8 @@ class ProblemBase
     Mesh* _mesh;
     Physics* _physics;
 
-    std::vector<double> _density;
-    std::vector<double> _stoppingPower;
+    std::vector<double> _density; /*! @brief: vector with patient densities */
+    std::vector<double> _stoppingPower; /*! @brief: vector with stopping powers*/
 
     ProblemBase() = delete;
 
@@ -37,6 +37,7 @@ class ProblemBase
      * @param density is vector with patient densities (at different spatial cells)
      */
     virtual VectorVector GetTotalXS( const Vector& energies ) = 0;
+    
     /**
      * @brief GetTotalXSE gives back vector of total cross sections for
      *        energies in vector energy

code/include/reconstructor.h

+/*!
+ * @file: reconstructor.h
+ * @brief: Class to create second order (in space) schemes for the advection solver.
+ *         This class is currently unused. But in the future, the second order capabilities of the code
+ *         may be stored here.
+ * @author: T. Xiao
+ */
+
 #ifndef RECONSTRUCTOR_H
 #define RECONSTRUCTOR_H
 
-#include "common/config.h"
-#include "common/typedef.h"
+#include <string>
+
+class Config;
 
 class Reconstructor
 {
@@ -13,7 +22,7 @@ class Reconstructor
      */
     Reconstructor( Config* settings );
 
-    /** Method 1: structured developing
+    /*! Method 1: structured developing
      * @brief Slope of angular flux psi inside a given cell
      * @param Omega fixed ordinate for flux computation
      * @param psiL left solution state

code/include/solvers/csdsnsolverfp.h

@@ -12,7 +12,7 @@ class CSDSNSolverFP : public SNSolver
     std::vector<double> _dose; /*! @brief: TODO */
 
     // Physics acess
-    Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
+    Vector _energies; /*! @brief: energy levels for CSD, length = _nEnergies */
     Vector _angle;    /*! @brief: angles for SN */
 
     std::vector<Matrix> _sigmaSE; /*!  @brief scattering cross section for all energies*/
@@ -21,17 +21,15 @@ class CSDSNSolverFP : public SNSolver
     Matrix _L;  /*!  @brief Laplace Beltrami Matrix */
     Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
 
-    double _alpha;
-    double _beta;
+    double _alpha; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, eq. (8)*/
+    double _beta; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, eq. (8)*/
 
-    Vector _xi1;
-    Vector _xi2;
-    Matrix _xi;
+    Matrix _xi;  /*!  @brief matrix of transport coefficients */
 
-    bool _RT;
+    bool _RT; /*!  @brief radiotherapy application (on/off), if true use crosssections + stopping powers from database  */
 
-    double _energyMin;
-    double _energyMax;
+    double _energyMin; /*!  @brief minimal energy in energy grid*/
+    double _energyMax; /*!  @brief maximal energy in energy grid*/
 
   public:
     /**

code/include/solvers/csdsolvertrafofp.h

@@ -21,20 +21,20 @@ class CSDSolverTrafoFP : public SNSolver
     Matrix _L;  /*!  @brief Laplace Beltrami Matrix */
     Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
 
-    double _alpha;
-    double _alpha2;
-    double _beta;
+    double _alpha; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
+    double _alpha2; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
+    double _beta; /*!  @brief  Coefficient of GFP operators (see Olbrant 2010, Appendix B)*/
 
+    Matrix _xi;  /*!  @brief matrix of transport coefficients */
     Vector _xi1;
     Vector _xi2;
-    Matrix _xi;
+    
+    unsigned _FPMethod; /*!  @brief Encodes different ways of computing coefficients alpha, alpha2 & beta, _FPMethod == 1, 2 ,3 stand for methods with increasing accuracy (see Olbrant 2010, Appendix B)*/
 
-    unsigned _FPMethod;
+    bool _RT; /*!  @brief radiotherapy application (on/off), if true use crosssections + stopping powers from database  */
 
-    bool _RT;
-
-    double _energyMin;
-    double _energyMax;
+    double _energyMin; /*!  @brief minimal energy in energy grid*/
+    double _energyMax; /*!  @brief maximal energy in energy grid*/
 
     void GenerateEnergyGrid( bool refinement );
 

code/include/solvers/csdsolvertrafofpsh2d.h

+#ifndef CSDSOLVERTRAFOFPSH2D_H
+#define CSDSOLVERTRAFOFPSH2D_H
+
+// externals
+#include "spdlog/spdlog.h"
+#include <mpi.h>
+
+#include "common/config.h"
+#include "common/io.h"
+#include "fluxes/numericalflux.h"
+#include "icru.h"
+#include "kernels/scatteringkernelbase.h"
+#include "problems/problembase.h"
+#include "quadratures/quadraturebase.h"
+#include "solvers/snsolver.h"
+#include "sphericalharmonics.h"
+
+class Physics;
+
+class CSDSolverTrafoFPSH2D : public SNSolver
+{
+  private:
+    std::vector<double> _dose; /*! @brief: TODO */
+
+    // Physics acess
+    Vector _energies; /*! @brief: energy levels for CSD, lenght = _nEnergies */
+    Vector _angle;    /*! @brief: angles for SN */
+
+    std::vector<Matrix> _sigmaSE; /*!  @brief scattering cross section for all energies*/
+    Vector _sigmaTE;              /*!  @brief total cross section for all energies*/
+
+    Matrix _L;  /*!  @brief Laplace Beltrami Matrix */
+    Matrix _IL; /*!  @brief Laplace Beltrami Matrix */
+
+    Matrix _O;
+    Matrix _S;
+    Matrix _M;
+
+    VectorVector _quadPoints;
+    VectorVector _quadPointsSphere;
+    Vector _weights;
+    Vector _mu;
+    Vector _phi;
+    Vector _wp;
+    Vector _wa;
+
+    double _alpha;
+    double _alpha2;
+    double _beta;
+
+    Vector _xi1;
+    Vector _xi2;
+    Matrix _xi;
+
+    unsigned _FPMethod;
+
+    bool _RT;
+
+    double _energyMin;
+    double _energyMax;
+
+    void GenerateEnergyGrid( bool refinement );
+
+  public:
+    /**
+     * @brief CSDSolverTrafoFP2D constructor
+     * @param settings stores all needed information
+     */
+    CSDSolverTrafoFPSH2D( Config* settings );
+    /**
+     * @brief Solve functions runs main time loop
+     */
+    virtual void Solve();
+    /**
+     * @brief Output solution to VTK file
+     */
+    virtual void Save() const;
+    virtual void Save( int currEnergy ) const;
+};
+
+#endif    // CSDSOLVERTRAFOFPSH2D_H

code/src/optimizers/newtonoptimizer.cpp

 
 #include "optimizers/newtonoptimizer.h"
 #include "common/config.h"
+#include "entropies/entropybase.h"
 #include "quadratures/quadraturebase.h"
 #include "toolboxes/errormessages.h"
 
@@ -15,6 +16,8 @@ NewtonOptimizer::NewtonOptimizer( Config* settings ) : OptimizerBase( settings )
     _epsilon          = settings->GetNewtonOptimizerEpsilon();
 }
 
+NewtonOptimizer::~NewtonOptimizer() { delete _quadrature; }
+
 double NewtonOptimizer::ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments ) {
     double result = 0.0;
 

code/src/optimizers/optimizerbase.cpp

 #include "optimizers/optimizerbase.h"
 #include "common/config.h"
+#include "entropies/entropybase.h"
 #include "optimizers/newtonoptimizer.h"
 
 OptimizerBase::OptimizerBase( Config* settings ) {
@@ -7,6 +8,8 @@ OptimizerBase::OptimizerBase( Config* settings ) {
     _settings = settings;
 }
 
+OptimizerBase::~OptimizerBase() { delete _entropy; }
+
 OptimizerBase* OptimizerBase::Create( Config* settings ) {
     switch( settings->GetOptimizerName() ) {
         case NEWTON: return new NewtonOptimizer( settings );

code/src/problems/checkerboard.cpp

@@ -2,10 +2,16 @@
 #include "common/config.h"
 #include "common/mesh.h"
 
+// ---- Checkerboard Sn ----
+//Constructor for Ckeckerboard case with Sn
 Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     _physics      = nullptr;
+
+    // Initialise crosssections to 1
     _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
     _totalXS      = Vector( _mesh->GetNumCells(), 1.0 );
+
+    // For absorption cells: set scattering XS to 0 and absorption to 10
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         if( isAbsorption( cellMids[j] ) ) {
@@ -36,6 +42,7 @@ VectorVector Checkerboard_SN::SetupIC() {
 }
 
 bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
+    // Check whether pos is inside absorbing squares
     std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
     std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
     for( unsigned k = 0; k < lbounds.size(); ++k ) {
@@ -50,6 +57,7 @@ bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
 }
 
 bool Checkerboard_SN::isSource( const Vector& pos ) const {
+    // Check whether pos is part of source region
     if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
         return true;
     else
@@ -58,10 +66,16 @@ bool Checkerboard_SN::isSource( const Vector& pos ) const {
 
 ////////////////////////////////////////////////////////////////////////////////////////////////////
 
+// ---- Checkerboard Pn ----
+//Constructor for checkerboard case with Pn
 Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     _physics      = nullptr;
+
+    //Initialise crosssections = 1 (scattering)
     _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
     _totalXS      = Vector( _mesh->GetNumCells(), 1.0 );
+
+    // for absorption regions change crosssections to all absorption
     auto cellMids = _mesh->GetCellMidPoints();
     for( unsigned j = 0; j < cellMids.size(); ++j ) {
         if( isAbsorption( cellMids[j] ) ) {
@@ -93,6 +107,7 @@ VectorVector Checkerboard_PN::SetupIC() {
 }
 
 bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
+    // Check whether pos is in absorption region
     std::vector<double> lbounds{ 1, 2, 3, 4, 5 };
     std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
     for( unsigned k = 0; k < lbounds.size(); ++k ) {
@@ -107,6 +122,7 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
 }
 
 bool Checkerboard_PN::isSource( const Vector& pos ) const {
+    // Check whether pos is in source region
     if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
         return true;
     else

code/src/problems/electronrt.cpp

@@ -2,6 +2,7 @@
 #include "common/config.h"
 #include "common/mesh.h"
 
+//Constructor: Legacy code, physics class is no longer used
 ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
     _physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" );
 }
@@ -9,37 +10,46 @@ ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings,
 ElectronRT::~ElectronRT() { delete _physics; }
 
 VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
 VectorVector ElectronRT::GetTotalXS( const Vector& energies ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
 }
 
 std::vector<Matrix> ElectronRT::GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return _physics->GetScatteringXS( energies, angles );
 }
 
 Vector ElectronRT::GetTotalXSE( const Vector& energies ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return _physics->GetTotalXSE( energies );
 }
 
 std::vector<VectorVector> ElectronRT::GetExternalSource( const Vector& energies ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
 }
 
 VectorVector ElectronRT::SetupIC() {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
     return VectorVector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
 }
 
 void ElectronRT::LoadXSH20( std::string fileSigmaS, std::string fileSigmaT ) {
-    // @TODO
+    // @TODO 
+    // Specified in subclasses 
 }
 
+
+// Default densities = 1, for patient files this is overwritten with "real" densities
 std::vector<double> ElectronRT::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }

code/src/problems/musclebonelung.cpp

@@ -22,11 +22,12 @@ VectorVector MuscleBoneLung::SetupIC() {
 }
 
 std::vector<double>  MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) { 
-    std::cout<<"Length of mesh "<< cellMidPoints.size() << std::endl;
     std::vector<double> densities ( 167, 1.04 ); //muscle layer
-    std::vector<double> bone( 167, 1.85  );