applications/solvers/solidMechanics/elasticPlasticNonLinULSolidFoam/calculatePointBoundaryWeights.H

   1 //void volPointInterpolation::makeBoundaryWeights() const
   2 //{
   3 // const labelList& ptc = boundaryPoints();
   4
   5 // Calculate the correction vectors and weighting factors
   6 //pointBoundaryWeightsPtr_ = new FieldField<Field, scalar>(ptc.size());
   7 //FieldField<Field, scalar>& w = *pointBoundaryWeightsPtr_;
   8 FieldField<Field, scalar> w(ptc.size());
   9
  10 {
  11   const labelListList& pf = mesh.pointFaces();
  12
  13   const volVectorField& centres = mesh.C();
  14
  15   const fvBoundaryMesh& bm = mesh.boundary();
  16
  17   pointScalarField volPointSumWeights
  18     (
  19      IOobject
  20      (
  21       "volPointSumWeights",
  22       mesh.polyMesh::instance(),
  23       mesh
  24       ),
  25      pMesh,
  26      dimensionedScalar("zero", dimless, 0)
  27      );
  28
  29   forAll (ptc, pointI)
  30     {
  31       const label curPoint = ptc[pointI];
  32
  33       const labelList& curFaces = pf[curPoint];
  34
  35       //w.hook(new scalarField(curFaces.size())); //philipc no hook function
  36       w.set
  37        (
  38         pointI,
  39         new scalarField(curFaces.size())
  40         );
  41
  42       scalarField& curWeights = w[pointI];
  43
  44       label nFacesAroundPoint = 0;
  45
  46       const vector& pointLoc = mesh.points()[curPoint];
  47
  48       // Go through all the faces
  49       forAll (curFaces, faceI)
  50         {
  51           if (!mesh.isInternalFace(curFaces[faceI]))
  52             {
  53               // This is a boundary face.  If not in the empty patch
  54               // or coupled calculate the extrapolation vector
  55               label patchID =
  56                 mesh.boundaryMesh().whichPatch(curFaces[faceI]);
  57
  58               if
  59                 (
  60                  !isA<emptyFvPatch>(bm[patchID])
  61                  && !(
  62                       bm[patchID].coupled()
  63                       //&& Pstream::parRun()
  64                       //&& !mesh.parallelData().cyclicParallel()
  65                       )
  66                  )
  67                 {
  68                   curWeights[nFacesAroundPoint] =
  69                     1.0/mag
  70                     (
  71                      pointLoc
  72                      - centres.boundaryField()[patchID]
  73                      [
  74                       bm[patchID].patch().whichFace(curFaces[faceI])
  75                       ]
  76                      );
  77
  78                   nFacesAroundPoint++;
  79                 }
  80             }
  81         }
  82
  83       // Reset the sizes of the local weights
  84       curWeights.setSize(nFacesAroundPoint);
  85
  86       // Collect the sum of weights for parallel correction
  87       volPointSumWeights[curPoint] += sum(curWeights);
  88     }
  89
  90   // Do parallel correction of weights
  91
  92   // Update coupled boundaries
  93   // Work-around for cyclic parallels.
  94   /*if (Pstream::parRun() && !mesh.parallelData().cyclicParallel())
  95     {
  96     forAll (volPointSumWeights.boundaryField(), patchI)
  97     {
  98     if (volPointSumWeights.boundaryField()[patchI].coupled())
  99     {
 100     volPointSumWeights.boundaryField()[patchI].initAddField();
 101     }
 102     }
 103
 104     forAll (volPointSumWeights.boundaryField(), patchI)
 105     {
 106     if (volPointSumWeights.boundaryField()[patchI].coupled())
 107     {
 108     volPointSumWeights.boundaryField()[patchI].addField
 109     (
 110     volPointSumWeights.internalField()
 111     );
 112     }
 113     }
 114     }*/
 115
 116   // Re-scale the weights for the current point
 117   forAll (ptc, pointI)
 118     {
 119       w[pointI] /= volPointSumWeights[ptc[pointI]];
 120     }
 121 }