rotateMesh performs faulty operations in some boundary conditions
Summary
When applied to cases containing specific boundary conditions defined for vectorial or tensorial fields, the utility rotateMesh performs a double rotation, due to inadequate behaviour when writing transformed points into the polyMesh directory. If multiple time directories are present, only the first one is affected.
Steps to reproduce
It is possible to reproduce the bug by setting an OpenFOAM case including specific boundary conditions for vectorial (or tensorial) fields that need to read geometric data from the mesh. For instance, in the simpleFoam tutorial pitzDaily, if we replace the inlet boundaryField condition with a flowRateInletVelocity condition equivalent to the original fixedValue condition, in the 0/U file, like:
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: v2106 |
| \\ / A nd | Website: www.openfoam.com |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class volVectorField;
object U;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
/*inlet
{
type fixedValue;
value uniform (10 0 0);
}*/
inlet
{
type flowRateInletVelocity;
volumetricFlowRate 2.54e-4;
}
outlet
{
type zeroGradient;
}
upperWall
{
type noSlip;
}
lowerWall
{
type noSlip;
}
frontAndBack
{
type empty;
}
}
// ************************************************************************* //And then perform a simple +45 degree rotation along the z-axis using the utility rotateMesh (could be any arbitrary rotation), the velocity vector defined in the inlet boundaryField for U is erroneously rotated twice, completing a 90 degree rotation.
This +45 degree rotation along the z-axis can be achieved by this Allrun shell script:
#!/bin/sh
cd "${0%/*}" || exit # Run from this directory
. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions # Tutorial run functions
#------------------------------------------------------------------------------
# Create 0.orig directory for backup
[ -d 0 ] && mv 0 0.orig
cp -r 0.orig/ 0/
runApplication blockMesh
runApplication rotateMesh "(1 0 0)" "(0.707107 0.707107 0)"
runApplication $(getApplication)
#------------------------------------------------------------------------------Example case
Find attached the example case described above.
What is the current bug behaviour?
The utility rotateMesh performs a double rotation to the velocity vector defined in the inlet boundaryField for U.
This is due to an inadequate behaviour in the rotateMesh algorithm, which writes the transformed points into the polyMesh directory before the vector and tensor fields are rotated. Vector or tensor boundary conditions that must read from the mesh points in order to write corresponding vector/tensor values are hence rotated twice.
What is the expected correct behavior?
The expected correct behaviour can be obtained by using the utility transformPoints instead, for reproducing the same rotation, with the rotateFields flag. The commands below can be employed and lead to the same +45 degree rotation, however without breaking the inlet boundaryField velocity vector:
$ transformPoints -rotate '((1 0 0) (0.707107 0.707107 0))' -rotateFields
or
$ transformPoints -rotate-angle '((0 0 1) 45)' -rotateFields
Relevant logs and/or images
Erroneous inlet boundaryField for the velocity after the rotation (produced by the rotateMesh utility):
inlet
{
type flowRateInletVelocity;
volumetricFlowRate constant 0.000254;
value nonuniform List<vector>
30
(
(3.730349363e-14 10 0)
(-2.069228318e-07 10 0)
(3.375077995e-14 10 0)
(1.80136758e-07 10 -8.429032341e-17)
(0 10 0)
(-1.568181549e-07 10 -7.337898766e-17)
(1.463164896e-07 10 -6.846501721e-17)
(-1.365181355e-07 10 0)
(1.273759462e-07 10 -1.192045185e-16)
(-1.188459393e-07 10 1.112217487e-16)
(1.108871617e-07 10 0)
(-1.231714002e-07 10 0)
(-3.463895837e-14 10 -1.891196916e-16)
(8.288824915e-08 10 0)
(-4.759729677e-08 10 0)
(-5.577976303e-09 10 2.088050669e-16)
(3.660640857e-08 10 0)
(-7.507386712e-09 10 -2.810307806e-16)
(-1.539641836e-08 10 0)
(1.261886595e-08 10 0)
(-1.449527609e-08 10 0)
(1.554063989e-08 10 2.077669184e-16)
(3.18776916e-09 10 -2.386615189e-16)
(1.318243825e-08 10 -2.741500913e-16)
(-1.766639723e-08 10 0)
(1.932702354e-08 10 0)
(-3.33014194e-08 10 0)
(6.375543737e-08 10 -1.193307603e-16)
(-5.858861751e-08 10 0)
(6.73006717e-08 10 0)
)
;Correct inlet boundaryField for the velocity after the rotation (produced by the transformPoints utility):
inlet
{
type flowRateInletVelocity;
volumetricFlowRate constant 0.000254;
value uniform (7.07107 7.07107 0);
}Environment information
- OpenFOAM version : v2106
- Operating system : Ubuntu 20.04.2 LTS
- Hardware info : Intel(R) Core(TM) i7-10510U CPU @ 1.80GHz
- Compiler : gcc 9.3.0
Possible fixes
This issue can be fixed by constructing the pointer-list of fvMesh objects inside the rotateMesh utility before the intermediate operation of points transformation is performed.
This can be achieved by inserting the include file createNamedMeshes.H before the first program loop. The relevant section of the code would read:
// ------------------------------------------------------------------------
#include "createNamedMeshes.H"
forAll(regionNames, regioni)
{
const word& regionName = regionNames[regioni];
const word& regionDir =
(
regionName == polyMesh::defaultRegion ? word::null : regionName
);
const fileName meshDir = regionDir/polyMesh::meshSubDir;
if (regionNames.size() > 1)
{
Info<< "region=" << regionName << nl;
}
pointIOField points
(
IOobject
(
"points",
runTime.findInstance(meshDir, "points"),
meshDir,
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
points = transform(rotT, points);
// Set the precision of the points data to 10
IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
Info<< "Writing points into directory "
<< runTime.relativePath(points.path()) << nl
<< endl;
points.write();
}
instantList timeDirs = timeSelector::select0(runTime, args);
forAll(timeDirs, timei)
{
// ...
}Find attached the rotateMesh.C C++ file with the proposed code, that patches this specific bug.
PS: I personally consider this solution to be naive and somewhat simplistic, where a more elegant and robust one would reorganise the rotateMesh algorithm to match the order of operations in the transformPoints utility (points transformation -> mesh objects construction -> fields rotation -> mesh writing). However, I could not find an example case in which the simple solution suggested above fails.