[doctorWho]

I am using OpenFOAM 4.x and found that applications/test/tensor make a floating point exception error.

in the below code, vector e = eigenValues(t6) makes the error

    tensor t6(1,0,-4,0,5,4,-4,4,3);
    //tensor t6(1,2,0,2,5,0,0,0,0);

    Info<< "tensor " << t6 << endl;
    vector e = eigenValues(t6);
    Info<< "eigenvalues " << e << endl;

problematic part is in the member function eigenValues(const tensor& t) of src/OpenFOAM/primitives/Tensor/tensor/tensor.C

in the code below, Q makes zero for the test matrix (1 0 -4 0 5 4 -4 4 3) which has eigen values (-3, 3, 9). So when checking mag(PPP/QQ -1), it leads to floating point exception error before going to the third conditional PPP > QQ.

The same test code works for OpenFOAM 2.3.0 but not work for OpenFOAM 4.x (also 3.0.x). I don't know much about the OpenFOAM code but it looks like a bug after version 3.0.x. Maybe Henry Weller can fix this behavior.

    // non-diagonal matrix
    else
    {
        // Coefficients of the characteristic polynmial
        // x^3 + a*x^2 + b*x + c = 0
        scalar a =
           - t.xx() - t.yy() - t.zz();

        scalar b =
            t.xx()*t.yy() + t.xx()*t.zz() + t.yy()*t.zz()
          - t.xy()*t.yx() - t.yz()*t.zy() - t.zx()*t.xz();

        scalar c =
          - t.xx()*t.yy()*t.zz()
          - t.xy()*t.yz()*t.zx() - t.xz()*t.zy()*t.yx()
          + t.xx()*t.yz()*t.zy() + t.yy()*t.zx()*t.xz() + t.zz()*t.xy()*t.yx();

        // Auxillary variables
        scalar aBy3 = a/3;

        scalar P = (a*a - 3*b)/9; // == -p_wikipedia/3
        scalar PPP = P*P*P;

        scalar Q = (2*a*a*a - 9*a*b + 27*c)/54; // == q_wikipedia/2
        scalar QQ = Q*Q;

        // Three identical roots
        if (mag(P) < SMALL && mag(Q) < SMALL)
        {
            return vector(- aBy3, - aBy3, - aBy3);
        }

        // Two identical roots and one distinct root
        else if (mag(PPP/QQ - 1) < SMALL)
        {
            scalar sqrtP = sqrt(P);
            scalar signQ = sign(Q);

            i = ii = signQ*sqrtP - aBy3;
            iii = - 2*signQ*sqrtP - aBy3;
        }