cellPointWeight.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-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 \*---------------------------------------------------------------------------*/
25 
26 #include "cellPointWeight.H"
27 #include "polyMesh.H"
29 
30 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
31 
32 int Foam::cellPointWeight::debug(debug::debugSwitch("cellPointWeight", 0));
33 
34 Foam::scalar Foam::cellPointWeight::tol(small);
35 
36 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
37 
39 (
40  const polyMesh& mesh,
41  const vector& position,
42  const label celli
43 )
44 {
45  if (debug)
46  {
47  Pout<< nl << "Foam::cellPointWeight::findTetrahedron" << nl
48  << "position = " << position << nl
49  << "celli = " << celli << endl;
50  }
51 
52  List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
53  (
54  mesh,
55  celli
56  );
57 
58  const scalar cellVolume = mesh.cellVolumes()[celli];
59 
60  forAll(cellTets, tetI)
61  {
62  const tetIndices& tetIs = cellTets[tetI];
63 
64  // Barycentric coordinates of the position
65  scalar det = tetIs.tet(mesh).pointToBarycentric(position, weights_);
66 
67  if (mag(det/cellVolume) > tol)
68  {
69  const scalar& u = weights_[0];
70  const scalar& v = weights_[1];
71  const scalar& w = weights_[2];
72 
73  if
74  (
75  (u + tol > 0)
76  && (v + tol > 0)
77  && (w + tol > 0)
78  && (u + v + w < 1 + tol)
79  )
80  {
81 
82  faceVertices_ = tetIs.faceTriIs(mesh);
83 
84  return;
85  }
86  }
87  }
88 
89  // A suitable point in a tetrahedron was not found, find the
90  // nearest.
91 
92  scalar minNearDist = vGreat;
93 
94  label nearestTetI = -1;
95 
96  forAll(cellTets, tetI)
97  {
98  const tetIndices& tetIs = cellTets[tetI];
99 
100  scalar nearDist = tetIs.tet(mesh).nearestPoint(position).distance();
101 
102  if (nearDist < minNearDist)
103  {
104  minNearDist = nearDist;
105 
106  nearestTetI = tetI;
107  }
108  }
109 
110  if (debug)
111  {
112  Pout<< "cellPointWeight::findTetrahedron" << nl
113  << " Tetrahedron search failed; using closest tet to point "
114  << position << nl
115  << " cell: "
116  << celli << nl
117  << endl;
118  }
119 
120 
121  const tetIndices& tetIs = cellTets[nearestTetI];
122 
123  // Barycentric coordinates of the position, ignoring if the
124  // determinant is suitable. If not, the return from barycentric
125  // to weights_ is safe.
126  weights_ = tetIs.tet(mesh).pointToBarycentric(position);
127 
128  faceVertices_ = tetIs.faceTriIs(mesh);
129 }
130 
131 
133 (
134  const polyMesh& mesh,
135  const vector& position,
136  const label facei
137 )
138 {
139  if (debug)
140  {
141  Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
142  << "position = " << position << nl
143  << "facei = " << facei << endl;
144  }
145 
146  List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
147  (
148  mesh,
149  facei,
150  mesh.faceOwner()[facei]
151  );
152 
153  const scalar faceAreaSqr = sqr(mesh.magFaceAreas()[facei]);
154 
155  forAll(faceTets, tetI)
156  {
157  const tetIndices& tetIs = faceTets[tetI];
158 
159  // Barycentric coordinates of the position
160  barycentric2D triWeights;
161  const scalar det =
162  tetIs.faceTri(mesh).pointToBarycentric(position, triWeights);
163 
164  if (0.25*mag(det)/faceAreaSqr > tol)
165  {
166  const scalar& u = triWeights[0];
167  const scalar& v = triWeights[1];
168 
169  if
170  (
171  (u + tol > 0)
172  && (v + tol > 0)
173  && (u + v < 1 + tol)
174  )
175  {
176  // Weight[0] is for the cell centre.
177  weights_[0] = 0;
178  weights_[1] = triWeights[0];
179  weights_[2] = triWeights[1];
180  weights_[3] = triWeights[2];
181 
182  faceVertices_ = tetIs.faceTriIs(mesh);
183 
184  return;
185  }
186  }
187  }
188 
189  // A suitable point in a triangle was not found, find the nearest.
190 
191  scalar minNearDist = vGreat;
192 
193  label nearestTetI = -1;
194 
195  forAll(faceTets, tetI)
196  {
197  const tetIndices& tetIs = faceTets[tetI];
198 
199  scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();
200 
201  if (nearDist < minNearDist)
202  {
203  minNearDist = nearDist;
204 
205  nearestTetI = tetI;
206  }
207  }
208 
209  if (debug)
210  {
211  Pout<< "cellPointWeight::findTriangle" << nl
212  << " Triangle search failed; using closest tri to point "
213  << position << nl
214  << " face: "
215  << facei << nl
216  << endl;
217  }
218 
219  const tetIndices& tetIs = faceTets[nearestTetI];
220 
221  // Barycentric coordinates of the position, ignoring if the
222  // determinant is suitable. If not, the return from barycentric
223  // to triWeights is safe.
224 
225  const barycentric2D triWeights =
226  tetIs.faceTri(mesh).pointToBarycentric(position);
227 
228  // Weight[0] is for the cell centre.
229  weights_[0] = 0;
230  weights_[1] = triWeights[0];
231  weights_[2] = triWeights[1];
232  weights_[3] = triWeights[2];
233 
234  faceVertices_ = tetIs.faceTriIs(mesh);
235 }
236 
237 
238 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
239 
241 (
242  const polyMesh& mesh,
243  const vector& position,
244  const label celli,
245  const label facei
246 )
247 :
248  celli_(celli)
249 {
250  if (facei < 0)
251  {
252  // Face data not supplied
253  findTetrahedron(mesh, position, celli);
254  }
255  else
256  {
257  // Face data supplied
258  findTriangle(mesh, position, facei);
259  }
260 }
261 
262 
263 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
static scalar tol
Tolerance used in calculating barycentric co-ordinates.
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
barycentric2D pointToBarycentric(const point &pt) const
Calculate the barycentric coordinates from the given point.
Definition: triangleI.H:259
static int debug
Debug switch.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: HashTable.H:59
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the face for.
Definition: tetIndicesI.H:113
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
triFace faceTriIs(const polyMesh &mesh) const
Return the indices corresponding to the tri on the face for.
Definition: tetIndicesI.H:64
dimensionedScalar det(const dimensionedSphericalTensor &dt)
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet.
Definition: tetIndicesI.H:98
Templated 2D Barycentric derived from VectorSpace. Has 3 components, one of which is redundant...
Definition: Barycentric2D.H:50
const scalarField & magFaceAreas() const
cellPointWeight(const polyMesh &mesh, const vector &position, const label celli, const label facei=-1)
Construct from components.
pointHit nearestPoint(const point &p) const
Return nearest point to p on triangle.
Definition: triangleI.H:667
int debugSwitch(const char *name, const int defaultValue=0)
Lookup debug switch or add default value.
Definition: debug.C:178
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
void findTetrahedron(const polyMesh &mesh, const vector &position, const label celli)
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
static const char nl
Definition: Ostream.H:260
void findTriangle(const polyMesh &mesh, const vector &position, const label facei)
barycentric pointToBarycentric(const point &pt) const
Calculate the barycentric coordinates from the given point.
Definition: tetrahedronI.H:264
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
pointHit nearestPoint(const point &p) const
Return nearest point to p on tetrahedron. Is p itself.
Definition: tetrahedronI.H:319
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
const scalarField & cellVolumes() const