pointMVCWeight.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-2018 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 "pointMVCWeight.H"
27 #include "polyMesh.H"
28 
29 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
30 
31 namespace Foam
32 {
33  defineDebugSwitchWithName(pointMVCWeight, "pointMVCWeight", 0);
34 }
35 
36 Foam::scalar Foam::pointMVCWeight::tol(small);
37 
38 
39 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
40 
42 (
43  const Map<label>& toLocal,
44  const face& f,
45  const DynamicList<point>& u,
46  const scalarField& dist,
47  scalarField& weights
48 ) const
49 {
50  weights.setSize(toLocal.size());
51  weights = 0.0;
52 
53  scalarField theta(f.size());
54 
55  // recompute theta, the theta computed previously are not robust
56  forAll(f, j)
57  {
58  label jPlus1 = f.fcIndex(j);
59  scalar l = mag(u[j] - u[jPlus1]);
60  theta[j] = 2.0*Foam::asin(l/2.0);
61  }
62 
63  scalar sumWeight = 0;
64  forAll(f, j)
65  {
66  label pid = toLocal[f[j]];
67  label jMin1 = f.rcIndex(j);
68  weights[pid] =
69  1.0
70  / dist[pid]
71  * (Foam::tan(theta[jMin1]/2.0) + Foam::tan(theta[j]/2.0));
72  sumWeight += weights[pid];
73  }
74 
75  if (sumWeight >= tol)
76  {
77  weights /= sumWeight;
78  }
79 }
80 
81 
83 (
84  const polyMesh& mesh,
85  const labelList& toGlobal,
86  const Map<label>& toLocal,
87  const vector& position,
88  const vectorField& uVec,
89  const scalarField& dist,
90  scalarField& weights
91 ) const
92 {
93  // Loop over all triangles of all polygons of cell to compute weights
95  DynamicList<scalar> theta(100);
96  DynamicList<point> u(100);
97 
98  const Foam::cell& cFaces = mesh.cells()[cellIndex_];
99 
100  forAll(cFaces, iter)
101  {
102  label facei = cFaces[iter];
103  const face& f = mesh.faces()[facei];
104 
105  // Pout<< "face:" << facei << " at:"
106  // << pointField(mesh.points(), f)
107  // << endl;
108 
109  // Collect the uVec for the face
110  forAll(f, j)
111  {
112  u(j) = uVec[toLocal[f[j]]];
113  }
114 
115  vector v(point::zero);
116  forAll(f, j)
117  {
118  label jPlus1 = f.fcIndex(j);
119  // Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
120 
121  vector temp = u[j] ^ u[jPlus1];
122 
123  scalar magTemp = mag(temp);
124 
125  if (magTemp < vSmall)
126  {
127  continue;
128  }
129 
130  temp /= magTemp;
131 
132  // Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1]
133  // << " temp:" << temp << endl;
134 
135  scalar l = min(mag(u[j] - u[jPlus1]), 2.0);
136  scalar angle = 2.0*Foam::asin(l/2.0);
137 
138  // Pout<< " j:" << j << " l:" << l << " angle:" << angle << endl;
139 
140  v += 0.5*angle*temp;
141  }
142 
143  scalar vNorm = mag(v);
144  v /= vNorm;
145 
146  if ((v & u[0]) < 0)
147  {
148  v = -v;
149  }
150 
151  // Pout<< " v:" << v << endl;
152 
153  // angles between edges
154  forAll(f, j)
155  {
156  label jPlus1 = f.fcIndex(j);
157  // Pout<< " uj:" << u[j] << " ujPlus1:" << u[jPlus1] << endl;
158 
159  vector n0 = u[j]^v;
160  n0 /= mag(n0);
161  vector n1 = u[jPlus1]^v;
162  n1 /= mag(n1);
163 
164  scalar l = min(mag(n0 - n1), 2.0);
165  // Pout<< " l:" << l << endl;
166  alpha(j) = 2.0*Foam::asin(l/2.0);
167 
168  vector temp = n0^n1;
169  if ((temp&v) < 0.0)
170  {
171  alpha[j] = -alpha[j];
172  }
173 
174  l = min(mag(u[j] - v), 2.0);
175  // Pout<< " l:" << l << endl;
176  theta(j) = 2.0*Foam::asin(l/2.0);
177  }
178 
179 
180  bool outlierFlag = false;
181  forAll(f, j)
182  {
183  if (mag(theta[j]) < tol)
184  {
185  outlierFlag = true;
186 
187  label pid = toLocal[f[j]];
188  weights[pid] += vNorm / dist[pid];
189  break;
190  }
191  }
192 
193  if (outlierFlag)
194  {
195  continue;
196  }
197 
198  scalar sum = 0.0;
199  forAll(f, j)
200  {
201  label jMin1 = f.rcIndex(j);
202  sum +=
203  1.0
204  / Foam::tan(theta[j])
205  * (Foam::tan(alpha[j]/2.0) + Foam::tan(alpha[jMin1]/2.0));
206  }
207 
208  // The special case when x lies on the polygon, handle it using 2D mvc.
209  // In the 2D case, alpha = theta
210  if (mag(sum) < tol)
211  {
212  // Calculate weights using face vertices only
213  calcWeights(toLocal, f, u, dist, weights);
214  return;
215  }
216 
217 
218  // Normal 3D case
219  forAll(f, j)
220  {
221  label pid = toLocal[f[j]];
222  label jMin1 = f.rcIndex(j);
223  weights[pid] +=
224  vNorm
225  / sum
226  / dist[pid]
227  / Foam::sin(theta[j])
228  * (Foam::tan(alpha[j]/2.0) + Foam::tan(alpha[jMin1]/2.0));
229  }
230  }
231 
232  // normalise weights
233  scalar sumWeight = sum(weights);
234 
235  if (mag(sumWeight) < tol)
236  {
237  return;
238  }
239  weights /= sumWeight;
240 }
241 
242 
243 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
244 
246 (
247  const polyMesh& mesh,
248  const vector& position,
249  const label cellIndex,
250  const label faceIndex
251 )
252 :
253  cellIndex_((cellIndex != -1) ? cellIndex : mesh.faceOwner()[faceIndex])
254 {
255  // Addressing - face vertices to local points and vice versa
256  const labelList& toGlobal = mesh.cellPoints()[cellIndex_];
257  Map<label> toLocal(2*toGlobal.size());
258  forAll(toGlobal, i)
259  {
260  toLocal.insert(toGlobal[i], i);
261  }
262 
263 
264  // Initialise weights
265  weights_.setSize(toGlobal.size());
266  weights_ = 0.0;
267 
268 
269  // Point-to-vertex vectors and distances
270  vectorField uVec(toGlobal.size());
271  scalarField dist(toGlobal.size());
272 
273  forAll(toGlobal, pid)
274  {
275  const point& pt = mesh.points()[toGlobal[pid]];
276 
277  uVec[pid] = pt-position;
278  dist[pid] = mag(uVec[pid]);
279 
280  // Special case: point is close to vertex
281  if (dist[pid] < tol)
282  {
283  weights_[pid] = 1.0;
284  return;
285  }
286  }
287 
288  // Project onto unit sphere
289  uVec /= dist;
290 
291 
292  if (faceIndex < 0)
293  {
294  // Face data not supplied
295  calcWeights
296  (
297  mesh,
298  toGlobal,
299  toLocal,
300  position,
301  uVec,
302  dist,
303 
304  weights_
305  );
306  }
307  else
308  {
309  DynamicList<point> u(100);
310  const face& f = mesh.faces()[faceIndex];
311  // Collect the uVec for the face
312  forAll(f, j)
313  {
314  u(j) = uVec[toLocal[f[j]]];
315  }
316 
317  // Calculate weights for face only
318  calcWeights(toLocal, f, u, dist, weights_);
319  }
320 }
321 
322 
323 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
defineDebugSwitchWithName(pointMVCWeight, "pointMVCWeight", 0)
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
const cellList & cells() const
static scalar tol
Tolerance used in calculating barycentric co-ordinates.
dimensionedScalar asin(const dimensionedScalar &ds)
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1003
bool insert(const Key &, const T &newElmt)
Insert a new hashedEntry.
Definition: HashTableI.H:80
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
pointMVCWeight(const polyMesh &mesh, const vector &position, const label celli, const label facei=-1)
Construct from components.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1041
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1028
pid_t pid()
Return the PID of this process.
Definition: POSIX.C:72
dimensionedScalar sin(const dimensionedScalar &ds)
void calcWeights(const Map< label > &toLocal, const face &f, const DynamicList< point > &u, const scalarField &dist, scalarField &weights) const
Calculate weights from single face&#39;s vertices only.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
void setSize(const label)
Reset size of List.
Definition: List.C:281
A cell is defined as a list of faces with extra functionality.
Definition: cell.H:56
dimensioned< scalar > mag(const dimensioned< Type > &)
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
dimensionedScalar tan(const dimensionedScalar &ds)
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
const labelListList & cellPoints() const
Namespace for OpenFOAM.