faceAreaInContact.C
Go to the documentation of this file.
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
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
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 // ************************************************************************* //
dimensioned< scalar > mag(const dimensioned< Type > &)
scalar mag(const pointField &) const
Magnitude of face area.
Definition: faceI.H:97
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void setSize(const label)
Reset size of List.
Definition: List.C:318
#define forAll(list, i)
Definition: UList.H:421
scalar areaInContact(const pointField &, const scalarField &v) const
Return area in contact, given the displacement in vertices.
Foam::intersection.
Definition: intersection.H:49