Skip to content
Snippets Groups Projects
Commit 8da80bf4 authored by Simon Praetorius's avatar Simon Praetorius Committed by Ansgar Burchardt
Browse files

changed value type to double in vector and matrix in eigenvalue computations...

changed value type to double in vector and matrix in eigenvalue computations in JacobiNQuadratureRule1D
parent 5f0cbf70
No related branches found
No related tags found
No related merge requests found
......@@ -111,21 +111,18 @@ namespace Dune {
{
using std::sqrt;
typedef DynamicVector<ct> Vector;
typedef DynamicMatrix<ct> Matrix;
// compute the degree of the needed jacobi polymomial
// compute the degree of the needed jacobi polynomial
const int n = degree/2 +1;
Matrix J(n,n,0);
DynamicMatrix<double> J(n,n,0);
J[0][0] = (ct) -alpha/(2 + alpha);
J[0][0] = -double(alpha)/(2 + alpha);
for(int i=1; i<n; ++i)
{
ct a_i = (ct) (2*i*(i + alpha)) /( (2*i + alpha - 1)*(2*i + alpha) );
ct c_i = (ct) (2*i*(i + alpha)) /( (2*i + alpha + 1)*(2*i + alpha) );
ct a_i = double(2*i*(i + alpha)) /( (2*i + alpha - 1)*(2*i + alpha) );
ct c_i = double(2*i*(i + alpha)) /( (2*i + alpha + 1)*(2*i + alpha) );
J[i][i] = (ct) - alpha*alpha /( (2*i + alpha + 2)*(2*i + alpha) );
J[i][i] = -double(alpha*alpha) /( (2*i + alpha + 2)*(2*i + alpha) );
J[i][i-1] = sqrt(a_i*c_i);
J[i-1][i] = J[i][i-1];
}
......@@ -135,7 +132,7 @@ namespace Dune {
DynamicMatrixHelp::eigenValuesNonSym(J, eigenValues, &eigenVectors);
ct mu = (ct) 1/(alpha + 1);
double mu = 1.0/(alpha + 1);
QuadratureRule<ct,1> quadratureRule;
quadratureRule.reserve(n);
......@@ -143,7 +140,7 @@ namespace Dune {
{
auto&& eV0 = eigenVectors[i][0];
ct weight = mu * eV0*eV0;
Vector node(1,0.5*eigenValues[i].real() + 0.5);
DynamicVector<ct> node(1,0.5*eigenValues[i].real() + 0.5);
// bundle the nodes and the weights
QuadraturePoint<ct,1> temp(node, weight);
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment