Commit 785ce1c6 authored by pia.stammer's avatar pia.stammer
Browse files

Added 2D phantom test case and changed densities according to ICRU database

parent 67ec4b07
Pipeline #117784 failed with stages
in 27 minutes and 39 seconds
......@@ -75,7 +75,8 @@ enum PROBLEM_NAME {
PROBLEM_LineSource_Pseudo_1D,
PROBLEM_LineSource_Pseudo_1D_Physics,
PROBLEM_AirCavity,
PROBLEM_MuscleBoneLung
PROBLEM_MuscleBoneLung,
PROBLEM_Phantom2D
};
inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_LineSource },
......@@ -84,6 +85,7 @@ inline std::map<std::string, PROBLEM_NAME> Problem_Map{ { "LINESOURCE", PROBLEM_
{ "WATERPHANTOM", PROBLEM_WaterPhantom },
{ "AIRCAVITY", PROBLEM_AirCavity },
{ "MUSCLEBONELUNG", PROBLEM_MuscleBoneLung },
{ "PHANTOM2D", PROBLEM_Phantom2D},
{ "LINESOURCE_PSEUDO_1D", PROBLEM_LineSource_Pseudo_1D },
{ "LINESOURCE_PSEUDO_1D_PHYSICS", PROBLEM_LineSource_Pseudo_1D_Physics } };
......
#ifndef PHANTOM2D
#define PHANTOM2D
#include "electronrt.h"
#include <fstream>
#include <numeric>
#include "common/config.h"
#include "common/io.h"
#include "common/mesh.h"
#include "interpolation.h"
class Phantom2D : public ElectronRT
{
private:
Phantom2D() = delete;
public:
Phantom2D( Config* settings, Mesh* mesh );
virtual ~Phantom2D();
virtual std::vector<VectorVector> GetExternalSource( const Vector& energies );
virtual VectorVector SetupIC();
std::vector<double> GetDensity( const VectorVector& cellMidPoints );
};
#endif // PHANTOM2D
\ No newline at end of file
......@@ -23,9 +23,9 @@ 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.05 ); //muscle layer
std::vector<double> bone( 167, 1.92 );
std::vector<double> lung ( 665, 0.26 );
std::vector<double> densities ( 167, 1.04 ); //muscle layer
std::vector<double> bone( 167, 1.85 );
std::vector<double> lung ( 665, 0.3 ); //maybe this is just the lung tissue and after that there should be air?
densities.insert(densities.end(),bone.begin(),bone.end());
densities.insert(densities.end(),lung.begin(),lung.end());
std::cout<<"Length of densities "<< densities.size() << std::endl;
......
#include "problems/phantom2d.h"
#include "common/config.h"
#include "common/mesh.h"
Phantom2D::Phantom2D( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {}
Phantom2D::~Phantom2D() { delete _physics; }
std::vector<VectorVector> Phantom2D::GetExternalSource( const Vector& energies ) {
return std::vector<VectorVector>( energies.size(), std::vector<Vector>( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 0.0 ) ) );
}
VectorVector Phantom2D::SetupIC() {
VectorVector psi( _mesh->GetNumCells(), Vector( _settings->GetNQuadPoints(), 1e-10 ) );
auto cellMids = _mesh->GetCellMidPoints();
double s = 0.1;
for( unsigned j = 0; j < cellMids.size(); ++j ) {
double x = cellMids[j][0];
psi[j][_settings->GetNQuadPoints() - 1] = 1.0 / ( s * sqrt( 2 * M_PI ) ) * std::exp( -x * x / ( 2 * s * s ) );
}
return psi;
}
std::vector<double> Phantom2D::GetDensity( const VectorVector& cellMidPoints ) {
std::string imageFile = "../tests/input/phantom.png";
std::string meshFile = _settings->GetMeshFile();
Matrix gsImage = createSU2MeshFromImage( imageFile, meshFile );
auto bounds = _mesh->GetBounds();
double xMin = bounds[0].first;
double xMax = bounds[0].second;
double yMin = bounds[1].first;
double yMax = bounds[1].second;
unsigned m = gsImage.rows();
unsigned n = gsImage.columns();
Vector x( m + 1 ), y( n + 1 );
for( unsigned i = 0; i < m + 1; ++i ) {
x[i] = static_cast<double>( i ) / static_cast<double>( m ) * ( xMax - xMin );
}
for( unsigned i = 0; i < n + 1; ++i ) y[i] = static_cast<double>( i ) / static_cast<double>( n ) * ( yMax - yMin );
Interpolation interp( x, y, gsImage );
std::vector<double> result( _mesh->GetNumCells(), 0.0 );
for( unsigned i = 0; i < _mesh->GetNumCells(); ++i ) {
result[i] = std::clamp( interp( cellMidPoints[i][0], cellMidPoints[i][1] ), 0.0, 1.0 );
}
}
\ No newline at end of file
......@@ -7,6 +7,7 @@
#include "problems/musclebonelung.h"
#include "problems/problembase.h"
#include "problems/waterphantom.h"
#include "problems/phantom2d.h"
ProblemBase::ProblemBase( Config* settings, Mesh* mesh ) {
_settings = settings;
......@@ -34,6 +35,7 @@ ProblemBase* ProblemBase::Create( Config* settings, Mesh* mesh ) {
case PROBLEM_AirCavity: return new AirCavity1D( settings, mesh );
case PROBLEM_MuscleBoneLung: return new MuscleBoneLung( settings, mesh );
case PROBLEM_WaterPhantom: return new WaterPhantom( settings, mesh );
case PROBLEM_Phantom2D: return new Phantom2D( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D: return new LineSource_SN_Pseudo1D( settings, mesh );
case PROBLEM_LineSource_Pseudo_1D_Physics: return new LineSource_SN_Pseudo1D_Physics( settings, mesh );
default: return new ElectronRT( settings, mesh ); // Use RadioTherapy as dummy
......
......@@ -94,6 +94,9 @@ CSDSolverTrafoFP::CSDSolverTrafoFP( Config* settings ) : SNSolver( settings ) {
std::cout<<"];"<<std::endl;
*/
}
_density = std::vector<double> ( _nCells, 1.0 );
//exit(EXIT_SUCCESS);
}
void CSDSolverTrafoFP::Solve() {
......
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