CFD Online Discussion Forums - isoSurface normal

I am not sure it is really a bug or (more likely) i am simply doing something wrong. So here are the steps i used to create the surface.

1) i created a dictionary name "alpha05Dict" and put it in /constant

Code:

FoamFile

{

version 2.0;

format ascii;

class dictionary;

location "constant";

object sampleDict;

}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //



isoField alpha1;

isoValue 0.5;

interpolate true;

//zone ABC; // Optional: zone only

//exposedPatchName fixedWalls; // Optional: zone only

//regularise false; // Optional: do not simplify

Then in createField.H i read the dictionary with

Code:

 IOdictionary alpha05Dict

    (

        IOobject

        (

            "alpha05Dict",

            runTime.constant(),

            mesh,

            IOobject::MUST_READ,

            IOobject::NO_WRITE

        )

    );

I created the surface by calling the appropiate function and then used the method update()

Code:

 sampledIsoSurface alpha05("alpha05",mesh,alpha05Dict);

 alpha05.update();

 

tmp<vectorField>  alpha05point = alpha05.points();    //points on the surface

faceList surfaceFace = alpha05.faces();      //faces on the surface

Then to calculate normals and centers i loop over the faces

Code:

  pointField pp(alpha05point );



 vectorField areavec(surfaceFace.size(),vector(0,0,0));

 vectorField centre(surfaceFace.size(),vector(0,0,0));

 vectorField areavec1(surfaceFace.size(),vector(0,0,0));

 vectorField centre1(surfaceFace.size(),vector(0,0,0));



forAll(surfaceFace,facei)

    {

      const face& ll = surfaceFace[facei];  //list of  point label of facei

      const pointField& pll = ll.points(alpha05point);



        areavec[facei] = 0.5*(pp[ll[0]] ^ pp[ll[1]]) + 0.5*(pp[ll[1]] ^ pp[ll[2]]) + 0.5*(pp[ll[2]] ^  pp[ll[0]]);

        centre[facei] = (1.0/3.0) * (pp[ll[0]] + pp[ll[1]] + pp[ll[2]]);



          areavec1[facei] = ll.normal(pll);

          centre1[facei] = ll.centre(pll);



    }

I calculated centers and vector area using standard formula for a triangle (centre, areavec) and using the in-bluid functions of the class face (centre1, areavec1). The should return the same values...
Then i pun some "Info" to print the calculated centers and normals.

Code:

Points3((0.028 0.0290354 0.0001) (0.028 0.0290354 3.53975e-06) (0.0280251 0.0290251 9.74915e-05)) //points on a particular face



areavec1(0 0 0)



areavec(-4.97389e-10 -1.20984e-09 0)



Centre1 (0 0 0)



Centre (0.0280084 0.029032 6.70104e-05)

This is just an example, but what i got is something like this for every face. .normal() and .centre() always return zero and i don't know why.
If i plot the vector area obtained with the explicit formula i have the problem mentioned in my previous post: the direction of the vectors seems more or less random (sometime inward, sometime outward) and there are vectors with non-zero z component (case is 2D). I noticed that in isoSurface.H is reported
"- does not do 2D correctly, creates non-flat iso surface.",
which probably explains why i get strange orientation with 2D. Despite this if I sum up the face area vector magnitude i get the total interfacial area with an error less than 0.002% for a fine mesh, so the triangulated isosurface seems to well represent the 0.5 contourn of the VOF function.
I simulated the same case 3D but i still have problem with the orientation of the vectors.

Any help will be really appreaciated.