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