Commit 67f0227b authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

encapsulated momentBasisIndexing in SphericalBasisClass. Added error catches...

encapsulated momentBasisIndexing in SphericalBasisClass. Added error catches for PN solver with wrong basis choice. Adapted MN Solver to monomial basis (constructor, output, solver. Scatterkernel works only for isotropic (not tested)). Adapted setup of initial conditions to changed basis. More encapsulating in solverbase. Shifting PrepareVolumeoutput from constructor to solverPreprocessing.


Former-commit-id: 70731b6e
parent a30ebff7
......@@ -82,16 +82,6 @@ class LineSource_PN : public LineSource
private:
LineSource_PN() = delete;
/**
* @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:
LineSource_PN( Config* settings, Mesh* mesh );
~LineSource_PN();
......
......@@ -74,12 +74,7 @@ class MNSolver : public Solver
void ComputeScatterMatrix();
// Helper
/*! @brief: Computes the radiative flux from the solution vector of the moment system */
void ComputeRadFlux();
/*! @brief : computes the global index of the moment corresponding to basis function (l,k)
* @param : degree l, it must hold: 0 <= l <=_nq
* @param : order k, it must hold: -l <=k <= l
* @returns : global index
*/
int GlobalIndex( int l, int k ) const;
};
#endif // MNSOLVER_H
......@@ -50,11 +50,11 @@ class Solver
std::vector<std::vector<unsigned>> _neighbors;
// slope related params
Reconstructor* _reconstructor; /*! @brief reconstructor object for high-order scheme */
unsigned _reconsOrder; /*! @brief reconstruction order (current: 1 & 2) */
VectorVector _psiDx; /*! @brief slope of solutions in X direction */
VectorVector _psiDy; /*! @brief slope of solutions in Y direction */
VectorVector _cellMidPoints; /*! @brief middle point locations of elements */
Reconstructor* _reconstructor; /*! @brief reconstructor object for high-order scheme */
unsigned _reconsOrder; /*! @brief reconstruction order (current: 1 & 2) */
VectorVector _psiDx; /*! @brief slope of solutions in X direction */
VectorVector _psiDy; /*! @brief slope of solutions in Y direction */
VectorVector _cellMidPoints; /*! @brief middle point locations of elements */
std::vector<std::vector<Vector>> _interfaceMidPoints; /*! @brief middle point locations of edges */
// Solution related members
......@@ -117,6 +117,10 @@ class Solver
void PrepareHistoryOutput();
/*! @brief Prints HistoryOutputFields to logger */
void PrintHistoryOutput( unsigned iteration );
/*! @brief Pre Solver Screen and Logger Output */
void DrawPreSolverOutput();
/*! @brief Post Solver Screen and Logger Output */
void DrawPostSolverOutput();
public:
/*! @brief Solver constructor
......
......@@ -15,7 +15,7 @@ class SphericalBase
{
public:
SphericalBase() {}
~SphericalBase() {}
virtual ~SphericalBase() {}
/*! @brief: Create a set of basis functions on the unit sphere defined in settings
* @param: Pointer to the config file
......@@ -35,7 +35,21 @@ class SphericalBase
*/
virtual Vector ComputeSphericalBasis( double x, double y, double z ) = 0;
/*! @brief : Return size of complete Basisvector */
virtual unsigned GetBasisSize() = 0;
/*! @brief: Return number of basis functions with degree equals to currDegree
* @param: currDegree must be smaller equals _LMaxDegree */
virtual unsigned GetCurrDegreeSize( unsigned currDegree ) = 0;
/*! @brief: Computes global index of basis vector depending on order k and degree l
* @param: l_degree = degree of polynomials l = 0,1,2,3,...
* @param: k_order = order of element of degree l. !ATTENTION. Requirements are different for monomials and harmonics! */
virtual unsigned GetGlobalIndexBasis( int l_degree, int k_order ) = 0;
protected:
/*! @brief: maximal (polynomial) degree of the spherical basis (this is "L" in the comments)*/
unsigned _LMaxDegree;
};
#endif // SPHERICALBASE_H
......@@ -39,13 +39,6 @@ class SphericalHarmonics : public SphericalBase
*/
Vector ComputeSphericalBasis( double x, double y, double z ) override;
/*! @brief : helper function to get the global index for given k and l of
* the basis function Y_k^l.
* @param : l_degree - current degree of basis function, 0 <= l <= L
* @param : k_order - current order of basis function, -l <= k <= l
*/
unsigned GlobalIdxBasis( unsigned l_degree, unsigned k_order );
/*! @brief : Computes an entire set of (komplex congjugate) P_l^k and stores
* it in the vector _assLegendreP
* @param : my = cos(theta) - spherical coordinate, -1 <= my <= 1
......@@ -56,10 +49,17 @@ class SphericalHarmonics : public SphericalBase
/*! @brief: Returns length of the basis, i.e. number of elements of the basis */
unsigned GetBasisSize() override;
private:
/*! @brief: maximal degree of the spherical harmonics basis (this is "L" in the comments)*/
unsigned _LMaxDegree;
/*! @brief: Return number of basis functions with degree equals to currDegree
* @param: currDegreeL = degree of polynomials that are counted */
unsigned GetCurrDegreeSize( unsigned currDegree ) override;
/*! @brief : helper function to get the global index for given k and l of
* the basis function Y_k^l.
* @param : l_degree - current degree of basis function, 0 <= l <= L
* @param : k_order - current order of basis function, -l <= k <= l */
unsigned GetGlobalIndexBasis( int l_degree, int k_order ) override;
private:
/*! @brief: coefficients for the computations of the basis
* length of _aParam, _bParam : L + (L*(L+1))/2
*/
......@@ -97,7 +97,7 @@ class SphericalHarmonics : public SphericalBase
/*! @brief: Computes the spherical harmonics basis function up to degree _LmaxDegree at
* polar coordinates (theta, psi) and stores the result in _YBasis;
* @param:
* @param: spherical coordinate phi
*/
void ComputeYBasis( const double phi );
};
......
......@@ -39,10 +39,18 @@ class SphericalMonomials : public SphericalBase
* @return: lenght of whole basis */
unsigned GetBasisSize() override;
private:
/*! @brief: maximal degree of the spherical monomial basis (this is "L" in the comments)*/
unsigned _LMaxDegree;
/*! @brief: Computes the amount of lin. independent monomials of degree currDegreeL and
* spatial dimension dim. len of a single oder: (currDegreeL + _spatialDim -1) over (currDegreeL)
* @param: currDegreeL = degree of polynomials that are counted
* @return: lenght of a single dimension */
unsigned GetCurrDegreeSize( unsigned currDegreeL ) override;
/*! @brief: Computes global index of basis vector depending on order k and degree l
* @param: l = degree of polynomials l = 0,1,2,3,...
* @param: k = order of element of degree l. 0 <=k <=GetCurrDegreeSize(l) */
unsigned GetGlobalIndexBasis( int l_degree, int k_order ) override;
private:
/*! @brief: Spatial dimension of the unit sphere (1,2,3) */
unsigned _spatialDim;
......@@ -52,11 +60,6 @@ class SphericalMonomials : public SphericalBase
*/
Vector _YBasis;
/*! @brief: Computes the amount of lin. independent monomials of degree degree and
* spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
* @return: lenght of a single dimension */
unsigned ComputeDimensionSize( unsigned degree );
/*! @brief: Function to compute factorial of n (n!) in recursive manner */
unsigned Factorial( unsigned n );
......
......@@ -34,7 +34,7 @@ CFL_NUMBER = 0.7
% Final time for simulation
TIME_FINAL = 0.3
% Maximal Moment degree
MAX_MOMENT_SOLVER = 1
MAX_MOMENT_SOLVER = 2
%
% ---- Entropy settings ----
%
......
......@@ -442,6 +442,12 @@ void Config::SetPostprocessing() {
}
}
// --- Solver setup ---
if( GetSolverName() == PN_SOLVER ) {
ErrorMessages::Error( "PN Solver only works with spherical harmonics basis.\nThis should be the default setting for option SPHERICAL_BASIS.",
CURRENT_FUNCTION );
}
// --- Output Postprocessing ---
// Volume Output Postprocessing
......
......@@ -2,7 +2,7 @@
#include "common/config.h"
#include "common/mesh.h"
#include "toolboxes/physics.h"
#include "toolboxes/sphericalbase.h"
#include <complex>
// ---- Linesource ----
......@@ -162,12 +162,6 @@ Vector LineSource_SN_Pseudo1D_Physics::GetTotalXSE( const Vector& energies ) { r
// ---- LineSource_PN ----
int LineSource_PN::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
return numIndicesPrevLevel + prevIndicesThisLevel;
}
LineSource_PN::LineSource_PN( Config* settings, Mesh* mesh ) : LineSource( settings, mesh ) {}
LineSource_PN::~LineSource_PN() {}
......@@ -184,7 +178,11 @@ std::vector<VectorVector> LineSource_PN::GetExternalSource( const Vector& /*ener
VectorVector LineSource_PN::SetupIC() {
// Compute number of equations in the system
int ntotalEquations = GlobalIndex( _settings->GetMaxMomentDegree(), _settings->GetMaxMomentDegree() ) + 1;
// In case of PN, spherical basis is per default SPHERICAL_HARMONICS
SphericalBase* tempBase = SphericalBase::Create( _settings );
unsigned ntotalEquations = tempBase->GetBasisSize();
delete tempBase; // Only temporally needed
VectorVector psi( _mesh->GetNumCells(), Vector( ntotalEquations, 0 ) ); // zero could lead to problems?
VectorVector cellMids = _mesh->GetCellMidPoints();
......
......@@ -49,9 +49,6 @@ MNSolver::MNSolver( Config* settings ) : Solver( settings ) {
// Initialize and Pre-Compute Moments at quadrature points
_moments = VectorVector( _nq, Vector( _nTotalEntries, 0.0 ) );
ComputeMoments();
// Solver output
PrepareVolumeOutput();
}
MNSolver::~MNSolver() {
......@@ -69,12 +66,6 @@ void MNSolver::ComputeScatterMatrix() {
}
}
int MNSolver::GlobalIndex( int l, int k ) const {
int numIndicesPrevLevel = l * l; // number of previous indices untill level l-1
int prevIndicesThisLevel = k + l; // number of previous indices in current level
return numIndicesPrevLevel + prevIndicesThisLevel;
}
void MNSolver::ComputeMoments() {
double my, phi;
......@@ -191,16 +182,13 @@ void MNSolver::FVMUpdate( unsigned idx_energy ) {
for( unsigned idx_cell = 0; idx_cell < _nCells; idx_cell++ ) {
// Dirichlet Boundaries stay
if( _boundaryCells[idx_cell] == BOUNDARY_TYPE::DIRICHLET ) continue;
for( unsigned idx_system = 0; idx_system < _nTotalEntries; idx_system++ ) {
_solNew[idx_cell][idx_system] = _sol[idx_cell][idx_system] -
( _dE / _areas[idx_cell] ) * _solNew[idx_cell][idx_system] /* cell averaged flux */
- _dE * _sol[idx_cell][idx_system] *
( _sigmaT[idx_energy][idx_cell] /* absorbtion influence */
+ _sigmaS[idx_energy][idx_cell] * _scatterMatDiag[idx_system] ); /* scattering influence */
}
_solNew[idx_cell][0] += _dE * _Q[0][idx_cell][0];
}
}
......@@ -231,12 +219,23 @@ void MNSolver::PrepareVolumeOutput() {
_outputFields[idx_group].resize( _nTotalEntries );
_outputFieldNames[idx_group].resize( _nTotalEntries );
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
}
}
}
else {
for( unsigned idx_l = 0; idx_l <= _LMaxDegree; idx_l++ ) {
unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
std::string( "u_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
}
}
}
break;
......@@ -246,12 +245,23 @@ void MNSolver::PrepareVolumeOutput() {
_outputFields[idx_group].resize( _nTotalEntries );
_outputFieldNames[idx_group].resize( _nTotalEntries );
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][GlobalIndex( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][GlobalIndex( idx_l, idx_k )] =
std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
if( _settings->GetSphericalBasisName() == SPHERICAL_HARMONICS ) {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
for( int idx_k = -idx_l; idx_k <= idx_l; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
}
}
}
else {
for( int idx_l = 0; idx_l <= (int)_LMaxDegree; idx_l++ ) {
unsigned maxOrder_k = _basis->GetCurrDegreeSize( idx_l );
for( unsigned idx_k = 0; idx_k < maxOrder_k; idx_k++ ) {
_outputFields[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )].resize( _nCells );
_outputFieldNames[idx_group][_basis->GetGlobalIndexBasis( idx_l, idx_k )] =
std::string( "alpha_" + std::to_string( idx_l ) + "^" + std::to_string( idx_k ) );
}
}
}
break;
......
......@@ -51,9 +51,6 @@ PNSolver::PNSolver( Config* settings ) : Solver( settings ) {
// Compute moments of initial condition
// TODO
// Solver output
PrepareVolumeOutput();
}
void PNSolver::IterPreprocessing() {
......@@ -80,8 +77,8 @@ void PNSolver::FluxUpdate() {
_mesh->ReconstructSlopesU( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
//_mesh->ComputeSlopes( _nTotalEntries, _solDx, _solDy, _sol ); // unstructured reconstruction
}
//Vector solL( _nTotalEntries );
//Vector solR( _nTotalEntries );
// Vector solL( _nTotalEntries );
// Vector solR( _nTotalEntries );
auto solL = _sol[2];
auto solR = _sol[2];
......
......@@ -21,9 +21,6 @@ SNSolver::SNSolver( Config* settings ) : Solver( settings ) {
ScatteringKernel* k = ScatteringKernel::CreateScatteringKernel( settings->GetKernelName(), _quadrature );
_scatteringKernel = k->GetScatteringKernel();
delete k;
// Solver output
PrepareVolumeOutput();
}
void SNSolver::IterPreprocessing() {
......
......@@ -101,46 +101,12 @@ Solver* Solver::Create( Config* settings ) {
}
void Solver::Solve() {
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
auto log = spdlog::get( "event" );
auto logCSV = spdlog::get( "tabular" );
std::string hLine = "--";
if( rank == 0 ) {
unsigned strLen = 10; // max width of one column
char paddingChar = ' ';
// Assemble Header for Screen Output
std::string lineToPrint = "| ";
std::string tmpLine = "------------";
for( unsigned idxFields = 0; idxFields < _settings->GetNScreenOutput(); idxFields++ ) {
std::string tmp = _screenOutputFieldNames[idxFields];
// --- Preprocessing ---
if( strLen > tmp.size() ) // Padding
tmp.insert( 0, strLen - tmp.size(), paddingChar );
else if( strLen < tmp.size() ) // Cutting
tmp.resize( strLen );
PrepareVolumeOutput();
lineToPrint += tmp + " |";
hLine += tmpLine;
}
log->info( "---------------------------- Solver Starts -----------------------------" );
log->info( "| The simulation will run for {} iterations.", _nEnergies );
log->info( hLine );
log->info( lineToPrint );
log->info( hLine );
std::string lineToPrintCSV = "";
for( int idxFields = 0; idxFields < _settings->GetNHistoryOutput() - 1; idxFields++ ) {
std::string tmp = _historyOutputFieldNames[idxFields];
lineToPrintCSV += tmp + ",";
}
lineToPrintCSV += _historyOutputFieldNames[_settings->GetNHistoryOutput() - 1];
logCSV->info( lineToPrintCSV );
}
DrawPreSolverOutput();
// Loop over energies (pseudo-time of continuous slowing down approach)
for( unsigned iter = 0; iter < _nEnergies; iter++ ) {
......@@ -154,7 +120,7 @@ void Solver::Solve() {
// --- Finite Volume Update ---
FVMUpdate( iter );
// --- Postprocessing ---
// --- Iter Postprocessing ---
IterPostprocessing();
// --- Solver Output ---
......@@ -164,11 +130,10 @@ void Solver::Solve() {
PrintHistoryOutput( iter );
PrintVolumeOutput( iter );
}
if( rank == 0 ) {
log->info( hLine );
log->info( "| Postprocessing screen output goes here." );
log->info( "--------------------------- Solver Finished ----------------------------" );
}
// --- Postprocessing ---
DrawPostSolverOutput();
}
void Solver::PrintVolumeOutput() const { ExportVTK( _settings->GetOutputFile(), _outputFields, _outputFieldNames, _mesh ); }
......@@ -427,3 +392,82 @@ void Solver::PrintHistoryOutput( unsigned iteration ) {
}
}
}
void Solver::DrawPreSolverOutput() {
// MPI
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
// Logger
auto log = spdlog::get( "event" );
auto logCSV = spdlog::get( "tabular" );
std::string hLine = "--";
if( rank == 0 ) {
unsigned strLen = 10; // max width of one column
char paddingChar = ' ';
// Assemble Header for Screen Output
std::string lineToPrint = "| ";
std::string tmpLine = "------------";
for( unsigned idxFields = 0; idxFields < _settings->GetNScreenOutput(); idxFields++ ) {
std::string tmp = _screenOutputFieldNames[idxFields];
if( strLen > tmp.size() ) // Padding
tmp.insert( 0, strLen - tmp.size(), paddingChar );
else if( strLen < tmp.size() ) // Cutting
tmp.resize( strLen );
lineToPrint += tmp + " |";
hLine += tmpLine;
}
log->info( "---------------------------- Solver Starts -----------------------------" );
log->info( "| The simulation will run for {} iterations.", _nEnergies );
log->info( hLine );
log->info( lineToPrint );
log->info( hLine );
std::string lineToPrintCSV = "";
for( int idxFields = 0; idxFields < _settings->GetNHistoryOutput() - 1; idxFields++ ) {
std::string tmp = _historyOutputFieldNames[idxFields];
lineToPrintCSV += tmp + ",";
}
lineToPrintCSV += _historyOutputFieldNames[_settings->GetNHistoryOutput() - 1];
logCSV->info( lineToPrintCSV );
}
}
void Solver::DrawPostSolverOutput() {
// MPI
int rank;
MPI_Comm_rank( MPI_COMM_WORLD, &rank );
// Logger
auto log = spdlog::get( "event" );
std::string hLine = "--";
if( rank == 0 ) {
unsigned strLen = 10; // max width of one column
char paddingChar = ' ';
// Assemble Header for Screen Output
std::string lineToPrint = "| ";
std::string tmpLine = "------------";
for( unsigned idxFields = 0; idxFields < _settings->GetNScreenOutput(); idxFields++ ) {
std::string tmp = _screenOutputFieldNames[idxFields];
if( strLen > tmp.size() ) // Padding
tmp.insert( 0, strLen - tmp.size(), paddingChar );
else if( strLen < tmp.size() ) // Cutting
tmp.resize( strLen );
lineToPrint += tmp + " |";
hLine += tmpLine;
}
log->info( hLine );
log->info( "| Postprocessing screen output goes here." );
log->info( "--------------------------- Solver Finished ----------------------------" );
}
}
......@@ -4,9 +4,9 @@
* @author S. Schotthöfer
*/
#include "toolboxes/sphericalbase.h"
#include "common/config.h"
#include "toolboxes/errormessages.h"
#include "toolboxes/sphericalbase.h"
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/sphericalmonomials.h"
......@@ -16,7 +16,7 @@ SphericalBase* SphericalBase::Create( Config* settings ) {
switch( name ) {
case SPHERICAL_HARMONICS: return new SphericalHarmonics( maxMomentDegree );
case SPHERICAL_MONIMIALS: return new SphericalMonomials( maxMomentDegree );
case SPHERICAL_MONOMIALS: return new SphericalMonomials( maxMomentDegree );
default: ErrorMessages::Error( "Creator for the chosen basis does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
}
......
#include "toolboxes/sphericalharmonics.h"
#include "toolboxes/errormessages.h"
#include <math.h>
SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
......@@ -15,7 +16,18 @@ SphericalHarmonics::SphericalHarmonics( unsigned L_degree ) {
_YBasis = Vector( GetBasisSize(), 0.0 );
}
unsigned SphericalHarmonics::GetBasisSize() { return GlobalIdxBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */ }
unsigned SphericalHarmonics::GetBasisSize() { return GetGlobalIndexBasis( _LMaxDegree, _LMaxDegree ) + 1; /* +1, since globalIdx computes indices */ }
unsigned SphericalHarmonics::GetCurrDegreeSize( unsigned currDegreeL ) { return 2 * currDegreeL + 1; }
unsigned SphericalHarmonics::GetGlobalIndexBasis( int l_degree, int k_order ) {
if( l_degree < 0 ) ErrorMessages::Error( "Negative polynomial degrees not supported.", CURRENT_FUNCTION );
if( k_order < -l_degree || k_order > l_degree )
ErrorMessages::Error( "Order k of spherical harmonics basis must be in bounds -l<= k <= l.", CURRENT_FUNCTION );
return (unsigned)( k_order + l_degree /*number of previous indices in current level*/ +
l_degree * l_degree ) /* number of previous indices untill level l-1 */;
}
Vector SphericalHarmonics::ComputeSphericalBasis( double my, double phi ) {
ComputeAssLegendrePoly( my );
......@@ -39,8 +51,6 @@ Vector SphericalHarmonics::ComputeSphericalBasis( double x, double y, double z )
return _YBasis;
}
unsigned SphericalHarmonics::GlobalIdxBasis( unsigned l_degree, unsigned k_order ) { return k_order + l_degree + l_degree * l_degree; }
std::vector<double> SphericalHarmonics::GetAssLegendrePoly( const double my ) {
ComputeAssLegendrePoly( my );
return _assLegendreP;
......@@ -95,7 +105,7 @@ void SphericalHarmonics::ComputeAssLegendrePoly( const double my ) {
void SphericalHarmonics::ComputeYBasis( const double phi ) {
for( unsigned l_idx = 0; l_idx <= _LMaxDegree; l_idx++ ) {
_YBasis[GlobalIdxBasis( l_idx, 0 )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, 0 )] * 0.5 * M_SQRT2; // M_SQRT2 = sqrt(2)
_YBasis[GetGlobalIndexBasis( l_idx, 0 )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, 0 )] * 0.5 * M_SQRT2; // M_SQRT2 = sqrt(2)
}
// helper constants
......@@ -116,8 +126,8 @@ void SphericalHarmonics::ComputeYBasis( const double phi ) {
c2 = c1;
c1 = c;
for( unsigned l_idx = k_idx; l_idx <= _LMaxDegree; l_idx++ ) {
_YBasis[GlobalIdxBasis( l_idx, -k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * s;
_YBasis[GlobalIdxBasis( l_idx, k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * c;
_YBasis[GetGlobalIndexBasis( l_idx, -k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * s;
_YBasis[GetGlobalIndexBasis( l_idx, k_idx )] = _assLegendreP[GlobalIdxAssLegendreP( l_idx, k_idx )] * c;
}
}
}
......@@ -5,6 +5,7 @@
*/
#include "toolboxes/sphericalmonomials.h"
#include "toolboxes/errormessages.h"
SphericalMonomials::SphericalMonomials( unsigned L_degree ) {
_LMaxDegree = L_degree;
......@@ -48,13 +49,25 @@ Vector SphericalMonomials::ComputeSphericalBasis( double x, double y, double z )
unsigned SphericalMonomials::GetBasisSize() {
unsigned basisLen = 0;
for( unsigned idx_degree = 0; idx_degree <= _LMaxDegree; idx_degree++ ) {