Commit f8118df4 authored by pia.stammer's avatar pia.stammer
Browse files

Added comments to problem files

parent a7458217
Pipeline #119960 canceled with stages
...@@ -6,13 +6,13 @@ ...@@ -6,13 +6,13 @@
class Checkerboard_SN : public ProblemBase class Checkerboard_SN : public ProblemBase
{ {
private: private:
Vector _scatteringXS; Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
Vector _totalXS; Vector _totalXS; /*! @brief Vector of total crosssections */
Checkerboard_SN() = delete; Checkerboard_SN() = delete;
bool isAbsorption( 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; bool isSource( const Vector& pos ) const; /*! @return True if pos is in source region, False otherwise */
public: public:
Checkerboard_SN( Config* settings, Mesh* mesh ); Checkerboard_SN( Config* settings, Mesh* mesh );
...@@ -27,13 +27,22 @@ class Checkerboard_SN : public ProblemBase ...@@ -27,13 +27,22 @@ class Checkerboard_SN : public ProblemBase
class Checkerboard_PN : public ProblemBase class Checkerboard_PN : public ProblemBase
{ {
private: private:
Vector _scatteringXS; Vector _scatteringXS; /*! @brief Vector of scattering crosssections */
Vector _totalXS; Vector _totalXS; /*! @brief Vector of total crosssections */
Checkerboard_PN() = delete; Checkerboard_PN() = delete;
bool isAbsorption( 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; 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; int GlobalIndex( int l, int k ) const;
public: public:
......
...@@ -16,8 +16,8 @@ class ProblemBase ...@@ -16,8 +16,8 @@ class ProblemBase
Mesh* _mesh; Mesh* _mesh;
Physics* _physics; Physics* _physics;
std::vector<double> _density; std::vector<double> _density; /*! @brief: vector with patient densities */
std::vector<double> _stoppingPower; std::vector<double> _stoppingPower; /*! @brief: vector with stopping powers*/
ProblemBase() = delete; ProblemBase() = delete;
...@@ -37,6 +37,7 @@ class ProblemBase ...@@ -37,6 +37,7 @@ class ProblemBase
* @param density is vector with patient densities (at different spatial cells) * @param density is vector with patient densities (at different spatial cells)
*/ */
virtual VectorVector GetTotalXS( const Vector& energies ) = 0; virtual VectorVector GetTotalXS( const Vector& energies ) = 0;
/** /**
* @brief GetTotalXSE gives back vector of total cross sections for * @brief GetTotalXSE gives back vector of total cross sections for
* energies in vector energy * energies in vector energy
......
...@@ -2,10 +2,16 @@ ...@@ -2,10 +2,16 @@
#include "common/config.h" #include "common/config.h"
#include "common/mesh.h" #include "common/mesh.h"
// ---- Checkerboard Sn ----
//Constructor for Ckeckerboard case with Sn
Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { Checkerboard_SN::Checkerboard_SN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr; _physics = nullptr;
// Initialise crosssections to 1
_scatteringXS = Vector( _mesh->GetNumCells(), 1.0 ); _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
_totalXS = 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(); auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) { for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( isAbsorption( cellMids[j] ) ) { if( isAbsorption( cellMids[j] ) ) {
...@@ -36,6 +42,7 @@ VectorVector Checkerboard_SN::SetupIC() { ...@@ -36,6 +42,7 @@ VectorVector Checkerboard_SN::SetupIC() {
} }
bool Checkerboard_SN::isAbsorption( const Vector& pos ) const { 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> lbounds{ 1, 2, 3, 4, 5 };
std::vector<double> ubounds{ 2, 3, 4, 5, 6 }; std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
for( unsigned k = 0; k < lbounds.size(); ++k ) { for( unsigned k = 0; k < lbounds.size(); ++k ) {
...@@ -50,6 +57,7 @@ bool Checkerboard_SN::isAbsorption( const Vector& pos ) const { ...@@ -50,6 +57,7 @@ bool Checkerboard_SN::isAbsorption( const Vector& pos ) const {
} }
bool Checkerboard_SN::isSource( 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 ) if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
return true; return true;
else else
...@@ -58,10 +66,16 @@ bool Checkerboard_SN::isSource( const Vector& pos ) const { ...@@ -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 ) { Checkerboard_PN::Checkerboard_PN( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr; _physics = nullptr;
//Initialise crosssections = 1 (scattering)
_scatteringXS = Vector( _mesh->GetNumCells(), 1.0 ); _scatteringXS = Vector( _mesh->GetNumCells(), 1.0 );
_totalXS = Vector( _mesh->GetNumCells(), 1.0 ); _totalXS = Vector( _mesh->GetNumCells(), 1.0 );
// for absorption regions change crosssections to all absorption
auto cellMids = _mesh->GetCellMidPoints(); auto cellMids = _mesh->GetCellMidPoints();
for( unsigned j = 0; j < cellMids.size(); ++j ) { for( unsigned j = 0; j < cellMids.size(); ++j ) {
if( isAbsorption( cellMids[j] ) ) { if( isAbsorption( cellMids[j] ) ) {
...@@ -93,6 +107,7 @@ VectorVector Checkerboard_PN::SetupIC() { ...@@ -93,6 +107,7 @@ VectorVector Checkerboard_PN::SetupIC() {
} }
bool Checkerboard_PN::isAbsorption( const Vector& pos ) const { 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> lbounds{ 1, 2, 3, 4, 5 };
std::vector<double> ubounds{ 2, 3, 4, 5, 6 }; std::vector<double> ubounds{ 2, 3, 4, 5, 6 };
for( unsigned k = 0; k < lbounds.size(); ++k ) { for( unsigned k = 0; k < lbounds.size(); ++k ) {
...@@ -107,6 +122,7 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const { ...@@ -107,6 +122,7 @@ bool Checkerboard_PN::isAbsorption( const Vector& pos ) const {
} }
bool Checkerboard_PN::isSource( 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 ) if( pos[0] >= 3 && pos[0] <= 4 && pos[1] >= 3 && pos[1] <= 4 )
return true; return true;
else else
......
...@@ -2,6 +2,7 @@ ...@@ -2,6 +2,7 @@
#include "common/config.h" #include "common/config.h"
#include "common/mesh.h" #include "common/mesh.h"
//Constructor: Legacy code, physics class is no longer used
ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) { ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" ); _physics = new Physics( settings->GetHydrogenFile(), settings->GetOxygenFile(), "../input/stopping_power.txt" );
} }
...@@ -9,37 +10,46 @@ ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, ...@@ -9,37 +10,46 @@ ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings,
ElectronRT::~ElectronRT() { delete _physics; } ElectronRT::~ElectronRT() { delete _physics; }
VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) { VectorVector ElectronRT::GetScatteringXS( const Vector& energies ) {
// @TODO // @TODO
// Specified in subclasses
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
} }
VectorVector ElectronRT::GetTotalXS( const Vector& energies ) { VectorVector ElectronRT::GetTotalXS( const Vector& energies ) {
// @TODO // @TODO
// Specified in subclasses
return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) ); return VectorVector( energies.size(), Vector( _mesh->GetNumCells(), 0.0 ) );
} }
std::vector<Matrix> ElectronRT::GetScatteringXSE( const Vector& energies, const Matrix& angles ) { std::vector<Matrix> ElectronRT::GetScatteringXSE( const Vector& energies, const Matrix& angles ) {
// @TODO // @TODO
// Specified in subclasses
return _physics->GetScatteringXS( energies, angles ); return _physics->GetScatteringXS( energies, angles );
} }
Vector ElectronRT::GetTotalXSE( const Vector& energies ) { Vector ElectronRT::GetTotalXSE( const Vector& energies ) {
// @TODO // @TODO
// Specified in subclasses
return _physics->GetTotalXSE( energies ); return _physics->GetTotalXSE( energies );
} }
std::vector<VectorVector> ElectronRT::GetExternalSource( const Vector& 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 ) ) ); return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
} }
VectorVector ElectronRT::SetupIC() { VectorVector ElectronRT::SetupIC() {
// @TODO // @TODO
// Specified in subclasses
return VectorVector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) ); return VectorVector( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
} }
void ElectronRT::LoadXSH20( std::string fileSigmaS, std::string fileSigmaT ) { 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 ); } std::vector<double> ElectronRT::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
...@@ -22,11 +22,12 @@ VectorVector MuscleBoneLung::SetupIC() { ...@@ -22,11 +22,12 @@ VectorVector MuscleBoneLung::SetupIC() {
} }
std::vector<double> MuscleBoneLung::GetDensity( const VectorVector& cellMidPoints ) { 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> densities ( 167, 1.04 ); //muscle layer
std::vector<double> bone( 167, 1.85 ); std::vector<double> bone( 167, 1.85 ); // bone layer
std::vector<double> lung ( 665, 0.3 ); //maybe this is just the lung tissue and after that there should be air? std::vector<double> lung ( 665, 0.3 ); // lung layer (maybe this is just the lung tissue and after that there should be air?)
//Concatenate layers
densities.insert(densities.end(),bone.begin(),bone.end()); densities.insert(densities.end(),bone.begin(),bone.end());
densities.insert(densities.end(),lung.begin(),lung.end()); densities.insert(densities.end(),lung.begin(),lung.end());
std::cout<<"Length of densities "<< densities.size() << std::endl;
return densities; } return densities; }
...@@ -18,6 +18,8 @@ ProblemBase::~ProblemBase() {} ...@@ -18,6 +18,8 @@ ProblemBase::~ProblemBase() {}
ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) { ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
auto name = settings->GetProblemName(); auto name = settings->GetProblemName();
// Choose problem type
switch( name ) { switch( name ) {
case PROBLEM_LineSource: { case PROBLEM_LineSource: {
if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER ) if( settings->GetSolverName() == PN_SOLVER || settings->GetSolverName() == MN_SOLVER )
...@@ -42,10 +44,13 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) { ...@@ -42,10 +44,13 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
} }
} }
// Default densities = 1
std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); } std::vector<double> ProblemBase::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
// Legacy code: Scattering crossection loaded from database ENDF with physics class -> later overwritten with ICRU data
VectorVector ProblemBase::GetScatteringXSE( const Vector& energies, const Vector& angles ) { return _physics->GetScatteringXS( energies, angles ); } VectorVector ProblemBase::GetScatteringXSE( const Vector& energies, const Vector& angles ) { return _physics->GetScatteringXS( energies, angles ); }
// Stopping powers from phyics class or default = -1
Vector ProblemBase::GetStoppingPower( const Vector& energies ) { Vector ProblemBase::GetStoppingPower( const Vector& energies ) {
if( _physics ) { if( _physics ) {
return _physics->GetStoppingPower( energies ); return _physics->GetStoppingPower( energies );
......
...@@ -15,4 +15,5 @@ VectorVector WaterPhantom::SetupIC() { ...@@ -15,4 +15,5 @@ VectorVector WaterPhantom::SetupIC() {
return psi; return psi;
} }
// Density of water = 1 everywhere
std::vector<double> WaterPhantom::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); } std::vector<double> WaterPhantom::GetDensity( const VectorVector& cellMidPoints ) { return std::vector<double>( cellMidPoints.size(), 1.0 ); }
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment