MPLICface.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) 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 "MPLICface.H"
27 
28 
29 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
30 
31 Foam::MPLICface::MPLICface(const bool unweighted)
32 :
33  unweighted_(unweighted)
34 {};
35 
36 
37 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
38 
40 (
41  const labelList& f,
42  const labelList& faceEdges,
43  const pointField& points,
44  const boolList& isEdgeCutOld,
45  boolList& isEdgeCut,
46  label& faceEdgei,
47  const UList<scalar>& pointsAlpha,
48  const UList<vector>& pointsU,
49  const label facei,
50  const scalar target,
51  const bool ow
52 )
53 {
54  // Clear all the storage
55  cutPoints_.clear();
56  subPoints_.clear();
57  cutEdges_.clear();
58  subPointsU_.clear();
59 
60  // Direction of face circulators
61  flipped_ = false;
62 
63  label fp;
64  if (faceEdgei == -1)
65  {
66  const UIndirectList<scalar> fAlphas(pointsAlpha, f);
67 
68  // Starting from this point on the face
69  fp = findMin(fAlphas);
70  flipped_ = ow ? false : true;
71  }
72  else
73  {
74  // Finding first index of the edge in the face point list
75  const label startPoint = faceEdgei;
76 
77  // Pick up starting point and decide direction of circulators
78  const label fwd = f.fcIndex(startPoint);
79  const label back = f.rcIndex(startPoint);
80  const label lEdgeP = f[fwd] == f[f.fcIndex(faceEdgei)] ? fwd : back;
81 
82  // Starting from this point on the face
83  fp = pointsAlpha[f[startPoint]] < pointsAlpha[f[lEdgeP]]
84  ? startPoint : lEdgeP;
85 
86  if
87  (
88  !(
89  (
90  lEdgeP == fwd
91  && pointsAlpha[f[startPoint]] < pointsAlpha[f[lEdgeP]]
92  )
93  || (
94  lEdgeP == back
95  && pointsAlpha[f[startPoint]] > pointsAlpha[f[lEdgeP]]
96  )
97  )
98  )
99  {
100  flipped_ = true;
101  }
102  }
103 
104  label nextFp, edgei;
105  if (flipped_)
106  {
107  nextFp = f.rcIndex(fp);
108  edgei = nextFp;
109  }
110  else
111  {
112  nextFp = f.fcIndex(fp);
113  edgei = fp;
114  }
115 
116  forAll(f, i)
117  {
118  // Edge points field value
119  const scalar& fl = pointsAlpha[f[fp]];
120  const scalar& fk = pointsAlpha[f[nextFp]];
121 
122  const point& pL = points[f[fp]];
123  const point& pK = points[f[nextFp]];
124 
125  const label edg = faceEdges[edgei];
126 
127  // Collect sub-point if the iso-value is bigger than target value
128  if (fl >= target && cutPoints_.size() > 0)
129  {
130  subPoints_.append(pL);
131  if (!unweighted_)
132  {
133  subPointsU_.append(pointsU[f[fp]]);
134  }
135  }
136 
137  // Cut the edge
138  if
139  (
140  !isEdgeCutOld[edg]
141  && ((fl < target && fk >= target) || (fl >= target && fk < target))
142  )
143  {
144  const scalar coeff = (target - fk)/(fl - fk);
145  const point cut = pK + coeff*(pL - pK);
146  if (!unweighted_)
147  {
148  subPointsU_.append
149  (
150  pointsU[f[nextFp]]
151  + coeff*(pointsU[f[fp]] - pointsU[f[nextFp]])
152  );
153  }
154  cutPoints_.append(cut);
155  subPoints_.append(cut);
156  cutEdges_.append(edg);
157  isEdgeCut[edg] = true;
158  }
159 
160  // Termination control
161  if (cutPoints_.size() == 2)
162  {
163  faceEdgei = edgei;
164  return 1;
165  }
166 
167  if (flipped_)
168  {
169  fp = f.rcIndex(fp);
170  nextFp = f.rcIndex(fp);
171  edgei = nextFp;
172  }
173  else
174  {
175  fp = f.fcIndex(fp);
176  nextFp = f.fcIndex(fp);
177  edgei = fp;
178  }
179  }
180 
181  return 0;
182 }
183 
184 
186 (
187  const UList<label>& f,
188  const UList<point>& points,
189  const UList<scalar>& pointsAlpha,
190  const UList<vector>& pointsU,
191  const scalar target,
192  const bool ow
193 )
194 {
195  // Clear all the storage
196  cutPoints_.clear();
197  subPoints_.clear();
198  subPointsU_.clear();
199 
200  // Direction of face circulators
201  flipped_ = ow ? false : true;
202 
203  const UIndirectList<scalar> fAlphas(pointsAlpha, f);
204 
205  // Starting from point with minimum alpha value
206  label fp = findMin(fAlphas);
207 
208  // Second point index
209  label nextFp = flipped_ ? f.rcIndex(fp) : f.fcIndex(fp);
210 
211  forAll(f, i)
212  {
213  // Edge points field value
214  const scalar& fl = pointsAlpha[f[fp]];
215  const scalar& fk = pointsAlpha[f[nextFp]];
216 
217  const point& pL = points[f[fp]];
218  const point& pK = points[f[nextFp]];
219 
220  // Collect sub-point if the iso-value is bigger than target value
221  if (fl >= target && cutPoints_.size() > 0)
222  {
223  subPoints_.append(pL);
224  if (!unweighted_)
225  {
226  subPointsU_.append(pointsU[f[fp]]);
227  }
228  }
229 
230  // Cut the edge
231  if ((fl < target && fk >= target) || (fl >= target && fk < target))
232  {
233  const scalar coeff = (target - fk)/(fl - fk);
234  const point cut = pK + coeff*(pL - pK);
235  if (!unweighted_)
236  {
237  subPointsU_.append
238  (
239  pointsU[f[nextFp]]
240  + coeff*(pointsU[f[fp]] - pointsU[f[nextFp]])
241  );
242  }
243  cutPoints_.append(cut);
244  subPoints_.append(cut);
245  }
246 
247  if (flipped_)
248  {
249  fp = f.rcIndex(fp);
250  nextFp = f.rcIndex(fp);
251  }
252  else
253  {
254  fp = f.fcIndex(fp);
255  nextFp = f.fcIndex(fp);
256  }
257  }
258 
259  // Simple cut
260  if (cutPoints_.size() == 2)
261  {
262  return 1;
263  }
264 
265  // Multiple cuts through the face
266  else if (cutPoints_.size() > 2)
267  {
268  return -1;
269  }
270 
271  // No cut
272  else
273  {
274  return 0;
275  }
276 }
277 
278 
279 // ************************************************************************* //
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
label rcIndex(const label i) const
Return the reverse circular index, i.e. the previous index.
Definition: UListI.H:65
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
label findMin(const ListType &, const label start=0)
Find index of min element (and less than given element).
MPLICface(const bool unweighted=true)
Construct empty.
Definition: MPLICface.C:31
label cutFace(const labelList &f, const labelList &faceEdges, const pointField &points, const boolList &isEdgeCutOld, boolList &isEdgeCut, label &faceEdgei, const UList< scalar > &pointsAlpha, const UList< vector > &pointsU, const label facei, const scalar target, const bool ow)
Function to cut for multi cut.
Definition: MPLICface.C:40