[Alucard]

Goodmorning i'm trying to run in parallel a code I've developed.
Almost everything works fine except this part (attached below) thatit doesn't works.
vtip is small function I've implemented by myself but it's not important right now.
Can someone please tell me what it's wrong with it?
What I do it's to start for a point P that satisfies the condition 0.01<alpha1<0.99 (alpha is 1 in one region of the domain ,0 in the other and varies as an hyperbolic tangent at the interface) and in a way I calulate 2 new points Pint et Pext that are located at alpha1=0.99 and 0.5 respectively and I calculate the two corresponding values of Cl Clext et Clint.
he problem (I think) is that when i run it in parallel the points Pint et Pext doesnt' necessary belong to the same subdomain of P.
Infact i also have a problem in running in in parallel cause it seems it doesn't like the interpolation step (i tried so to use a taylor series to calculate the Cl values as:
C=C(P)+grad(C)*deltax where deltax it's the difference vector between Pint (ou Pext) and the center of the cell Pint?Pext belong to.

Thank you in advance

///////INTERPOLATION STEP///////////////

dictionary interpolationDict = mesh.solutionDict().subDict("interpolationSchemes");
autoPtr<interpolation<scalar> > ClInterpolator = interpolation<scalar>::New(interpolationDict, C_liq);


Omega=0.;

forAll (mesh.cellCentres(),cellI)
{ 
    if(alpha1[cellI]<=0.99 && alpha1[cellI]>=0.01)  
    {

         vector p = mesh.cellCentres()[cellI];
         vector N = -nHat[cellI];
         double delta_plus;
         double delta_minus;
         double delta_tot;
         double coordinate=W.value()*Foam::log((1.-alpha1[cellI])/alpha1[cellI]);
         delta_tot=-W.value()*Foam::log((1.-0.99)/0.99);

         delta_minus=delta_tot+coordinate;
         delta_plus=-coordinate;

         vector Pext = p + delta_plus*N;
         vector Pint = p - delta_minus*N; 

         vector pos_int(Pint.x(),Pint.y(),Pint.z());
         vector pos_ext(Pext.x(),Pext.y(),Pext.z()); 
       

         label cellJ = mesh.findCell(pos_int);
         label cellK = mesh.findCell(pos_ext);

         scalar Clint = ClInterpolator->interpolate(pos_int, cellJ);
         scalar Clext = ClInterpolator->interpolate(pos_ext, cellK);
         scalar deltaCl=(Clint-Clext)/(Clint*(1.-kp0.value())); 
  

       Omega[cellI]=vtip(deltaCl, kp0.value(), 2.4e-07, liqSlope.value(),20.,20.e-06, D_liq.value());

    
    }    
}