surfaceIteratorPLIC.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) 2019-2020 DLR
9 -------------------------------------------------------------------------------
10 License
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 "surfaceIteratorPLIC.H"
29 #include "scalarMatrices.H"
30 
31 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
32 
34 (
35  const fvMesh& mesh,
36  const scalar tol
37 )
38 :
39  mesh_(mesh),
40  cutCell_(mesh_),
41  surfCellTol_(tol)
42 {}
43 
44 
45 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
46 
48 (
49  const label celli,
50  const scalar alpha1,
51  const scalar tol,
52  const label maxIter,
53  vector normal
54 )
55 {
56  if (mag(normal) == 0)
57  {
59  << "normal length is zero in cell: " << celli << nl
60  << "try increasing nCorrectors" << endl;
61 
62  return sign(alpha1-0.5);
63  }
64 
65  normal.normalise();
66 
67  // Finding cell vertex extrema values
68  const labelList& pLabels = mesh_.cellPoints(celli);
69  scalarField fvert(pLabels.size());
70  forAll(pLabels, pi)
71  {
72  fvert[pi] = (mesh_.points()[pLabels[pi]] - mesh_.C()[celli]) & (normal);
73  }
74 
75  const labelList order(Foam::sortedOrder(fvert));
76  scalar f1 = fvert[order.first()];
77  scalar f2 = fvert[order.last()];
78 
79 
80  //Handling special case where method is handed an almost full or empty cell
81  if (alpha1 < tol)
82  {
83  return -1;
84  }
85  else if (1 - alpha1 < tol)
86  {
87  return 1;
88  }
89 
90  // Finding the two vertices inbetween which the isovalue giving alpha1 lies
91  label L1 = 0;
92  label L2 = fvert.size() - 1;
93  scalar a1 = 1;
94  scalar a2 = 0;
95  scalar L3, f3, a3;
96 
97  while (L2 - L1 > 1)
98  {
99  L3 = round(0.5*(L1 + L2));
100  f3 = fvert[order[L3]];
101  cutCell_.calcSubCell(celli, f3, normal);
102  a3 = cutCell_.VolumeOfFluid();
103  if (a3 > alpha1)
104  {
105  L1 = L3; f1 = f3; a1 = a3;
106  }
107  else if (a3 < alpha1)
108  {
109  L2 = L3; f2 = f3; a2 = a3;
110  }
111  else
112  {
113  return cutCell_.calcSubCell(celli, f3, normal);
114  }
115  }
116 
117  if (mag(f1 - f2) < 10*SMALL)
118  {
119  return cutCell_.calcSubCell(celli, f1, normal);
120  }
121 
122  if (mag(a1 - a2) < tol)
123  {
124  return cutCell_.calcSubCell(celli, 0.5*(f1 + f2), normal);
125  }
126  // Now we know that a(f) = alpha1 is to be found on the f interval
127  // [f1, f2], i.e. alpha1 will be in the interval [a2,a1]
128 
129 
130  // Finding coefficients in 3 deg polynomial alpha(f) from 4 solutions
131 
132  // Finding 2 additional points on 3 deg polynomial
133  f3 = f1 + (f2 - f1)/3;
134  cutCell_.calcSubCell(celli, f3, normal);
135  a3 = cutCell_.VolumeOfFluid();
136 
137  scalar f4 = f1 + 2*(f2 - f1)/3;
138  cutCell_.calcSubCell(celli, f4, normal);
139  scalar a4 = cutCell_.VolumeOfFluid();
140 
141  // Building and solving Vandermonde matrix equation
142  scalarField a(4), f(4), C(4), fOld(4);
143  {
144  a[0] = a1, a[1] = a3, a[2] = a4, a[3] = a2;
145  fOld[0] = f1, fOld[1] = f3, fOld[2] = f4, fOld[3] = f2;
146  f[0] = 0, f[1] = (f3-f1)/(f2-f1), f[2] = (f4-f1)/(f2-f1), f[3] = 1;
148  forAll(f, i)
149  {
150  forAll(f, j)
151  {
152  M[i][j] = pow(f[i], 3 - j);
153  }
154  }
155  // C holds the 4 polynomial coefficients
156  C = a;
157  LUsolve(M, C);
158  }
159 
160  // Finding root with Newton method
161 
162  f3 = f[1]; a3 = a[1];
163  label nIter = 0;
164  scalar res = mag(a3 - alpha1);
165  while (res > tol && nIter < 10*maxIter)
166  {
167  f3 -= (C[0]*pow3(f3) + C[1]*sqr(f3) + C[2]*f3 + C[3] - alpha1)
168  /(3*C[0]*sqr(f3) + 2*C[1]*f3 + C[2]);
169  a3 = C[0]*pow3(f3) + C[1]*sqr(f3) + C[2]*f3 + C[3];
170  res = mag(a3 - alpha1);
171  nIter++;
172  }
173  // Scaling back to original range
174  f3 = f3*(f2 - f1) + f1;
175 
176  //Check result
177  label status = cutCell_.calcSubCell(celli, f3, normal);
178  const scalar VOF = cutCell_.VolumeOfFluid();
179  res = mag(VOF - alpha1);
180 
181  if (res > tol)
182  {
183  }
184  else
185  {
186  return status;
187  }
188 
189  // If tolerance not met use the secant method with f3 as a hopefully very
190  // good initial guess to crank res the last piece down below tol
191 
192  scalar x2 = f3;
193  scalar g2 = VOF - alpha1;
194  scalar x1 = max(1e-3*(f2 - f1), 100*SMALL);
195  x1 = max(x1, f1);
196  x1 = min(x1, f2);
197  cutCell_.calcSubCell(celli, x1,normal);
198  scalar g1 = cutCell_.VolumeOfFluid() - alpha1;
199 
200  nIter = 0;
201  scalar g0(0), x0(0);
202  while (res > tol && nIter < maxIter && g1 != g2)
203  {
204  x0 = (x2*g1 - x1*g2)/(g1 - g2);
205  status = cutCell_.calcSubCell(celli, x0, normal);
206  g0 = cutCell_.VolumeOfFluid() - alpha1;
207  res = mag(g0);
208  x2 = x1; g2 = g1;
209  x1 = x0; g1 = g0;
210  nIter++;
211  }
212 
213  return status;
214 }
215 
216 
217 // ************************************************************************* //
Foam::LUsolve
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
Definition: scalarMatricesTemplates.C:215
surfaceIteratorPLIC.H
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
alpha1
const volScalarField & alpha1
Definition: setRegionFluidFields.H:8
Foam::Vector::normalise
Vector< Cmpt > & normalise()
Normalise the vector by its magnitude.
Definition: VectorI.H:123
Foam::surfaceIteratorPLIC::vofCutCell
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
Definition: surfaceIteratorPLIC.C:48
Foam::sign
dimensionedScalar sign(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:166
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
C
volScalarField & C
Definition: readThermalProperties.H:102
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::Field< scalar >
Foam::pow3
dimensionedScalar pow3(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:89
Foam::pow
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Definition: dimensionedScalar.C:75
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::fvMesh
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:85
Foam::SquareMatrix< scalar >
Foam::sqr
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Definition: dimensionedSymmTensor.C:51
Foam::constant::mathematical::pi
constexpr scalar pi(M_PI)
Foam::nl
constexpr char nl
Definition: Ostream.H:404
f
labelList f(nPoints)
Foam::Vector< scalar >
Foam::List< label >
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
M
#define M(I)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::surfaceIteratorPLIC::surfaceIteratorPLIC
surfaceIteratorPLIC(const fvMesh &mesh, const scalar tol)
Construct from fvMesh and a scalarField.
Definition: surfaceIteratorPLIC.C:34
scalarMatrices.H
Foam::sortedOrder
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
Foam::C
Graphite solid properties.
Definition: C.H:50
WarningInFunction
#define WarningInFunction
Report a warning using Foam::Warning.
Definition: messageStream.H:328