faceAreaInContact.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 Foam::scalar Foam::face::areaInContact
33 (
34  const pointField& meshPoints,
35  const scalarField& v
36 ) const
37 {
38  // Calculate area in contact given displacement of vertices relative to
39  // the face plane. Positive displacement is above the face (no contact);
40  // negative is in contact
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
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
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 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
Foam::intersection.
Definition: intersection.H:49
scalar mag(const pointField &) const
Return scalar magnitude.
Definition: faceI.H:97
scalar areaInContact(const pointField &, const scalarField &v) const
Return area in contact, given the displacement in vertices.
void setSize(const label)
Reset size of List.
Definition: List.C:281
dimensioned< scalar > mag(const dimensioned< Type > &)