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 }