FreeFOAM The Cross-Platform CFD Toolkit
Hosted by SourceForge:
Get FreeFOAM at SourceForge.net.
            Fast, secure and Free Open Source software downloads

faceAreaInContact.C

Go to the documentation of this file.
00001 /*---------------------------------------------------------------------------*\
00002   =========                 |
00003   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
00004    \\    /   O peration     |
00005     \\  /    A nd           | Copyright (C) 1991-2010 OpenCFD Ltd.
00006      \\/     M anipulation  |
00007 -------------------------------------------------------------------------------
00008 License
00009     This file is part of OpenFOAM.
00010 
00011     OpenFOAM is free software: you can redistribute it and/or modify it
00012     under the terms of the GNU General Public License as published by
00013     the Free Software Foundation, either version 3 of the License, or
00014     (at your option) any later version.
00015 
00016     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
00017     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
00019     for more details.
00020 
00021     You should have received a copy of the GNU General Public License
00022     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
00023 
00024 \*---------------------------------------------------------------------------*/
00025 
00026 #include "face.H"
00027 #include <OpenFOAM/scalarField.H>
00028 
00029 
00030 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
00031 
00032 // Calculate area in contact given displacement of vertices relative to
00033 // the face plane. Positive displacement is above the face (no contact);
00034 // negative is in contact
00035 Foam::scalar Foam::face::areaInContact
00036 (
00037     const pointField& meshPoints,
00038     const scalarField& v
00039 ) const
00040 {
00041     // Assemble the vertex values
00042     const labelList& labels = *this;
00043 
00044     scalarField vertexValue(labels.size());
00045 
00046     forAll (labels, i)
00047     {
00048         vertexValue[i] = v[labels[i]];
00049     }
00050 
00051 
00052     // Loop through vertexValue. If all greater that 0 return 0 (no contact);
00053     // if all less than zero return 1
00054     // all zeros is assumed to be in contact.
00055 
00056     bool allPositive = true;
00057     bool allNegative = true;
00058 
00059     forAll (vertexValue, vI)
00060     {
00061         if (vertexValue[vI] > 0)
00062         {
00063             allNegative = false;
00064         }
00065         else
00066         {
00067             allPositive = false;
00068         }
00069     }
00070 
00071     if (allPositive)
00072     {
00073         return 0.0;
00074     }
00075 
00076     if (allNegative)
00077     {
00078         return 1.0;
00079     }
00080 
00081     // There is a partial contact.
00082     // Algorithm:
00083     // Go through all edges. if both vertex values for the edge are
00084     // positive, discard. If one is positive and one is negative,
00085     // create a point and start the edge with it. If both are
00086     // negative, add the edge into the new face.  When finished,
00087     // calculate area of new face and return relative area (0<x<1)
00088 
00089     // Dimension new point list to max possible size
00090     const labelList& faceLabels = *this;
00091 
00092     pointField newFacePoints(2*size());
00093     label nNewFacePoints = 0;
00094 
00095     for (label vI = 0; vI < size() - 1; vI++)
00096     {
00097         if (vertexValue[vI] <= 0)
00098         {
00099             // This is a point in contact
00100             newFacePoints[nNewFacePoints] = meshPoints[faceLabels[vI]];
00101             nNewFacePoints++;
00102         }
00103 
00104         if
00105         (
00106             (vertexValue[vI] > 0 && vertexValue[vI + 1] < 0)
00107          || (vertexValue[vI] < 0 && vertexValue[vI + 1] > 0)
00108         )
00109         {
00110             // Edge intersection. Calculate intersection point and add to list
00111             point intersection =
00112                 meshPoints[faceLabels[vI]]
00113               + vertexValue[vI]/(vertexValue[vI + 1] - vertexValue[vI])
00114                 *(meshPoints[faceLabels[vI]] - meshPoints[faceLabels[vI + 1]]);
00115 
00116             newFacePoints[nNewFacePoints] = intersection;
00117             nNewFacePoints++;
00118         }
00119     }
00120 
00121     // Do last point by hand
00122     if (vertexValue[size() - 1] <= 0)
00123     {
00124         // This is a point in contact
00125         newFacePoints[nNewFacePoints] = meshPoints[faceLabels[size() - 1]];
00126         nNewFacePoints++;
00127     }
00128 
00129     if
00130     (
00131         (vertexValue[size() - 1] > 0 && vertexValue[0] < 0)
00132      || (vertexValue[size() - 1] < 0 && vertexValue[0] > 0)
00133     )
00134     {
00135         // Edge intersection. Calculate intersection point and add to list
00136         point intersection =
00137             meshPoints[faceLabels[size() - 1]]
00138           + vertexValue[size() - 1]/(vertexValue[0] - vertexValue[size() - 1])
00139             *(meshPoints[faceLabels[size() - 1]] - meshPoints[faceLabels[0]]);
00140 
00141         newFacePoints[nNewFacePoints] = intersection;
00142         nNewFacePoints++;
00143     }
00144 
00145     newFacePoints.setSize(nNewFacePoints);
00146 
00147     // Make a labelList for the sub-face (points are ordered!)
00148     labelList sfl(newFacePoints.size());
00149 
00150     forAll (sfl, sflI)
00151     {
00152         sfl[sflI] = sflI;
00153     }
00154 
00155     // Calculate relative area
00156     return face(sfl).mag(newFacePoints)/(mag(meshPoints) + VSMALL);
00157 }
00158 
00159 
00160 // ************************ vim: set sw=4 sts=4 et: ************************ //
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines