// Do "normal" interpolation
//volPointInterpolation::interpolate(vf, pf);

//- interpolation is done just before this file
pointVectorField& pf = pointDU;


// Do the correction
//GeometricField<Type, pointPatchField, pointMesh>  pfCorr
/*pointVectorField pfCorr
(
 IOobject
 (
  //  "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
  "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
  //vf.instance(),
  DU,
  pMesh,
  IOobject::NO_READ,
  IOobject::NO_WRITE
  ),
 pMesh,
 //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
 dimensionedVector("zero", pf.dimensions(), vector::zero),
 pf.boundaryField().types()
 );*/

pointVectorField pfCorr 
(
 IOobject
 (
  "pointDUcorr",
  runTime.timeName(),
  mesh,
  IOobject::NO_READ,
  IOobject::NO_WRITE
  ),
 pMesh,
 dimensionedVector("vector", dimLength, vector::zero),
 "calculated"
 );

//const labelList& ptc = boundaryPoints();
#include "findBoundaryPoints.H"

//const FieldField<Field, vector>& extraVecs = extrapolationVectors(); ///*********
#include "calculateExtrapolationVectors.H"

//const FieldField<Field, scalar>& w = pointBoundaryWeights(); ///*********
#include "calculatePointBoundaryWeights.H"

//Info << "w: " << w << endl;

const labelListList& PointFaces = mesh.pointFaces();

//const fvBoundaryMesh& bm = mesh.boundary(); // declared in findBoundaryPoints.H

forAll (ptc, pointI)
{
  const label curPoint = ptc[pointI];
  
  const labelList& curFaces = PointFaces[curPoint];
  
  label fI = 0;
  
  // Go through all the faces
  forAll (curFaces, faceI)
    {
      if (!mesh.isInternalFace(curFaces[faceI]))
	{
	  // This is a boundary face.  If not in the empty patch
	  // or coupled calculate the extrapolation vector
	  label patchID =
	    mesh.boundaryMesh().whichPatch(curFaces[faceI]);
	  
	  if
	    (
	     !isA<emptyFvPatch>(mesh.boundary()[patchID])
	     && !mesh.boundary()[patchID].coupled()
	     )
	    {
	      label faceInPatchID =
		bm[patchID].patch().whichFace(curFaces[faceI]);
	      
	      pfCorr[curPoint] +=
		w[pointI][fI]*
		(
		 extraVecs[pointI][fI]
		 & gradDU.boundaryField()[patchID][faceInPatchID]
		 );
	      
	      fI++;
	    }
	}
    }
}

// Update coupled boundaries
/*forAll (pfCorr.boundaryField(), patchI)
{
  if (pfCorr.boundaryField()[patchI].coupled())
    {
      pfCorr.boundaryField()[patchI].initAddField();
    }
    }*/

 /*forAll (pfCorr.boundaryField(), patchI)
{
  if (pfCorr.boundaryField()[patchI].coupled())
    {
      pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
    }
    }*/


  //Info << "pfCorr: " << pfCorr << endl;
 pfCorr.correctBoundaryConditions();

//pfCorr.write();

//Info << "pf: " << pf << endl;
pf += pfCorr;

//- this was already commented
// Replace extrapolated values in pf
//     forAll (ptc, pointI)
//     {
//         pf[ptc[pointI]] = pfCorr[ptc[pointI]];
//     }

pf.correctBoundaryConditions();

//pf.write();