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