Skip to content

GitLab

Commit 3555e048 authored by steffen.schotthoefer's avatar steffen.schotthoefer
Browse files

newton optimizer works. unit tests for it work. Quadrature Gauss legendre must... 

newton optimizer works. unit tests for it work. Quadrature Gauss legendre must be on full sphere. Header clean up
parent a83d955c
Pipeline #97699 passed with stages
in 39 minutes and 35 seconds
Expand all Hide whitespace changes
Inline Side-by-side
code/include/io.h
View file @ 3555e048
#ifndef IO_H
#define IO_H
#include <chrono>
#include <filesystem>
#include <iostream>
#include "settings/typedef.h"
#include <string>
#include <vector>
#include <mpi.h>
#include <omp.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
#include <vtkPointDataToCellData.h>
#include <vtkQuad.h>
#include <vtkSmartPointer.h>
#include <vtkTriangle.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridWriter.h>
#include <Python.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>
#include "mesh.h"
#include "settings/config.h"
using vtkPointsSP = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridSP = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkTriangleSP = vtkSmartPointer<vtkTriangle>;
using vtkCellArraySP = vtkSmartPointer<vtkCellArray>;
using vtkDoubleArraySP = vtkSmartPointer<vtkDoubleArray>;
using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
using vtkCellDataToPointDataSP = vtkSmartPointer<vtkCellDataToPointData>;
using vtkPointDataToCellDataSP = vtkSmartPointer<vtkPointDataToCellData>;
// Forward Declarations
class Config;
class Mesh;
// class Matrix;
void ExportVTK( const std::string fileName,
const std::vector<std::vector<std::vector<double>>>& results,
......
code/include/kernels/scatteringkernelbase.h
View file @ 3555e048
#ifndef SCATTERINGKERNELBASE_CPP
#define SCATTERINGKERNELBASE_CPP
#include "quadratures/quadraturebase.h"
#include "settings/globalconstants.h"
#include "settings/typedef.h"
// Forward Declaration
class QuadratureBase;
class ScatteringKernel
{
private:
......
code/include/mesh.h
View file @ 3555e048
#ifndef MESH_H
#define MESH_H
#include <algorithm>
#include <mpi.h>
#include <omp.h>
#include <vector>
#include "blaze/math/CompressedMatrix.h"
#include "metis.h"
#include "parmetis.h"
#include "settings/globalconstants.h"
#include "settings/typedef.h"
#include "toolboxes/errormessages.h"
#include "reconstructor.h"
#include <vector>
class Mesh
{
......@@ -132,7 +123,6 @@ class Mesh
void ReconstructSlopesS( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
void ReconstructSlopesU( unsigned nq, VectorVector& psiDerX, VectorVector& psiDerY, const VectorVector& psi ) const;
};
#endif // MESH_H
code/include/optimizers/newtonoptimizer.h
View file @ 3555e048
......@@ -23,11 +23,14 @@ class NewtonOptimizer : public OptimizerBase
grad = <mXm*eta*'(alpha*m)> */
void ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian );
double ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments );
QuadratureBase* _quadrature; /*! @brief: used quadrature */ // THis is memory doubling! Try to use a pointer.
unsigned _nq; /*! @brief: number of quadrature points */
Vector _weights; /*! @brief quadrature weights, dim(_weights) = (_nq) */
VectorVector _quadPointsSphere; /*! @brief (my,phi), dim(_quadPoints) = (_nq,2) */
double _epsilon; /*! @brief: Termination criterion for newton optimizer */
unsigned short _maxIterations; /*! @brief: Max iterations of the newton solver */
double _alpha; /*! @brief: Newton Step Size */
unsigned short _maxLineSearches; /*! @brief: Max amount of line searches for Newton Algo */
......
code/include/quadratures/qgausslegendretensorized.h
View file @ 3555e048
......@@ -16,7 +16,7 @@ class QGaussLegendreTensorized : public QuadratureBase
virtual ~QGaussLegendreTensorized() {}
inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature"; }
inline void SetNq() override { _nq = pow( GetOrder(), 2 ); }
inline void SetNq() override { _nq = 2 * pow( GetOrder(), 2 ); }
void SetPointsAndWeights() override;
void SetConnectivity() override;
......
code/include/quadratures/quadraturebase.h
View file @ 3555e048
......@@ -69,7 +69,7 @@ class QuadratureBase
Vector _weights; /*! @brief weights of the gridpoints of the quadrature */
VectorVectorU _connectivity; /*! @brief connectivity of the gripoints of the quadrature */
VectorVector _pointsSphere; /*! @brief gridpoints of the quadrature in spherical cordinates */
VectorVector _pointsSphere; /*! @brief (my,phi)gridpoints of the quadrature in spherical cordinates */
};
#endif // QUADRATURE_H
code/include/settings/config.h
View file @ 3555e048
......@@ -82,9 +82,12 @@ class Config
KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
// Optimizer
OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
double _optimizerEpsilon; /*!< @brief termination criterion epsilon for Newton Optmizer */
unsigned short _maxIterNewtonOptimizer; /*!< @brief Maximal Number of newton iterations */
OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
double _optimizerEpsilon; /*!< @brief termination criterion epsilon for Newton Optmizer */
unsigned short _newtonIter; /*!< @brief Maximal Number of newton iterations */
double _newtonStepSize; /*!< @brief Stepsize factor for newton optimizer */
unsigned short _newtonLineSearchIter; /*!< @brief Maximal Number of line search iterations for newton optimizer */
bool _newtonFastMode; /*!< @brief If true, we skip the NewtonOptimizer for quadratic entropy and assign alpha = u */
// --- Parsing Functionality and Initializing of Options ---
/*!
......@@ -228,10 +231,14 @@ class Config
bool inline IsCSD() const { return _csd; }
unsigned inline GetMaxMomentDegree() { return _maxMomentDegree; }
// Optimizer
// Optimizer
OPTIMIZER_NAME inline GetOptimizerName() const { return _entropyOptimizerName; }
double inline GetNewtonOptimizerEpsilon() const { return _optimizerEpsilon; }
unsigned inline GetMaxIterNewtonOptimizer() const { return _maxIterNewtonOptimizer; }
unsigned inline GetNewtonIter() const { return _newtonIter; }
double inline GetNewtonStepSize() const { return _newtonStepSize; }
unsigned inline GetMaxLineSearches() const { return _newtonLineSearchIter; }
bool inline GetNewtonFastMode() const { return _newtonFastMode; }
// Boundary Conditions
BOUNDARY_TYPE GetBoundaryType( std::string nameMarker ) const; /*! @brief Get Boundary Type of given marker */
......
code/include/toolboxes/textprocessingtoolbox.h
View file @ 3555e048
......@@ -8,11 +8,11 @@
#ifndef TEXTPROCESSINGTOOLBOX_H
#define TEXTPROCESSINGTOOLBOX_H
#include "settings/typedef.h"
#include <iostream>
#include <string>
#include <vector>
#include "settings/typedef.h"
namespace TextProcessingToolbox {
/*!
......
code/input/linesource.su2
View file @ 3555e048
This diff is collapsed.
code/src/io.cpp
View file @ 3555e048
......@@ -2,6 +2,42 @@
#include "toolboxes/errormessages.h"
#include "toolboxes/textprocessingtoolbox.h"
//#include <chrono>
//#include <filesystem>
#include <iostream>
//#include <string>
#include <mpi.h>
#include <omp.h>
#include <vtkCellArray.h>
#include <vtkCellData.h>
#include <vtkCellDataToPointData.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>
//#include <vtkPointDataToCellData.h>
#include <vtkQuad.h>
#include <vtkSmartPointer.h>
#include <vtkTriangle.h>
#include <vtkUnstructuredGrid.h>
#include <vtkUnstructuredGridWriter.h>
#include <Python.h>
#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
#include <numpy/arrayobject.h>
#include "mesh.h"
#include "settings/config.h"
using vtkPointsSP = vtkSmartPointer<vtkPoints>;
using vtkUnstructuredGridSP = vtkSmartPointer<vtkUnstructuredGrid>;
using vtkTriangleSP = vtkSmartPointer<vtkTriangle>;
using vtkCellArraySP = vtkSmartPointer<vtkCellArray>;
using vtkDoubleArraySP = vtkSmartPointer<vtkDoubleArray>;
using vtkUnstructuredGridWriterSP = vtkSmartPointer<vtkUnstructuredGridWriter>;
using vtkCellDataToPointDataSP = vtkSmartPointer<vtkCellDataToPointData>;
// using vtkPointDataToCellDataSP = vtkSmartPointer<vtkPointDataToCellData>;
void ExportVTK( const std::string fileName,
const std::vector<std::vector<std::vector<double>>>& results,
const std::vector<std::string> fieldNames,
......
code/src/kernels/isotropic.cpp
View file @ 3555e048
#include "kernels/isotropic.h"
#include "quadratures/quadraturebase.h"
Isotropic::Isotropic( QuadratureBase* quad ) : ScatteringKernel( quad ) {}
......
code/src/kernels/isotropic1D.cpp
View file @ 3555e048
#include "kernels/isotropic1D.h"
#include "quadratures/quadraturebase.h"
Isotropic1D::Isotropic1D( QuadratureBase* quad ) : ScatteringKernel( quad ) {}
......
code/src/main.cpp
View file @ 3555e048
......@@ -12,6 +12,7 @@
#include <string>
// ----
#include "optimizers/optimizerbase.h"
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qmontecarlo.h"
#include "solvers/sphericalharmonics.h"
......@@ -22,108 +23,8 @@ double testFunc2( double x, double y, double z ) { return x + y + z; }
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
wchar_t* program = Py_DecodeLocale( argv[0], NULL );
Py_SetProgramName( program );
// QGaussLegendreTensorized quad( 8 );
// int maxMomentDegree = 8;
//
//// test quadrature;
// std::cout << " Integration test:" << quad.Integrate( testFunc2 ) << "\n";
//
// SphericalHarmonics testBase( maxMomentDegree );
//
// double my, phi, w, resTest = 0.0;
// Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
//// 9 basis moments if degree = 2
//// unsigned nTotalMoments = moment.size();
////
// Vector results( moment.size(), 0.0 );
//
// int idx_basis = 0, idx_basisM = 0;
//
// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
// my = quad.GetPointsSphere()[idx_quad][0];
// phi = quad.GetPointsSphere()[idx_quad][1];
// w = quad.GetWeights()[idx_quad];
// moment = testBase.ComputeSphericalBasis( my, phi );
//
// double funcEval = testFunc( my, phi );
// std::cout << "f(my= " << my << ", phi= " << phi / M_PI << "*pi) = " << funcEval << " w " << w << "\n";
// resTest += w * funcEval;
//
// // std::cout << "my, phi " << my << " , " << phi << " ::\t mom5 " << moment[5] << " mom7 \t" << moment[7] << "\n";
// // results[5] += w * moment[5] * moment[5];
// // results[7] += w * moment[7] * moment[7];
// //
// for( int idx_l = 0; idx_l <= maxMomentDegree; idx_l++ ) {
// idx_basis = testBase.GlobalIdxBasis( idx_l, 0 );
//
// results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
//
// for( int idx_k = 1; idx_k <= idx_l; idx_k++ ) {
// idx_basis = testBase.GlobalIdxBasis( idx_l, idx_k );
// idx_basisM = testBase.GlobalIdxBasis( idx_l, -idx_k );
//
// results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
// results[idx_basisM] += w * moment[idx_basisM] * moment[idx_basisM];
// }
// }
//}
// std::cout << " testInt " << resTest << "\n";
// std::cout << " Squared Integration:\n " << results << "\n";
//
// Normalized integration
// Vector resultsNormal( moment.size(), 0.0 );
//
// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
// x = quad.GetPoints()[idx_quad][0];
// y = quad.GetPoints()[idx_quad][1];
// z = quad.GetPoints()[idx_quad][2];
// w = quad.GetWeights()[idx_quad];
// moment = testBase.ComputeSphericalBasis( x, y, z );
//
// for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
//
// resultsNormal[idx_sys] += w * moment[idx_sys] * moment[idx_sys] / results[idx_sys];
// // std::cout << idx_quad << ": " << results[0] << "\n";
// }
// }
// std::cout << "moment integration:\n " << resultsNormal << "\n";
// Vector results2( moment.size(), 0.0 );
//
// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
// x = quad.GetPoints()[idx_quad][0];
// y = quad.GetPoints()[idx_quad][1];
// z = quad.GetPoints()[idx_quad][2];
// w = quad.GetWeights()[idx_quad];
// moment = testBase.ComputeSphericalBasis( x, y, z );
//
// for( unsigned idx_sys = 1; idx_sys < 9; idx_sys++ ) {
// results2[idx_sys] += w * moment[idx_sys - 1] * moment[idx_sys];
// // std::cout << idx_quad << ": " << results[0] << "\n";
// }
// }
// std::cout << "moment integration:\n " << results2 << "\n";
//
// std::ofstream myfile;
// myfile.open( "assLegendre.csv" );
//
// std::vector<double> results3;
// for( double dMy = -1.0; dMy <= 1.0; dMy += 0.01 ) {
// results3 = testBase.GetAssLegendrePoly( dMy );
// myfile << dMy;
// // std::cout << dMy;
// for( unsigned idx_basis = 0; idx_basis < results3.size(); idx_basis++ ) {
// myfile << "," << results3[idx_basis];
// // std::cout << " | " << results3[idx_basis];
// }
// // std::cout << "\n";
// myfile << "\n";
//}
// myfile.close();
//
// wchar_t* program = Py_DecodeLocale( argv[0], NULL );
// Py_SetProgramName( program );
std::string filename = ParseArguments( argc, argv );
......@@ -144,7 +45,7 @@ int main( int argc, char** argv ) {
solver->Solve();
solver->Save();
if( Py_IsInitialized() ) Py_Finalize();
// if( Py_IsInitialized() ) Py_Finalize();
MPI_Finalize();
return EXIT_SUCCESS;
}
code/src/mesh.cpp
View file @ 3555e048
#include "mesh.h"
#include <algorithm>
#include <mpi.h>
#include <omp.h>
#include "metis.h"
#include "parmetis.h"
#include "toolboxes/errormessages.h"
#include "reconstructor.h"
Mesh::Mesh( std::vector<Vector> nodes,
std::vector<std::vector<unsigned>> cells,
std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries )
......
code/src/optimizers/newtonoptimizer.cpp
View file @ 3555e048
......@@ -2,21 +2,34 @@
#include "optimizers/newtonoptimizer.h"
#include "quadratures/quadraturebase.h"
#include "settings/config.h"
#include "toolboxes/errormessages.h"
NewtonOptimizer::NewtonOptimizer( Config* settings ) : OptimizerBase( settings ) {
_quadrature = QuadratureBase::CreateQuadrature( settings->GetQuadName(), settings->GetQuadOrder() );
_nq = _quadrature->GetNq();
_weights = _quadrature->GetWeights();
_quadPointsSphere = _quadrature->GetPointsSphere();
_maxIterations = settings->GetMaxIterNewtonOptimizer();
_alpha = 0.1; // settings->GetNewtonStepSize();
_maxLineSearches = 1000; // settings->GetMaxLineSearches();
_maxIterations = settings->GetNewtonIter();
_alpha = settings->GetNewtonStepSize();
_maxLineSearches = settings->GetMaxLineSearches();
_epsilon = settings->GetNewtonOptimizerEpsilon();
}
double NewtonOptimizer::ComputeObjFunc( Vector& alpha, Vector& sol, VectorVector& moments ) {
double result = 0.0;
// Integrate
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
result += _entropy->EntropyDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad];
}
result -= dot( alpha, sol );
return result;
}
void NewtonOptimizer::ComputeGradient( Vector& alpha, Vector& sol, VectorVector& moments, Vector& grad ) {
// Reset Vector
for( unsigned idx_sys; idx_sys < grad.size(); idx_sys++ ) {
for( unsigned idx_sys = 0; idx_sys < grad.size(); idx_sys++ ) {
grad[idx_sys] = 0.0;
}
......@@ -29,14 +42,18 @@ void NewtonOptimizer::ComputeGradient( Vector& alpha, Vector& sol, VectorVector&
void NewtonOptimizer::ComputeHessian( Vector& alpha, VectorVector& moments, Matrix& hessian ) {
// Reset Matrix
for( unsigned idx_Row; idx_Row < moments.size(); idx_Row++ ) {
for( unsigned idx_Col; idx_Col < moments.size(); idx_Col++ ) {
unsigned nSize = alpha.size();
for( unsigned idx_Row = 0; idx_Row < nSize; idx_Row++ ) {
for( unsigned idx_Col = 0; idx_Col < nSize; idx_Col++ ) {
hessian( idx_Row, idx_Col ) = 0.0;
// if( idx_Col == idx_Row ) hessian( idx_Row, idx_Col ) = 1.0;
}
}
// Integrate
for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
hessian +=
outer( moments[idx_quad], moments[idx_quad] ) * ( _entropy->EntropyHessianDual( dot( alpha, moments[idx_quad] ) ) * _weights[idx_quad] );
}
......@@ -50,27 +67,22 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments
* u = current "moment solution"
*/
// if we have quadratic entropy, then alpha = u;
if( _settings->GetEntropyName() == QUADRATIC ) {
// if we have quadratic entropy, then alpha = u;
if( _settings->GetEntropyName() == QUADRATIC && _settings->GetNewtonFastMode() ) {
lambda = sol;
return;
}
// Start Newton Algorithm
//
unsigned nSize = lambda.size();
// int maxRefinements = 1000;
// unsigned nTotal = _nTotalForRef[refLevel];
unsigned nSize = lambda.size();
Vector grad( nSize, 0.0 );
// check if initial guess is good enough
ComputeGradient( lambda, sol, moments, grad );
// Gradient( g, lambda, u, refLevel );
if( norm( grad ) < _settings->GetNewtonOptimizerEpsilon() ) {
if( norm( grad ) < _epsilon ) {
return;
}
......@@ -79,11 +91,11 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments
// calculate initial Hessian and gradient
Vector dlambda = -grad;
ComputeHessian( lambda, moments, H );
// invert Hessian
invert( H );
if( _maxIterations == 1 ) {
lambda = lambda - _alpha * H * grad;
return;
......@@ -110,9 +122,10 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments
ComputeGradient( lambdaNew, sol, moments, dlambdaNew );
}
// Line Search
int lineSearchCounter = 0;
// Line Search
while( norm( dlambda ) < norm( dlambdaNew ) || !std::isfinite( norm( dlambdaNew ) ) ) {
stepSize *= 0.5;
......@@ -120,23 +133,19 @@ void NewtonOptimizer::Solve( Vector& lambda, Vector& sol, VectorVector& moments
ComputeGradient( lambdaNew, sol, moments, dlambdaNew );
// Check if FONC is fullfilled
if( norm( dlambdaNew ) < _settings->GetNewtonOptimizerEpsilon() ) {
if( norm( dlambdaNew ) < _epsilon ) {
lambda = lambdaNew;
return;
}
else if( ++lineSearchCounter > _maxLineSearches ) {
//_log->error( "[closure] Newton needed too many refinement steps!" );
std::cout << "[closure] Newton needed too many refinement steps!";
exit( EXIT_FAILURE );
ErrorMessages::Error( "Newton needed too many refinement steps!", CURRENT_FUNCTION );
}
}
lambda = lambdaNew;
if( norm( dlambdaNew ) < _settings->GetNewtonOptimizerEpsilon() ) {
if( norm( dlambdaNew ) < _epsilon ) {
lambda = lambdaNew;
return;
}
}
// _log->error( "[closure] Newton did not converge!" );
std::cout << "[closure] Newton did not converge!";
exit( EXIT_FAILURE );
ErrorMessages::Error( "Newton did not converge!", CURRENT_FUNCTION );
}
code/src/problems/checkerboard.cpp
View file @ 3555e048
#include "problems/checkerboard.h"
#include "mesh.h"
#include "settings/config.h"
Checkerboard::Checkerboard( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
_physics = nullptr;
......
code/src/problems/electronrt.cpp
View file @ 3555e048
#include "problems/electronrt.h"
#include "mesh.h"
#include "settings/config.h"
ElectronRT::ElectronRT( Config* settings, Mesh* mesh ) : ProblemBase( settings, mesh ) {
// @TODO
......
code/src/problems/linesource.cpp
View file @ 3555e048
#include "problems/linesource.h"
#include "mesh.h"
#include "settings/config.h"
// ---- LineSource_SN ----
......
code/src/problems/waterphantom.cpp
View file @ 3555e048
#include "problems/waterphantom.h"
#include "mesh.h"
#include "settings/config.h"
WaterPhantom::WaterPhantom( Config* settings, Mesh* mesh ) : ElectronRT( settings, mesh ) {
// @TODO get pointer to correct physics class
......
code/src/quadratures/qgausslegendretensorized.cpp
View file @ 3555e048
......@@ -46,8 +46,8 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
phi[i] = ( i + 0.5 ) * M_PI / _order;
}
unsigned range = std::floor( _order / 2.0 ); // comment (steffen): why do we only need half of the points:
//=> In 2D we would count everything twice. (not wrong with scaling
// unsigned range = std::floor( _order / 2.0 ); // comment (steffen): why do we only need half of the points:
//=> In 2D we would count everything twice. (not wrong with scaling
// resize points and weights
_points.resize( _nq );
......@@ -62,16 +62,16 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
_weights.resize( _nq );
// transform tensorized (x,y,z)-grid to spherical grid points
for( unsigned j = 0; j < range; ++j ) {
for( unsigned j = 0; j < _order; ++j ) {
for( unsigned i = 0; i < 2 * _order; ++i ) {
_points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::cos( phi[i] );
_points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::sin( phi[i] );
_points[j * ( 2 * _order ) + i][2] = nodes1D[j];
_pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j];
_pointsSphere[j * ( 2 * _order ) + i][1] = phi[i];
_pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j]; // my
_pointsSphere[j * ( 2 * _order ) + i][1] = phi[i]; // phi
_weights[j * ( 2 * _order ) + i] = 2.0 * M_PI / _order * weights1D[j];
_weights[j * ( 2 * _order ) + i] = M_PI / _order * weights1D[j];
}
}
}
......
Markdown is supported
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

Imprint