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-------------------------------------------------------------------------------
10License
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// ************************************************************************* //
#define M(I)
const volScalarField & alpha1
Graphite solid properties.
Definition: C.H:53
T & first()
Return the first element of the list.
Definition: UListI.H:202
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
Vector< Cmpt > & normalise(const scalar tol=ROOTVSMALL)
Inplace normalise the vector by its magnitude.
Definition: VectorI.H:123
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:91
Finds the cutValue that matches the volume fraction.
label vofCutCell(const label celli, const scalar alpha1, const scalar tol, const label maxIter, vector normal)
Finds matching cutValue for the given value fraction.
dynamicFvMesh & mesh
#define WarningInFunction
Report a warning using Foam::Warning.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar sign(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedScalar pow3(const dimensionedScalar &ds)
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
void LUsolve(scalarSquareMatrix &matrix, List< Type > &source)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333