findCellPointFaceTet.H
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-2020 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 Description
25  find the tetrahedron in which the position is.
26  while searching for the tet, store the tet
27  closest to the position.
28  This way, if position is not inside any tet,
29  choose the closest one.
30 
31 \*---------------------------------------------------------------------------*/
32 
33 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
34 
35 namespace Foam
36 {
37 
38 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
39 
40 template<class Type>
41 bool interpolationCellPointFace<Type>::findTet
42 (
43  const vector& position,
44  const label nFace,
45  vector tetPoints[],
46  label tetLabelCandidate[],
47  label tetPointLabels[],
48  scalar phi[],
49  scalar phiCandidate[],
50  label& closestFace,
51  scalar& minDistance
52 ) const
53 {
54  bool foundTet = false;
55 
56  const labelList& thisFacePoints = this->mesh_.faces()[nFace];
57  tetPoints[2] = this->mesh_.faceCentres()[nFace];
58 
59  label pointi = 0;
60 
61  while (pointi < thisFacePoints.size() && !foundTet)
62  {
63  label nextPointLabel = (pointi + 1) % thisFacePoints.size();
64 
65  tetPointLabels[0] = thisFacePoints[pointi];
66  tetPointLabels[1] = thisFacePoints[nextPointLabel];
67 
68  tetPoints[0] = this->mesh_.points()[tetPointLabels[0]];
69  tetPoints[1] = this->mesh_.points()[tetPointLabels[1]];
70 
71  bool inside = true;
72  scalar dist = 0.0;
73 
74  for (label n=0; n<4; n++)
75  {
76  label p1 = (n + 1) % 4;
77  label p2 = (n + 2) % 4;
78  label p3 = (n + 3) % 4;
79 
80  vector referencePoint, faceNormal;
81  referencePoint = tetPoints[p1];
82 
83  faceNormal =
84  (tetPoints[p3] - tetPoints[p1])
85  ^ (tetPoints[p2] - tetPoints[p1]);
86 
87  faceNormal /= mag(faceNormal);
88 
89  // correct normal to point into the tet
90  vector v0 = tetPoints[n] - referencePoint;
91  scalar correct = v0 & faceNormal;
92  if (correct < 0)
93  {
94  faceNormal = -faceNormal;
95  }
96 
97  vector v1 = position - referencePoint + small*faceNormal;
98  scalar rightSide = v1 & faceNormal;
99 
100  // since normal is inwards, inside the tet
101  // is defined as all dot-products being positive
102  inside = inside && (rightSide >= 0);
103 
104  scalar phiLength = (position - referencePoint) & faceNormal;
105 
106  scalar maxLength =
107  max(vSmall, (tetPoints[n] - referencePoint) & faceNormal);
108 
109  phi[n] = phiLength/maxLength;
110 
111  dist += phi[n];
112  }
113 
114  if (!inside)
115  {
116  if (mag(dist - 1.0) < minDistance)
117  {
118  minDistance = mag(dist - 1.0);
119  closestFace = nFace;
120 
121  for (label i=0; i<4; i++)
122  {
123  phiCandidate[i] = phi[i];
124  }
125 
126  tetLabelCandidate[0] = tetPointLabels[0];
127  tetLabelCandidate[1] = tetPointLabels[1];
128  }
129  }
130 
131  foundTet = inside;
132 
133  pointi++;
134  }
135 
136  if (foundTet)
137  {
138  closestFace = nFace;
139  }
140 
141  return foundTet;
142 }
143 
144 
145 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
146 
147 } // End namespace Foam
148 
149 // ************************************************************************* //
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Namespace for OpenFOAM.