-
Notifications
You must be signed in to change notification settings - Fork 60
/
packingLimiter.H
57 lines (48 loc) · 2.19 KB
/
packingLimiter.H
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
/*---------------------------------------------------------------------------*\
Copyright (C) 2015 Cyrille Bonamy, Julien Chauchat, Tian-Jian Hsu
and contributors
License
This file is part of SedFOAM.
SedFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
SedFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with SedFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
if (packingLimiter)
{
//since the particle pressure near the alphaMax is extremely large
//we need to limit the particle concentration to alphaMax-smallnumber
// Calculating exceeding volume fractions
volScalarField alphaEx(max(alpha - alphaMax+scalar(0.005), scalar(0)));
// Finding neighbouring cells of the whole domain
labelListList neighbour = mesh.cellCells();
scalarField cellVolumes(mesh.cellVolumes());
forAll(alphaEx, celli)
{
// Finding the labels of the neighbouring cells
labelList neighbourCell = neighbour[celli];
// Initializing neighbouring cells contribution
scalar neighboursEx = 0.0;
forAll(neighbourCell, cellj)
{
labelList neighboursNeighbour = neighbour[neighbourCell[cellj]];
scalar neighboursNeighbourCellVolumes = 0.0;
forAll(neighboursNeighbour, cellk)
{
neighboursNeighbourCellVolumes +=
cellVolumes[neighboursNeighbour[cellk]];
}
neighboursEx +=
alphaEx[neighbourCell[cellj]]*cellVolumes[celli]
/neighboursNeighbourCellVolumes;
}
alpha[celli] += neighboursEx - alphaEx[celli];
}
alpha.correctBoundaryConditions();
}