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 "defineDebugSwitch.H"
28 #include "polyMesh.H"
30 
31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
32 
33 namespace Foam
34 {
35  defineTypeNameAndDebug(cellPointWeight, 0);
36 }
37 
38 Foam::scalar Foam::cellPointWeight::tol(small);
39 
40 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
41 
43 (
44  const polyMesh& mesh,
45  const vector& position,
46  const label celli
47 )
48 {
49  if (debug)
50  {
51  Pout<< nl << "Foam::cellPointWeight::findTetrahedron" << nl
52  << "position = " << position << nl
53  << "celli = " << celli << endl;
54  }
55 
57  (
58  mesh,
59  celli
60  );
61 
62  const scalar cellVolume = mesh.cellVolumes()[celli];
63 
64  forAll(cellTets, tetI)
65  {
66  const tetIndices& tetIs = cellTets[tetI];
67 
68  // Barycentric coordinates of the position
69  scalar det = tetIs.tet(mesh).pointToBarycentric(position, weights_);
70 
71  if (mag(det/cellVolume) > tol)
72  {
73  const scalar& u = weights_[0];
74  const scalar& v = weights_[1];
75  const scalar& w = weights_[2];
76 
77  if
78  (
79  (u + tol > 0)
80  && (v + tol > 0)
81  && (w + tol > 0)
82  && (u + v + w < 1 + tol)
83  )
84  {
85 
86  faceVertices_ = tetIs.faceTriIs(mesh);
87 
88  return;
89  }
90  }
91  }
92 
93  // A suitable point in a tetrahedron was not found, find the
94  // nearest.
95 
96  scalar minNearDist = vGreat;
97 
98  label nearestTetI = -1;
99 
100  forAll(cellTets, tetI)
101  {
102  const tetIndices& tetIs = cellTets[tetI];
103 
104  scalar nearDist = tetIs.tet(mesh).nearestPoint(position).distance();
105 
106  if (nearDist < minNearDist)
107  {
108  minNearDist = nearDist;
109 
110  nearestTetI = tetI;
111  }
112  }
113 
114  if (debug)
115  {
116  Pout<< "cellPointWeight::findTetrahedron" << nl
117  << " Tetrahedron search failed; using closest tet to point "
118  << position << nl
119  << " cell: "
120  << celli << nl
121  << endl;
122  }
123 
124 
125  const tetIndices& tetIs = cellTets[nearestTetI];
126 
127  // Barycentric coordinates of the position, ignoring if the
128  // determinant is suitable. If not, the return from barycentric
129  // to weights_ is safe.
130  weights_ = tetIs.tet(mesh).pointToBarycentric(position);
131 
132  faceVertices_ = tetIs.faceTriIs(mesh);
133 }
134 
135 
137 (
138  const polyMesh& mesh,
139  const vector& position,
140  const label facei
141 )
142 {
143  if (debug)
144  {
145  Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
146  << "position = " << position << nl
147  << "facei = " << facei << endl;
148  }
149 
151  (
152  mesh,
153  facei,
154  mesh.faceOwner()[facei]
155  );
156 
157  const scalar faceAreaSqr = sqr(mesh.magFaceAreas()[facei]);
158 
159  forAll(faceTets, tetI)
160  {
161  const tetIndices& tetIs = faceTets[tetI];
162 
163  // Barycentric coordinates of the position
164  barycentric2D triWeights;
165  const scalar det =
166  tetIs.faceTri(mesh).pointToBarycentric(position, triWeights);
167 
168  if (0.25*mag(det)/faceAreaSqr > tol)
169  {
170  const scalar& u = triWeights[0];
171  const scalar& v = triWeights[1];
172 
173  if
174  (
175  (u + tol > 0)
176  && (v + tol > 0)
177  && (u + v < 1 + tol)
178  )
179  {
180  // Weight[0] is for the cell centre.
181  weights_[0] = 0;
182  weights_[1] = triWeights[0];
183  weights_[2] = triWeights[1];
184  weights_[3] = triWeights[2];
185 
186  faceVertices_ = tetIs.faceTriIs(mesh);
187 
188  return;
189  }
190  }
191  }
192 
193  // A suitable point in a triangle was not found, find the nearest.
194 
195  scalar minNearDist = vGreat;
196 
197  label nearestTetI = -1;
198 
199  forAll(faceTets, tetI)
200  {
201  const tetIndices& tetIs = faceTets[tetI];
202 
203  scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();
204 
205  if (nearDist < minNearDist)
206  {
207  minNearDist = nearDist;
208 
209  nearestTetI = tetI;
210  }
211  }
212 
213  if (debug)
214  {
215  Pout<< "cellPointWeight::findTriangle" << nl
216  << " Triangle search failed; using closest tri to point "
217  << position << nl
218  << " face: "
219  << facei << nl
220  << endl;
221  }
222 
223  const tetIndices& tetIs = faceTets[nearestTetI];
224 
225  // Barycentric coordinates of the position, ignoring if the
226  // determinant is suitable. If not, the return from barycentric
227  // to triWeights is safe.
228 
229  const barycentric2D triWeights =
230  tetIs.faceTri(mesh).pointToBarycentric(position);
231 
232  // Weight[0] is for the cell centre.
233  weights_[0] = 0;
234  weights_[1] = triWeights[0];
235  weights_[2] = triWeights[1];
236  weights_[3] = triWeights[2];
237 
238  faceVertices_ = tetIs.faceTriIs(mesh);
239 }
240 
241 
242 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
243 
245 (
246  const polyMesh& mesh,
247  const vector& position,
248  const label celli,
249  const label facei
250 )
251 :
252  celli_(celli)
253 {
254  if (facei < 0)
255  {
256  // Face data not supplied
257  findTetrahedron(mesh, position, celli);
258  }
259  else
260  {
261  // Face data supplied
262  findTriangle(mesh, position, facei);
263  }
264 }
265 
266 
267 // ************************************************************************* //
#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
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
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
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
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
defineTypeNameAndDebug(combustionModel, 0)
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
Macro definitions for debug switches.
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
Namespace for OpenFOAM.
const scalarField & cellVolumes() const
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.