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) 1991-2009 OpenCFD Ltd.
   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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "face.H"
  28 #include "scalarField.H"
  29
  30
  31 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
  32
  33 // Calculate area in contact given displacement of vertices relative to
  34 // the face plane. Positive displacement is above the face (no contact);
  35 // negative is in contact
  36 Foam::scalar Foam::face::areaInContact
  37 (
  38     const pointField& meshPoints,
  39     const scalarField& v
  40 ) const
  41 {
  42     // Assemble the vertex values
  43     const labelList& labels = *this;
  44
  45     scalarField vertexValue(labels.size());
  46
  47     forAll (labels, i)
  48     {
  49         vertexValue[i] = v[labels[i]];
  50     }
  51
  52
  53     // Loop through vertexValue. If all greater that 0 return 0 (no contact);
  54     // if all less than zero return 1
  55     // all zeros is assumed to be in contact.
  56
  57     bool allPositive = true;
  58     bool allNegative = true;
  59
  60     forAll (vertexValue, vI)
  61     {
  62         if (vertexValue[vI] > 0)
  63         {
  64             allNegative = false;
  65         }
  66         else
  67         {
  68             allPositive = false;
  69         }
  70     }
  71
  72     if (allPositive)
  73     {
  74         return 0.0;
  75     }
  76
  77     if (allNegative)
  78     {
  79         return 1.0;
  80     }
  81
  82     // There is a partial contact.
  83     // Algorithm:
  84     // Go through all edges. if both vertex values for the edge are
  85     // positive, discard. If one is positive and one is negative,
  86     // create a point and start the edge with it. If both are
  87     // negative, add the edge into the new face.  When finished,
  88     // calculate area of new face and return relative area (0<x<1)
  89
  90     // Dimension new point list to max possible size
  91     const labelList& faceLabels = *this;
  92
  93     pointField newFacePoints(2*size());
  94     label nNewFacePoints = 0;
  95
  96     for (label vI = 0; vI < size() - 1; vI++)
  97     {
  98         if (vertexValue[vI] <= 0)
  99         {
 100             // This is a point in contact
 101             newFacePoints[nNewFacePoints] = meshPoints[faceLabels[vI]];
 102             nNewFacePoints++;
 103         }
 104
 105         if
 106         (
 107             (vertexValue[vI] > 0 && vertexValue[vI + 1] < 0)
 108          || (vertexValue[vI] < 0 && vertexValue[vI + 1] > 0)
 109         )
 110         {
 111             // Edge intersection. Calculate intersection point and add to list
 112             point intersection =
 113                 meshPoints[faceLabels[vI]]
 114               + vertexValue[vI]/(vertexValue[vI + 1] - vertexValue[vI])
 115                 *(meshPoints[faceLabels[vI]] - meshPoints[faceLabels[vI + 1]]);
 116
 117             newFacePoints[nNewFacePoints] = intersection;
 118             nNewFacePoints++;
 119         }
 120     }
 121
 122     // Do last point by hand
 123     if (vertexValue[size() - 1] <= 0)
 124     {
 125         // This is a point in contact
 126         newFacePoints[nNewFacePoints] = meshPoints[faceLabels[size() - 1]];
 127         nNewFacePoints++;
 128     }
 129
 130     if
 131     (
 132         (vertexValue[size() - 1] > 0 && vertexValue[0] < 0)
 133      || (vertexValue[size() - 1] < 0 && vertexValue[0] > 0)
 134     )
 135     {
 136         // Edge intersection. Calculate intersection point and add to list
 137         point intersection =
 138             meshPoints[faceLabels[size() - 1]]
 139           + vertexValue[size() - 1]/(vertexValue[0] - vertexValue[size() - 1])
 140             *(meshPoints[faceLabels[size() - 1]] - meshPoints[faceLabels[0]]);
 141
 142         newFacePoints[nNewFacePoints] = intersection;
 143         nNewFacePoints++;
 144     }
 145
 146     newFacePoints.setSize(nNewFacePoints);
 147
 148     // Make a labelList for the sub-face (points are ordered!)
 149     labelList sfl(newFacePoints.size());
 150
 151     forAll (sfl, sflI)
 152     {
 153         sfl[sflI] = sflI;
 154     }
 155
 156     // Calculate relative area
 157     return face(sfl).mag(newFacePoints)/(mag(meshPoints) + VSMALL);
 158 }
 159
 160
 161 // ************************************************************************* //