help to initialize alpha1
 

Hello, i wrote a function who initializa alpha1, based in sample from funkysetFields, but it is giving me floating point error when access to the center of the face.

Any idea??

//// inicializamos alpha1 antes de comenzar los cálculos

const pointField& pp = mesh.points();

forAll(mesh.cells(),cellI) {

const cell& cFaces = mesh.cells()[cellI];
scalar temp1=0;
scalar temp2=0;
scalar temp3=0;
scalar count=0;
forAll(cFaces, i)
{
const labelList& f = mesh.faces()[cFaces[i]];
label nPoints = f.size();
scalar maxz=-1000;
scalar minz=1000;
Info<< "cara leida " << endl;
for (label id = 0; id < nPoints; id++)
{
if (pp[f[id]].component(2) >= maxz) maxz = pp[f[id]].component(2);
if (pp[f[id]].component(2) <= minz) minz = pp[f[id]].component(2);
}
scalar projecZ = maxz-minz;
Info<< "projeccion leida " << endl;
vector nn = mesh.Cf()[mesh.cells()[cellI][cFaces[i]]];
Info<< "vector centros " << nn << endl;
Info<< "calculo proj1 " << projecZ << " z face " << nn.component(2) << endl;
temp1 += (nn.component(2) + 0.5 * projecZ); Info<< "temp1 " << temp1 << endl;
temp2 += (nn.component(2) - 0.5 * projecZ); Info<< "temp2 " << temp2 << endl;
temp3 += (0.5 - nn.component(2) / (projecZ+0.00000001)); Info<< "temp3 " << temp3 << endl;
count++;
}
Info<< "leidas caras " << endl;
temp1 /= count; temp2 /= count; temp3 /= count;
if (temp1 < 0 ) alpha1[cellI] = 1;
else if (temp2 > 0) alpha1[cellI] = 0;
else alpha1[cellI] = temp3;
Info<< "calculo alpha1 " << endl;
}

Info<< "escribo alpha1 " << endl;

alpha1.write();

Now it is working perfect, the solution was;

//// initia alpha1 before loop time

const pointField& pp = mesh.points();

forAll(mesh.cells(),cellI) {

const cell& cFaces = mesh.cells()[cellI];
scalar temp1=0;
scalar temp2=0;
scalar temp3=0;
scalar count=0;
scalar alphatemp=0;
forAll(cFaces, i)
{
const labelList& f = mesh.faces()[cFaces[i]];
label nPoints = f.size();
scalar maxz=-1000;
scalar minz=1000;

for (label id = 0; id < nPoints; id++)
{
if (pp[f[id]].component(2) >= maxz) maxz = pp[f[id]].component(2);
if (pp[f[id]].component(2) <= minz) minz = pp[f[id]].component(2);
}
scalar projecZ = maxz-minz;

vector nn = mesh.Cf()[cFaces[i]];

temp1 = (nn.component(2) + 0.5 * projecZ); //Info<< "temp1 " << temp1 << endl;
temp2 = (nn.component(2) - 0.5 * projecZ); //Info<< "temp2 " << temp2 << endl;
temp3 = (0.5 - nn.component(2) / (projecZ+0.00000001)); //Info<< "temp3 " << temp3 << endl;
if (temp1 < 0 ) alphatemp += 1;
else if (temp2 > 0) alphatemp += 0;
else alphatemp += temp3;
count++;
}

alpha1[cellI] = alphatemp/count;
}

Info<< "writing alpha1 " << endl;

alpha1.write();

So if you include this file before the loop you will avoid to use funkysetfields.

Pablo