Commit 774d030d authored by jannick.wolters's avatar jannick.wolters
Browse files

added mesh doc comments #25; added tensorized gauss legendre quadrature #34

parent d87451bd
Pipeline #91038 failed with stages
in 22 minutes and 12 seconds
......@@ -52,18 +52,57 @@ class Mesh
std::vector<std::pair<BOUNDARY_TYPE, std::vector<unsigned>>> boundaries );
~Mesh();
unsigned GetDim() const;
unsigned GetNumCells() const;
unsigned GetNumNodes() const;
unsigned GetNumNodesPerCell() const;
inline unsigned GetDim() const { return _dim; }
inline unsigned GetNumCells() const { return _numCells; }
inline unsigned GetNumNodes() const { return _numNodes; }
inline unsigned GetNumNodesPerCell() const { return _numNodesPerCell; }
/**
* @brief Returns all node coordinates
* @return dimension: numNodes x dim
*/
const std::vector<Vector>& GetNodes() const;
/**
* @brief Returns the mid point coordinates of each cell
* @return dimension: numCells x dim
*/
const std::vector<Vector>& GetCellMidPoints() const;
/**
* @brief Returns all node IDs that construct up each cell
* @return dimension: numCells x numNodes
*/
const std::vector<std::vector<unsigned>>& GetCells() const;
/**
* @brief Returns the cell area of each cell
* @return dimension: numCells
*/
const std::vector<double>& GetCellAreas() const;
/**
* @brief Return the color/ID of the mesh partition
* @return dimension: numCells
*/
const std::vector<unsigned>& GetPartitionIDs() const;
/**
* @brief Returns the neighbor cell IDs for every cell
* @return dimension: numCells x numNodes
*/
const std::vector<std::vector<unsigned>>& GetNeighbours() const;
/**
* @brief Returns the edge length scaled normal vectors of each cell
* @return dimension: numCells x numNodes x dim
*/
const std::vector<std::vector<Vector>>& GetNormals() const;
/**
* @brief Returns the boundary enum for each cell. BOUNDARY_TYPE::NONE is the default.
* @return dimension: numCells
*/
const std::vector<BOUNDARY_TYPE>& GetBoundaryTypes() const;
/**
......
......@@ -5,6 +5,10 @@
class QGaussLegendreTensorized : public QuadratureBase
{
private:
double Pythag( const double a, const double b );
std::pair<Vector, Matrix> ComputeEigenValTriDiagMatrix( const Matrix& mat );
public:
QGaussLegendreTensorized( unsigned order );
virtual ~QGaussLegendreTensorized() {}
......
......@@ -11,5 +11,5 @@ PROBLEM = CHECKERBOARD
% ---- Boundary Conditions ----
BC_NEUMANN = ( void )
QUAD_TYPE = LEBEDEV
QUAD_ORDER = 35
QUAD_TYPE = GAUSS_LEGENDRE_TENSORIZED
QUAD_ORDER = 10
......@@ -6,7 +6,6 @@
#include "settings/config.h"
int main( int argc, char** argv ) {
MPI_Init( &argc, &argv );
std::string filename = "default.cfg";
......
......@@ -320,11 +320,6 @@ void Mesh::ComputeSlopes( unsigned nq, VectorVector& psiDerX, VectorVector& psiD
}
}
unsigned Mesh::GetDim() const { return _dim; }
unsigned Mesh::GetNumCells() const { return _numCells; }
unsigned Mesh::GetNumNodes() const { return _numNodes; }
unsigned Mesh::GetNumNodesPerCell() const { return _numNodesPerCell; }
const std::vector<Vector>& Mesh::GetNodes() const { return _nodes; }
const std::vector<Vector>& Mesh::GetCellMidPoints() const { return _cellMidPoints; }
const std::vector<std::vector<unsigned>>& Mesh::GetCells() const { return _cells; }
......
......@@ -7,13 +7,60 @@ QGaussLegendreTensorized::QGaussLegendreTensorized( unsigned order ) : Quadratur
SetConnectivity();
}
void QGaussLegendreTensorized::SetPointsAndWeights() { // TODO
// Compute Points
// Nq random points on the sphere.
_points = VectorVector( GetNq() );
void QGaussLegendreTensorized::SetPointsAndWeights() {
Vector nodes( _order ), weights( _order );
// Compute Weights
_weights = Vector( GetNq(), 4.0 * M_PI / GetNq() );
// construct companion matrix
Matrix CM( _order, _order, 0.0 );
for( unsigned i = 0; i < _order - 1; ++i ) {
CM( i + 1, i ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
CM( i, i + 1 ) = std::sqrt( 1 / ( 4 - 1 / std::pow( static_cast<double>( i + 1 ), 2 ) ) );
}
auto evSys = ComputeEigenValTriDiagMatrix( CM );
for( unsigned i = 0; i < _order; ++i ) {
if( std::fabs( evSys.first[i] ) < 1e-15 )
nodes[i] = 0;
else
nodes[i] = evSys.first[i];
weights[i] = 2 * std::pow( evSys.second( 0, i ), 2 );
}
for( unsigned i = 0; i < _order; ++i ) {
nodes[i] = ( nodes[i] + 1.0 ) * 0.5;
}
std::vector<unsigned> p( nodes.size() );
std::iota( p.begin(), p.end(), 0 );
std::sort( p.begin(), p.end(), [&]( unsigned i, unsigned j ) { return nodes[i] < nodes[j]; } );
Vector sorted_nodes( static_cast<unsigned>( p.size() ) ), sorted_weights( static_cast<unsigned>( p.size() ) );
std::transform( p.begin(), p.end(), sorted_nodes.begin(), [&]( unsigned i ) { return nodes[i]; } );
std::transform( p.begin(), p.end(), sorted_weights.begin(), [&]( unsigned i ) { return weights[i]; } );
nodes = sorted_nodes;
weights = sorted_weights;
Vector phi( 2 * _order );
for( unsigned i = 0; i < 2 * _order; ++i ) {
phi[i] = ( i + 0.5 ) * M_PI / _order;
}
unsigned range = std::floor( _order / 2.0 );
_points.resize( _nq );
for( auto& p : _points ) {
p.resize( 3 );
}
_weights.resize( _nq );
for( unsigned j = 0; j < range; ++j ) {
for( unsigned i = 0; i < 2 * _order; ++i ) {
_points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes[j] * nodes[j] ) * std::cos( phi[i] );
_points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes[j] * nodes[j] ) * std::sin( phi[i] );
_points[j * ( 2 * _order ) + i][2] = nodes[j];
_weights[j * ( 2 * _order ) + i] = 2.0 * M_PI / _order * weights[j];
}
}
}
void QGaussLegendreTensorized::SetConnectivity() { // TODO
......@@ -21,3 +68,74 @@ void QGaussLegendreTensorized::SetConnectivity() { // TODO
VectorVectorU connectivity;
_connectivity = connectivity;
}
std::pair<Vector, Matrix> QGaussLegendreTensorized::ComputeEigenValTriDiagMatrix( const Matrix& mat ) {
unsigned n = mat.rows();
Vector d( n, 0.0 ), e( n, 0.0 );
Matrix z( n, n, 0.0 );
for( unsigned i = 0; i < n; ++i ) {
d[i] = mat( i, i );
z( i, i ) = 1.0;
i == 0 ? e[i] = 0.0 : e[i] = mat( i, i - 1 );
}
int m, l, iter, i, k;
m = l = iter = i = k = 0;
double s, r, p, g, f, dd, c, b;
s = r = p = g = f = dd = c = b = 0.0;
const double eps = std::numeric_limits<double>::epsilon();
for( i = 1; i < static_cast<int>( n ); i++ ) e[i - 1] = e[i];
e[n - 1] = 0.0;
for( l = 0; l < static_cast<int>( n ); l++ ) {
iter = 0;
do {
for( m = l; m < static_cast<int>( n ) - 1; m++ ) {
dd = std::fabs( d[m] ) + std::fabs( d[m + 1] );
if( std::fabs( e[m] ) <= eps * dd ) break;
}
if( m != l ) {
if( iter++ == 30 ) ErrorMessages::Error( "Solving the tridiagonal matrix took too many iterations!", CURRENT_FUNCTION );
g = ( d[l + 1] - d[l] ) / ( 2.0 * e[l] );
r = Pythag( g, 1.0 );
g = d[m] - d[l] + e[l] / ( g + std::copysign( r, g ) );
s = c = 1.0;
p = 0.0;
for( i = m - 1; i >= l; i-- ) {
f = s * e[i];
b = c * e[i];
e[i + 1] = ( r = Pythag( f, g ) );
if( r == 0.0 ) {
d[i + 1] -= p;
e[m] = 0.0;
break;
}
s = f / r;
c = g / r;
g = d[i + 1] - p;
r = ( d[i] - g ) * s + 2.0 * c * b;
d[i + 1] = g + ( p = s * r );
g = c * r - b;
for( k = 0; k < static_cast<int>( n ); k++ ) {
f = z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 );
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) + 1 ) =
s * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) + c * f;
z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) =
c * z( static_cast<unsigned>( k ), static_cast<unsigned>( i ) ) - s * f;
}
}
if( r == 0.0 && i >= l ) continue;
d[l] -= p;
e[l] = g;
e[m] = 0.0;
}
} while( m != l );
}
return std::make_pair( d, z );
}
double QGaussLegendreTensorized::Pythag( const double a, const double b ) {
double absa = std::fabs( a ), absb = std::fabs( b );
return ( absa > absb ? absa * std::sqrt( 1.0 + ( absb / absa ) * ( absb / absa ) )
: ( absb == 0.0 ? 0.0 : absb * std::sqrt( 1.0 + ( absa / absb ) * ( absa / absb ) ) ) );
}
......@@ -4,10 +4,10 @@
#include <vector>
std::vector<QUAD_NAME> quadraturenames = { QUAD_MonteCarlo, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
std::vector<QUAD_NAME> quadraturenames = { QUAD_MonteCarlo, QUAD_GaussLegendreTensorized, QUAD_LevelSymmetric, QUAD_Lebedev, QUAD_LDFESA };
std::vector<std::vector<int>> quadratureorders = {
{ 4, 5, 6, 7 }, // Monte Carlo
{}, // Gauss Legendre not working right now
{ 4, 6, 8, 10 }, // Gauss Legendre
{ 2, 4, 6, 8, 10, 12, 14, 16, 18, 20 }, // Available Orders for LevelSymmetric
{ 3, 5, 7, 9, 11, 13, 15, 17, 19, 21, 23, 25, 27, 29, 31, 35,
41, 47, 53, 59, 65, 71, 77, 83, 89, 95, 101, 107, 113, 119, 125, 131 }, // Available orders for Lebedev
......@@ -24,7 +24,9 @@ TEST_CASE( "Quadrature weights sum to 4*pi.", "[quadrature]" ) {
bool lowAccuracyTesting = false;
for( auto quadraturename : quadraturenames ) {
lowAccuracyTesting = false;
if( quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA ) lowAccuracyTesting = true;
if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev ||
quadraturename == QUAD_LDFESA )
lowAccuracyTesting = true;
for( auto quadratureorder : quadratureorders[quadraturename] ) {
QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
......@@ -45,7 +47,9 @@ TEST_CASE( "Quadrature points are on the unit sphere.", "[quadrature]" ) {
bool lowAccuracyTesting = false;
for( auto quadraturename : quadraturenames ) {
lowAccuracyTesting = false;
if( quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA ) lowAccuracyTesting = true;
if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev ||
quadraturename == QUAD_LDFESA )
lowAccuracyTesting = true;
for( auto quadratureorder : quadratureorders[quadraturename] ) {
QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
......@@ -92,7 +96,9 @@ TEST_CASE( "Integrate a constant function.", "[quadrature]" ) {
bool lowAccuracyTesting = false;
for( auto quadraturename : quadraturenames ) {
lowAccuracyTesting = false;
if( quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev || quadraturename == QUAD_LDFESA ) lowAccuracyTesting = true;
if( quadraturename == QUAD_GaussLegendreTensorized || quadraturename == QUAD_LevelSymmetric || quadraturename == QUAD_Lebedev ||
quadraturename == QUAD_LDFESA )
lowAccuracyTesting = true;
for( auto quadratureorder : quadratureorders[quadraturename] ) {
QuadratureBase* Q = QuadratureBase::CreateQuadrature( quadraturename, quadratureorder );
......
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