setAlphaField.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) 2016-2017 DHI
9 Copyright (C) 2017-2020 OpenCFD Ltd.
10 Copyright (C) 2017-2020 German Aerospace Center (DLR)
11 Copyright (C) 2020 Johan Roenby
12-------------------------------------------------------------------------------
13License
14 This file is part of OpenFOAM.
15
16 OpenFOAM is free software: you can redistribute it and/or modify it
17 under the terms of the GNU General Public License as published by
18 the Free Software Foundation, either version 3 of the License, or
19 (at your option) any later version.
20
21 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
22 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
23 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
24 for more details.
25
26 You should have received a copy of the GNU General Public License
27 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
28
29Application
30 setAlphaField
31
32Description
33 Uses cutCellIso to create a volume fraction field from either a cylinder,
34 a sphere or a plane.
35
36 Original code supplied by Johan Roenby, DHI (2016)
37 Modification Henning Scheufler, DLR (2019)
38
39\*---------------------------------------------------------------------------*/
40
41#include "fvCFD.H"
42
43#include "triSurface.H"
44#include "triSurfaceTools.H"
45
46#include "implicitFunction.H"
47
48#include "cutCellIso.H"
49#include "cutFaceIso.H"
50#include "OBJstream.H"
51
52
53// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
54
55void isoFacesToFile
56(
57 const DynamicList<List<point>>& faces,
58 const word& fileName
59)
60{
61 // Writing isofaces to OBJ file for inspection in paraview
62 OBJstream os(fileName + ".obj");
63
64 if (Pstream::parRun())
65 {
66 // Collect points from all the processors
67 List<DynamicList<List<point>>> allProcFaces(Pstream::nProcs());
68 allProcFaces[Pstream::myProcNo()] = faces;
69 Pstream::gatherList(allProcFaces);
70
71 if (Pstream::master())
72 {
73 Info<< "Writing file: " << fileName << endl;
74
75 for (const DynamicList<List<point>>& procFaces : allProcFaces)
76 {
77 for (const List<point>& facePts : procFaces)
78 {
79 os.write(face(identity(facePts.size())), facePts);
80 }
81 }
82 }
83 }
84 else
85 {
86 Info<< "Writing file: " << fileName << endl;
87
88 for (const List<point>& facePts : faces)
89 {
90 os.write(face(identity(facePts.size())), facePts);
91 }
92 }
93}
94
95void setAlpha
96(
97 volScalarField& alpha1,
98 DynamicList<List<point>>& facePts,
99 scalarField& f,
100 const bool writeOBJ
101)
102{
103 const fvMesh& mesh = alpha1.mesh();
104 cutCellIso cutCell(mesh, f);
105 cutFaceIso cutFace(mesh, f);
106
107 forAll(alpha1, cellI)
108 {
109 cutCell.calcSubCell(cellI, 0.0);
110
111 alpha1[cellI] = max(min(cutCell.VolumeOfFluid(), 1), 0);
112
113 if (writeOBJ && (mag(cutCell.faceArea()) >= 1e-14))
114 {
115 facePts.append(cutCell.facePoints());
116 }
117 }
118
119 // Setting boundary alpha1 values
120 forAll(mesh.boundary(), patchi)
121 {
122 if (mesh.boundary()[patchi].size() > 0)
123 {
124 const label start = mesh.boundary()[patchi].patch().start();
125 scalarField& alphap = alpha1.boundaryFieldRef()[patchi];
126 const scalarField& magSfp = mesh.magSf().boundaryField()[patchi];
127
128 forAll(alphap, patchFacei)
129 {
130 const label facei = patchFacei + start;
131 cutFace.calcSubFace(facei, 0.0);
132 alphap[patchFacei] =
133 mag(cutFace.subFaceArea())/magSfp[patchFacei];
134 }
135 }
136 }
137}
138
139
140int main(int argc, char *argv[])
141{
142 argList::addNote
143 (
144 "Uses cutCellIso to create a volume fraction field from an "
145 "implicit function."
146 );
147
148 argList::addOption
149 (
150 "dict",
151 "file",
152 "Alternative setAlphaFieldDict dictionary"
153 );
154 #include "addRegionOption.H"
155 #include "setRootCase.H"
156 #include "createTime.H"
157 #include "createNamedMesh.H"
158
159 const word dictName("setAlphaFieldDict");
161
162 IOdictionary setAlphaFieldDict(dictIO);
163
164 Info<< "Reading " << setAlphaFieldDict.name() << endl;
165
166 const word fieldName = setAlphaFieldDict.get<word>("field");
167 const bool invert = setAlphaFieldDict.getOrDefault("invert", false);
168 const bool writeOBJ = setAlphaFieldDict.getOrDefault("writeOBJ", false);
169
170 Info<< "Reading field " << fieldName << nl << endl;
172 (
173 IOobject
174 (
175 fieldName,
176 runTime.timeName(),
177 mesh,
178 IOobject::MUST_READ,
179 IOobject::AUTO_WRITE
180 ),
181 mesh
182 );
183
184 autoPtr<implicitFunction> func = implicitFunction::New
185 (
186 setAlphaFieldDict.get<word>("type"),
187 setAlphaFieldDict
188 );
189
190 scalarField f(mesh.nPoints(), Zero);
191
192 forAll(f, pi)
193 {
194 f[pi] = func->value(mesh.points()[pi]);
195 };
196
197 DynamicList<List<point>> facePts;
198
199 setAlpha(alpha1, facePts, f, writeOBJ);
200
201 if (writeOBJ)
202 {
203 isoFacesToFile(facePts, fieldName + "0");
204 }
205
206 ISstream::defaultPrecision(18);
207
208 if (invert)
209 {
210 alpha1 = scalar(1) - alpha1;
211 }
212
213 Info<< "Writing new alpha field " << alpha1.name() << endl;
214 alpha1.write();
215
216 const scalarField& alpha = alpha1.internalField();
217
218 Info<< "sum(alpha*V):" << gSum(mesh.V()*alpha)
219 << ", 1-max(alpha1): " << 1 - gMax(alpha)
220 << " min(alpha1): " << gMin(alpha) << endl;
221
222 Info<< "\nEnd\n" << endl;
223
224 return 0;
225}
226
227
228// ************************************************************************* //
Y[inertIndex] max(0.0)
const volScalarField & alpha1
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
OBJstream os(runTime.globalPath()/outputName)
const word dictName("faMeshDefinition")
constexpr scalar pi(M_PI)
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of a point.
Definition: meshTools.C:203
Type gSum(const FieldField< Field, Type > &f)
labelList identity(const label len, label start=0)
Return an identity map of the given length with (map[i] == i)
Definition: labelList.C:38
void func(FieldField< Field, Type > &f, const FieldField< Field, Type > &f1, const FieldField< Field, Type > &f2)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
labelList invert(const label len, const labelUList &map)
Create an inverse one-to-one mapping.
Definition: ListOps.C:36
Type gMin(const FieldField< Field, Type > &f)
Type gMax(const FieldField< Field, Type > &f)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
volScalarField & alpha
IOobject dictIO
volScalarField & e
Definition: createFields.H:11
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333