src/OpenFOAM/meshes/meshShapes/face/faceAreaInContact.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of OpenFOAM.
  10
  11     OpenFOAM is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by
  13     the Free Software Foundation, either version 3 of the License, or
  14     (at your option) any later version.
  15
  16     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
  17     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  18     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
  19     for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "face.H"
  27 #include "scalarField.H"
  28
  29
  30 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  31
  32 // Calculate area in contact given displacement of vertices relative to
  33 // the face plane. Positive displacement is above the face (no contact);
  34 // negative is in contact
  35 Foam::scalar Foam::face::areaInContact
  36 (
  37     const pointField& meshPoints,
  38     const scalarField& v
  39 ) const
  40 {
  41     // Assemble the vertex values
  42     const labelList& labels = *this;
  43
  44     scalarField vertexValue(labels.size());
  45
  46     forAll(labels, i)
  47     {
  48         vertexValue[i] = v[labels[i]];
  49     }
  50
  51
  52     // Loop through vertexValue. If all greater that 0 return 0 (no contact);
  53     // if all less than zero return 1
  54     // all zeros is assumed to be in contact.
  55
  56     bool allPositive = true;
  57     bool allNegative = true;
  58
  59     forAll(vertexValue, vI)
  60     {
  61         if (vertexValue[vI] > 0)
  62         {
  63             allNegative = false;
  64         }
  65         else
  66         {
  67             allPositive = false;
  68         }
  69     }
  70
  71     if (allPositive)
  72     {
  73         return 0.0;
  74     }
  75
  76     if (allNegative)
  77     {
  78         return 1.0;
  79     }
  80
  81     // There is a partial contact.
  82     // Algorithm:
  83     // Go through all edges. if both vertex values for the edge are
  84     // positive, discard. If one is positive and one is negative,
  85     // create a point and start the edge with it. If both are
  86     // negative, add the edge into the new face.  When finished,
  87     // calculate area of new face and return relative area (0<x<1)
  88
  89     // Dimension new point list to max possible size
  90     const labelList& faceLabels = *this;
  91
  92     pointField newFacePoints(2*size());
  93     label nNewFacePoints = 0;
  94
  95     for (label vI = 0; vI < size() - 1; vI++)
  96     {
  97         if (vertexValue[vI] <= 0)
  98         {
  99             // This is a point in contact
 100             newFacePoints[nNewFacePoints] = meshPoints[faceLabels[vI]];
 101             nNewFacePoints++;
 102         }
 103
 104         if
 105         (
 106             (vertexValue[vI] > 0 && vertexValue[vI + 1] < 0)
 107          || (vertexValue[vI] < 0 && vertexValue[vI + 1] > 0)
 108         )
 109         {
 110             // Edge intersection. Calculate intersection point and add to list
 111             point intersection =
 112                 meshPoints[faceLabels[vI]]
 113               + vertexValue[vI]/(vertexValue[vI + 1] - vertexValue[vI])
 114                 *(meshPoints[faceLabels[vI]] - meshPoints[faceLabels[vI + 1]]);
 115
 116             newFacePoints[nNewFacePoints] = intersection;
 117             nNewFacePoints++;
 118         }
 119     }
 120
 121     // Do last point by hand
 122     if (vertexValue[size() - 1] <= 0)
 123     {
 124         // This is a point in contact
 125         newFacePoints[nNewFacePoints] = meshPoints[faceLabels[size() - 1]];
 126         nNewFacePoints++;
 127     }
 128
 129     if
 130     (
 131         (vertexValue[size() - 1] > 0 && vertexValue[0] < 0)
 132      || (vertexValue[size() - 1] < 0 && vertexValue[0] > 0)
 133     )
 134     {
 135         // Edge intersection. Calculate intersection point and add to list
 136         point intersection =
 137             meshPoints[faceLabels[size() - 1]]
 138           + vertexValue[size() - 1]/(vertexValue[0] - vertexValue[size() - 1])
 139             *(meshPoints[faceLabels[size() - 1]] - meshPoints[faceLabels[0]]);
 140
 141         newFacePoints[nNewFacePoints] = intersection;
 142         nNewFacePoints++;
 143     }
 144
 145     newFacePoints.setSize(nNewFacePoints);
 146
 147     // Make a labelList for the sub-face (points are ordered!)
 148     labelList sfl(newFacePoints.size());
 149
 150     forAll(sfl, sflI)
 151     {
 152         sfl[sflI] = sflI;
 153     }
 154
 155     // Calculate relative area
 156     return face(sfl).mag(newFacePoints)/(mag(meshPoints) + VSMALL);
 157 }
 158
 159
 160 // ************************************************************************* //