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