From 54dfcf5046f070081595f7df8707ef94d90155d2 Mon Sep 17 00:00:00 2001
From: Johan Roenby <johan.roenby@gmail.com>
Date: Mon, 21 Dec 2020 09:20:32 +0000
Subject: [PATCH] BUG: setAlphaField: fix incompatibilities with BCs (#1962)
Co-authored-by: Johan Roenby <johan.roenby@gmail.com>
Co-authored-by: Henning Scheufler <Henning.Scheufler@dlr.de>
---
.../setAlphaField/setAlphaField.C | 81 +++++++++++--------
1 file changed, 48 insertions(+), 33 deletions(-)
diff --git a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
index 82e37fac648..c340b3f224c 100644
--- a/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
+++ b/applications/utilities/preProcessing/setAlphaField/setAlphaField.C
@@ -7,7 +7,8 @@
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 DHI
Copyright (C) 2017-2020 OpenCFD Ltd.
- Copyright (c) 2017-2020, German Aerospace Center (DLR)
+ Copyright (C) 2017-2020 German Aerospace Center (DLR)
+ Copyright (C) 2020 Johan Roenby
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -45,6 +46,7 @@ Description
#include "implicitFunction.H"
#include "cutCellIso.H"
+#include "cutFaceIso.H"
#include "OBJstream.H"
@@ -90,6 +92,50 @@ void isoFacesToFile
}
}
+void setAlpha
+(
+ volScalarField& alpha1,
+ DynamicList<List<point>>& facePts,
+ scalarField& f,
+ const bool writeOBJ
+)
+{
+ const fvMesh& mesh = alpha1.mesh();
+ cutCellIso cutCell(mesh, f);
+ cutFaceIso cutFace(mesh, f);
+
+ forAll(alpha1, cellI)
+ {
+ cutCell.calcSubCell(cellI, 0.0);
+
+ alpha1[cellI] = max(min(cutCell.VolumeOfFluid(), 1), 0);
+
+ if (writeOBJ && (mag(cutCell.faceArea()) >= 1e-14))
+ {
+ facePts.append(cutCell.facePoints());
+ }
+ }
+
+ // Setting boundary alpha1 values
+ forAll(mesh.boundary(), patchi)
+ {
+ if (mesh.boundary()[patchi].size() > 0)
+ {
+ const label start = mesh.boundary()[patchi].patch().start();
+ scalarField& alphap = alpha1.boundaryFieldRef()[patchi];
+ const scalarField& magSfp = mesh.magSf().boundaryField()[patchi];
+
+ forAll(alphap, patchFacei)
+ {
+ const label facei = patchFacei + start;
+ cutFace.calcSubFace(facei, 0.0);
+ alphap[patchFacei] =
+ mag(cutFace.subFaceArea())/magSfp[patchFacei];
+ }
+ }
+ }
+}
+
int main(int argc, char *argv[])
{
@@ -143,38 +189,9 @@ int main(int argc, char *argv[])
f[pi] = func->value(mesh.points()[pi]);
};
- cutCellIso cutCell(mesh, f);
-
DynamicList<List<point>> facePts;
- DynamicList<triSurface> surface;
-
- surfaceScalarField cellToCellDist(mag(mesh.delta()));
-
- forAll(alpha1, cellI)
- {
- label cellStatus = cutCell.calcSubCell(cellI, 0.0);
-
- if (cellStatus == -1)
- {
- alpha1[cellI] = 1;
- }
- else if (cellStatus == 1)
- {
- alpha1[cellI] = 0;
- }
- else if (cellStatus == 0)
- {
- if (mag(cutCell.faceArea()) != 0)
- {
- alpha1[cellI] = max(min(cutCell.VolumeOfFluid(), 1), 0);
- if (writeOBJ && (mag(cutCell.faceArea()) >= 1e-14))
- {
- facePts.append(cutCell.facePoints());
- }
- }
- }
- }
+ setAlpha(alpha1, facePts, f, writeOBJ);
if (writeOBJ)
{
@@ -188,8 +205,6 @@ int main(int argc, char *argv[])
alpha1 = scalar(1) - alpha1;
}
- alpha1.correctBoundaryConditions();
-
Info<< "Writing new alpha field " << alpha1.name() << endl;
alpha1.write();
--
GitLab