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 | Copyright (C) 2011-2013 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 // * * * * * * * * * * * * Private 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 faceList& pFaces = mesh.faces();
59  const scalar cellVolume = mesh.cellVolumes()[cellI];
60 
61  forAll(cellTets, tetI)
62  {
63  const tetIndices& tetIs = cellTets[tetI];
64 
65  const face& f = pFaces[tetIs.face()];
66 
67  // Barycentric coordinates of the position
68  scalar det = tetIs.tet(mesh).barycentric(position, weights_);
69 
70  if (mag(det/cellVolume) > tol)
71  {
72  const scalar& u = weights_[0];
73  const scalar& v = weights_[1];
74  const scalar& w = weights_[2];
75 
76  if
77  (
78  (u + tol > 0)
79  && (v + tol > 0)
80  && (w + tol > 0)
81  && (u + v + w < 1 + tol)
82  )
83  {
84  faceVertices_[0] = f[tetIs.faceBasePt()];
85  faceVertices_[1] = f[tetIs.facePtA()];
86  faceVertices_[2] = f[tetIs.facePtB()];
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  const face& f = pFaces[tetIs.face()];
128 
129  // Barycentric coordinates of the position, ignoring if the
130  // determinant is suitable. If not, the return from barycentric
131  // to weights_ is safe.
132  tetIs.tet(mesh).barycentric(position, weights_);
133 
134  faceVertices_[0] = f[tetIs.faceBasePt()];
135  faceVertices_[1] = f[tetIs.facePtA()];
136  faceVertices_[2] = f[tetIs.facePtB()];
137 }
138 
139 
141 (
142  const polyMesh& mesh,
143  const vector& position,
144  const label faceI
145 )
146 {
147  if (debug)
148  {
149  Pout<< "\nbool Foam::cellPointWeight::findTriangle" << nl
150  << "position = " << position << nl
151  << "faceI = " << faceI << endl;
152  }
153 
154  List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
155  (
156  mesh,
157  faceI,
158  mesh.faceOwner()[faceI]
159  );
160 
161  const scalar faceAreaSqr = magSqr(mesh.faceAreas()[faceI]);
162 
163  const face& f = mesh.faces()[faceI];
164 
165  forAll(faceTets, tetI)
166  {
167  const tetIndices& tetIs = faceTets[tetI];
168 
169  List<scalar> triWeights(3);
170 
171  // Barycentric coordinates of the position
172  scalar det = tetIs.faceTri(mesh).barycentric(position, triWeights);
173 
174  if (0.25*mag(det)/faceAreaSqr > tol)
175  {
176  const scalar& u = triWeights[0];
177  const scalar& v = triWeights[1];
178 
179  if
180  (
181  (u + tol > 0)
182  && (v + tol > 0)
183  && (u + v < 1 + tol)
184  )
185  {
186  // Weight[0] is for the cell centre.
187  weights_[0] = 0;
188  weights_[1] = triWeights[0];
189  weights_[2] = triWeights[1];
190  weights_[3] = triWeights[2];
191 
192  faceVertices_[0] = f[tetIs.faceBasePt()];
193  faceVertices_[1] = f[tetIs.facePtA()];
194  faceVertices_[2] = f[tetIs.facePtB()];
195 
196  return;
197  }
198  }
199  }
200 
201  // A suitable point in a triangle was not found, find the nearest.
202 
203  scalar minNearDist = VGREAT;
204 
205  label nearestTetI = -1;
206 
207  forAll(faceTets, tetI)
208  {
209  const tetIndices& tetIs = faceTets[tetI];
210 
211  scalar nearDist = tetIs.faceTri(mesh).nearestPoint(position).distance();
212 
213  if (nearDist < minNearDist)
214  {
215  minNearDist = nearDist;
216 
217  nearestTetI = tetI;
218  }
219  }
220 
221  if (debug)
222  {
223  Pout<< "cellPointWeight::findTriangle" << nl
224  << " Triangle search failed; using closest tri to point "
225  << position << nl
226  << " face: "
227  << faceI << nl
228  << endl;
229  }
230 
231  const tetIndices& tetIs = faceTets[nearestTetI];
232 
233  // Barycentric coordinates of the position, ignoring if the
234  // determinant is suitable. If not, the return from barycentric
235  // to triWeights is safe.
236 
237  List<scalar> triWeights(3);
238 
239  tetIs.faceTri(mesh).barycentric(position, triWeights);
240 
241  // Weight[0] is for the cell centre.
242  weights_[0] = 0;
243  weights_[1] = triWeights[0];
244  weights_[2] = triWeights[1];
245  weights_[3] = triWeights[2];
246 
247  faceVertices_[0] = f[tetIs.faceBasePt()];
248  faceVertices_[1] = f[tetIs.facePtA()];
249  faceVertices_[2] = f[tetIs.facePtB()];
250 }
251 
252 
253 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
254 
256 (
257  const polyMesh& mesh,
258  const vector& position,
259  const label cellI,
260  const label faceI
261 )
262 :
263  cellI_(cellI),
264  weights_(4),
265  faceVertices_(3)
266 {
267  if (faceI < 0)
268  {
269  // Face data not supplied
270  findTetrahedron(mesh, position, cellI);
271  }
272  else
273  {
274  // Face data supplied
275  findTriangle(mesh, position, faceI);
276  }
277 }
278 
279 
280 // ************************************************************************* //
const vectorField & faceAreas() const
int debugSwitch(const char *name, const int defaultValue=0)
Lookup debug switch or add default value.
Definition: debug.C:165
static int debug
Debug switch.
void findTetrahedron(const polyMesh &mesh, const vector &position, const label cellI)
dimensioned< scalar > mag(const dimensioned< Type > &)
labelList f(nPoints)
label faceBasePt() const
Return the face base point.
Definition: tetIndicesI.H:42
dimensioned< scalar > magSqr(const dimensioned< Type > &)
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
scalar barycentric(const point &pt, List< scalar > &bary) const
Calculate the barycentric coordinates of the given.
Definition: tetrahedronI.H:317
tetPointRef tet(const polyMesh &mesh) const
Return the geometry corresponding to this tet from the.
Definition: tetIndicesI.H:66
pointHit nearestPoint(const point &p) const
Return nearest point to p on triangle.
Definition: triangleI.H:657
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
cellPointWeight(const polyMesh &mesh, const vector &position, const label cellI, const label faceI=-1)
Construct from components.
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
static scalar tol
Tolerance used in calculating barycentric co-ordinates.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
#define forAll(list, i)
Definition: UList.H:421
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1073
scalar distance() const
Return distance to hit.
Definition: PointHit.H:139
void findTriangle(const polyMesh &mesh, const vector &position, const label faceI)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the.
Definition: tetIndicesI.H:109
label facePtA() const
Return face point A.
Definition: tetIndicesI.H:48
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1060
pointHit nearestPoint(const point &p) const
Return nearest point to p on tetrahedron. Is p itself.
Definition: tetrahedronI.H:362
const scalarField & cellVolumes() const
Info<< "Finished reading KIVA file"<< endl;cellShapeList cellShapes(nPoints);labelList cellZoning(nPoints,-1);const cellModel &hex=*(cellModeller::lookup("hex"));labelList hexLabels(8);label activeCells=0;labelList pointMap(nPoints);forAll(pointMap, i){pointMap[i]=i;}for(label i=0;i< nPoints;i++){if(f[i] > 0.0){hexLabels[0]=i;hexLabels[1]=i1tab[i];hexLabels[2]=i3tab[i1tab[i]];hexLabels[3]=i3tab[i];hexLabels[4]=i8tab[i];hexLabels[5]=i1tab[i8tab[i]];hexLabels[6]=i3tab[i1tab[i8tab[i]]];hexLabels[7]=i3tab[i8tab[i]];cellShapes[activeCells]=cellShape(hex, hexLabels);edgeList edges=cellShapes[activeCells].edges();forAll(edges, ei){if(edges[ei].mag(points)< SMALL){label start=pointMap[edges[ei].start()];while(start!=pointMap[start]){start=pointMap[start];}label end=pointMap[edges[ei].end()];while(end!=pointMap[end]){end=pointMap[end];}label minLabel=min(start, end);pointMap[start]=pointMap[end]=minLabel;}}cellZoning[activeCells]=idreg[i];activeCells++;}}cellShapes.setSize(activeCells);cellZoning.setSize(activeCells);forAll(cellShapes, celli){cellShape &cs=cellShapes[celli];forAll(cs, i){cs[i]=pointMap[cs[i]];}cs.collapse();}label bcIDs[11]={-1, 0, 2, 4,-1, 5,-1, 6, 7, 8, 9};const label nBCs=12;const word *kivaPatchTypes[nBCs]={&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&wallPolyPatch::typeName,&symmetryPolyPatch::typeName,&wedgePolyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&polyPatch::typeName,&symmetryPolyPatch::typeName,&oldCyclicPolyPatch::typeName};enum patchTypeNames{PISTON, VALVE, LINER, CYLINDERHEAD, AXIS, WEDGE, INFLOW, OUTFLOW, PRESIN, PRESOUT, SYMMETRYPLANE, CYCLIC};const char *kivaPatchNames[nBCs]={"piston","valve","liner","cylinderHead","axis","wedge","inflow","outflow","presin","presout","symmetryPlane","cyclic"};List< SLList< face > > pFaces[nBCs]
Definition: readKivaGrid.H:235
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:73
label facePtB() const
Return face point B.
Definition: tetIndicesI.H:54
scalar barycentric(const point &pt, List< scalar > &bary) const
Calculate the barycentric coordinates of the given.
Definition: triangleI.H:268
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
label face() const
Return the face.
Definition: tetIndicesI.H:36