Commit dddb70f0 authored by Ansgar Burchardt's avatar Ansgar Burchardt

[!130] [bugfix] changed value-type to double in vector and matrix in JacobiNQuadratureRule1D

Merge branch 'issue/quadrature_rules_ct' into 'master'

ref:core/dune-geometry

### Summary

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

### Details

When instantiating JacobiNQuadratureRule1D with ct=float, the function
DynamicMatrixHelp::eigenValuesNonSym is not found since float and double are
mixed in the argument types. The reason is that the matrix was filled with ct
and the vector was filled with complex<double>.

The bug is introduced in [!127]

See merge request [!130]

  [!127]: gitlab.dune-project.org/NoneNone/merge_requests/127
  [!130]: gitlab.dune-project.org/core/dune-geometry/merge_requests/130
parents 5f0cbf70 8da80bf4
Pipeline #21907 passed with stage
in 3 minutes and 15 seconds
......@@ -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);
......
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