findCellPointFaceTriangle.H
Go to the documentation of this file.
1 /*
2  find the triangle in which the position is,
3  when the position is known to be on the face
4 */
5 
6 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
7 
8 namespace Foam
9 {
10 
11 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
12 
13 template<class Type>
14 bool interpolationCellPointFace<Type>::findTriangle
15 (
16  const vector& position,
17  const label nFace,
18  label tetPointLabels[],
19  scalar phi[]
20 ) const
21 {
22  bool foundTriangle = false;
23  vector tetPoints[3];
24  const labelList& facePoints = this->mesh_.faces()[nFace];
25  tetPoints[2] = this->mesh_.faceCentres()[nFace];
26 
27  label pointi = 0;
28 
29  while (pointi < facePoints.size() && !foundTriangle)
30  {
31  // The triangle is constructed from face center and one
32  // face edge
33  label nextPointLabel = (pointi + 1) % facePoints.size();
34 
35  tetPointLabels[0] = facePoints[pointi];
36  tetPointLabels[1] = facePoints[nextPointLabel];
37 
38  tetPoints[0] = this->mesh_.points()[tetPointLabels[0]];
39  tetPoints[1] = this->mesh_.points()[tetPointLabels[1]];
40 
41  vector fc = (tetPoints[0] + tetPoints[1] + tetPoints[2])/3.0;
42 
43  vector newPos = position + small*(fc-position);
44 
45  // calculate triangle edge vectors and triangle face normal
46  // the 'i':th edge is opposite node i
47  vector edge[3], normal[3];
48  for (label i=0; i<3; i++)
49  {
50  label ip0 = (i+1) % 3;
51  label ipp = (i+2) % 3;
52  edge[i] = tetPoints[ipp]-tetPoints[ip0];
53  }
54 
55  vector triangleFaceNormal = edge[1] ^ edge[2];
56 
57  // calculate edge normal (pointing inwards)
58  for (label i=0; i<3; i++)
59  {
60  normal[i] = triangleFaceNormal ^ edge[i];
61  normal[i] /= mag(normal[i]) + vSmall;
62  }
63 
64  // check if position is inside triangle
65  bool inside = true;
66  for (label i=0; i<3; i++)
67  {
68  label ip = (i+1) % 3;
69  inside = inside && (((newPos - tetPoints[ip]) & edge[i]) >= 0);
70  }
71 
72  if (inside)
73  {
74  foundTriangle = true;
75 
76  // calculate phi-values
77  for (label i=0; i<3; i++)
78  {
79  label ip = (i+1) % 3;
80  scalar phiMax = max(vSmall, normal[i] & edge[ip]);
81  scalar phiLength = (position-tetPoints[ip]) & normal[i];
82  phi[i] = phiLength/phiMax;
83  }
84  }
85 
86  pointi++;
87  }
88 
89  return foundTriangle;
90 }
91 
92 
93 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
94 
95 } // End namespace Foam
96 
97 // ************************************************************************* //
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< scalar > mag(const dimensioned< Type > &)
Namespace for OpenFOAM.