Did you know that TwoPhaseComputeNormals3D makes the computations three times slower?! It is nice to have precise normals but not for such costs… It would be perfect if you also provided naive fast alternative.
Do you need the normal’s during the simulation? If not, you could always implement a naive fast alternative and do it in post processing.
Well, as far as I know free surface should use them to compute which interface cells come from the positive half-space and thus need to be replaced, right?
I have implemented the following functional, but I had to modify Palabos source (I do not like when I have to modify the Palabos source code directly 
).
Let me know, if you spot any mistake 
template<typename T, template<typename U> class D>
class TwoPhaseComputeNormals3D : public BoxProcessingFunctional3D {
public:
TwoPhaseComputeNormals3D(){}
virtual TwoPhaseComputeNormals3D<T,D>* clone() const {
return new TwoPhaseComputeNormals3D<T,D>(*this);
}
virtual void processGenericBlocks(Box3D domain, std::vector<AtomicBlock3D*> atomicBlocks) {
FreeSurfaceProcessorParam3D<T,D> ps(atomicBlocks);
for (plint iX=domain.x0; iX<=domain.x1; ++iX) {
for (plint iY=domain.y0; iY<=domain.y1; ++iY) {
for (plint iZ=domain.z0; iZ<=domain.z1; ++iZ) {
T xdxW = getWaterRatio(iX+1,iY,iZ,ps);
T x_dxW = getWaterRatio(iX-1,iY,iZ,ps);
T ydyW = getWaterRatio(iX,iY+1,iZ,ps);
T y_dyW = getWaterRatio(iX,iY-1,iZ,ps);
T zdzW = getWaterRatio(iX,iY,iZ+1,ps);
T z_dzW = getWaterRatio(iX,iY,iZ-1,ps);
T x = 0.5 * (x_dxW - xdxW);
T y = 0.5 * (y_dyW - ydyW);
T z = 0.5 * (z_dzW - zdzW);
T l = sqrt(x*x + y*y + z*z);
if (l < 1e-10) {
ps.setNormal(iX, iY, iZ, Array<T,3>(0,0,0));
} else {
ps.setNormal(iX, iY, iZ, Array<T,3>(x/l, y/l, z/l));
}
}}}
}
T getWaterRatio(plint ix, plint iy, plint iz, FreeSurfaceProcessorParam3D<T,D>& params) {
using namespace twoPhaseFlag;
if (params.flag(ix,iy,iz) == fluid) {
return 1.0;
} else if (params.flag(ix,iy,iz) == interface) {
return params.volumeFraction(ix,iy,iz) / params.getDensity(ix,iy,iz);
} else {
return 0;
}
}
virtual void getTypeOfModification (std::vector<modif::ModifT>& modified) const {
std::fill(modified.begin(), modified.end(), modif::nothing);
modified[0] = modif::nothing; // Fluid.
modified[1] = modif::nothing; // rhoBar.
modified[2] = modif::nothing; // j.
modified[3] = modif::nothing; // Mass.
modified[4] = modif::nothing; // Volume fraction.
modified[5] = modif::nothing; // Flag-status.
modified[6] = modif::staticVariables; // Normal.
modified[7] = modif::nothing; // Interface-lists.
modified[8] = modif::nothing; // Curvature.
modified[9] = modif::nothing; // Outside density.
}
};