applications/solvers/solidMechanics/elasticNonLinULSolidFoam/performEdgeCorrectedVolPointInterpolation.H

   1 // Do "normal" interpolation
   2 //volPointInterpolation::interpolate(vf, pf);
   3
   4 //- interpolation is done just before this file
   5 pointVectorField& pf = pointDU;
   6
   7
   8 // Do the correction
   9 //GeometricField<Type, pointPatchField, pointMesh>  pfCorr
  10 /*pointVectorField pfCorr
  11 (
  12  IOobject
  13  (
  14   //  "edgeCorrectedVolPointInterpolate(" + vf.name() + ")Corr",
  15   "edgeCorrectedVolPointInterpolate(" + DU.name() + ")Corr",
  16   //vf.instance(),
  17   DU,
  18   pMesh,
  19   IOobject::NO_READ,
  20   IOobject::NO_WRITE
  21   ),
  22  pMesh,
  23  //dimensioned<Type>("zero", pf.dimensions(), pTraits<Type>::zero),
  24  dimensionedVector("zero", pf.dimensions(), vector::zero),
  25  pf.boundaryField().types()
  26  );*/
  27
  28 pointVectorField pfCorr
  29 (
  30  IOobject
  31  (
  32   "pointDUcorr",
  33   runTime.timeName(),
  34   mesh,
  35   IOobject::NO_READ,
  36   IOobject::NO_WRITE
  37   ),
  38  pMesh,
  39  dimensionedVector("vector", dimLength, vector::zero),
  40  "calculated"
  41  );
  42
  43 //const labelList& ptc = boundaryPoints();
  44 #include "findBoundaryPoints.H"
  45
  46 //const FieldField<Field, vector>& extraVecs = extrapolationVectors(); ///*********
  47 #include "calculateExtrapolationVectors.H"
  48
  49 //const FieldField<Field, scalar>& w = pointBoundaryWeights(); ///*********
  50 #include "calculatePointBoundaryWeights.H"
  51
  52 //Info << "w: " << w << endl;
  53
  54 const labelListList& PointFaces = mesh.pointFaces();
  55
  56 //const fvBoundaryMesh& bm = mesh.boundary(); // declared in findBoundaryPoints.H
  57
  58 forAll (ptc, pointI)
  59 {
  60   const label curPoint = ptc[pointI];
  61
  62   const labelList& curFaces = PointFaces[curPoint];
  63
  64   label fI = 0;
  65
  66   // Go through all the faces
  67   forAll (curFaces, faceI)
  68     {
  69       if (!mesh.isInternalFace(curFaces[faceI]))
  70         {
  71           // This is a boundary face.  If not in the empty patch
  72           // or coupled calculate the extrapolation vector
  73           label patchID =
  74             mesh.boundaryMesh().whichPatch(curFaces[faceI]);
  75
  76           if
  77             (
  78              !isA<emptyFvPatch>(mesh.boundary()[patchID])
  79              && !mesh.boundary()[patchID].coupled()
  80              )
  81             {
  82               label faceInPatchID =
  83                 bm[patchID].patch().whichFace(curFaces[faceI]);
  84
  85               pfCorr[curPoint] +=
  86                 w[pointI][fI]*
  87                 (
  88                  extraVecs[pointI][fI]
  89                  & gradDU.boundaryField()[patchID][faceInPatchID]
  90                  );
  91
  92               fI++;
  93             }
  94         }
  95     }
  96 }
  97
  98 // Update coupled boundaries
  99 /*forAll (pfCorr.boundaryField(), patchI)
 100 {
 101   if (pfCorr.boundaryField()[patchI].coupled())
 102     {
 103       pfCorr.boundaryField()[patchI].initAddField();
 104     }
 105     }*/
 106
 107  /*forAll (pfCorr.boundaryField(), patchI)
 108 {
 109   if (pfCorr.boundaryField()[patchI].coupled())
 110     {
 111       pfCorr.boundaryField()[patchI].addField(pfCorr.internalField());
 112     }
 113     }*/
 114
 115
 116   //Info << "pfCorr: " << pfCorr << endl;
 117  pfCorr.correctBoundaryConditions();
 118
 119 //pfCorr.write();
 120
 121 //Info << "pf: " << pf << endl;
 122 pf += pfCorr;
 123
 124 //- this was already commented
 125 // Replace extrapolated values in pf
 126 //     forAll (ptc, pointI)
 127 //     {
 128 //         pf[ptc[pointI]] = pfCorr[ptc[pointI]];
 129 //     }
 130
 131 pf.correctBoundaryConditions();
 132
 133 //pf.write();