diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index 5d58e43fbab2bcc51bd13d6a0bd007b90362780e..98d15c770fc7ae23397b0cc362fe31a9eee0d0fa 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -1903,6 +1903,10 @@ __global__ void interpolateStaticToRotating(
real dx);
__global__ void interpolateRotatingToStatic(
+ real *distributionsStatic,
+ real *distributionsRotating,
+ unsigned int numberOfLBNodesStatic,
+ unsigned int numberOfLBNodesRotating,
unsigned int numberOfInterfaceNodes,
const unsigned int *indicesStatic,
unsigned int *indicesRotatingCell,
@@ -1925,6 +1929,8 @@ __global__ void interpolateRotatingToStatic(
real angularVelocityX,
real angularVelocityY,
real angularVelocityZ,
+ real omegaRotating,
+ bool isEvenTimestep,
real dx);
__global__ void updateGlobalCoordinates(
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
new file mode 100644
index 0000000000000000000000000000000000000000..c4e66a9ee3df0ad44eeda31f27b16a78beddb32d
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateRotatingToStatic.cu
@@ -0,0 +1,321 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids 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.
+//
+// VirtualFluids 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 VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file scaleCF_compressible.cu
+//! \ingroup GPU/GridScaling
+//! \author Martin Schoenherr, Anna Wellmann
+//=======================================================================================
+
+#include "DataTypes.h"
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+#include "LBM/GPUHelperFunctions/ChimeraTransformation.h"
+#include "LBM/GPUHelperFunctions/ScalingUtilities.h"
+#include "LBM/GPUHelperFunctions/GridTraversion.h"
+#include "LBM/GPUHelperFunctions/CoordinateTransformation.h"
+
+#include <lbm/MacroscopicQuantities.h>
+#include <lbm/refinement/Coefficients.h>
+#include <lbm/refinement/InterpolationCF.h>
+
+using namespace vf::lbm;
+using namespace vf::gpu;
+
+__device__ __inline__ void calculateRotationalForces(real &forceX, real &forceY, real &forceZ, real angularVelocityX,
+ real angularVelocityY, real angularVelocityZ, real coordSourceXlocal,
+ real coordSourceYlocal, real coordSourceZlocal, real vx1, real vx2,
+ real vx3)
+{
+ forceX = -((-(angularVelocityY * angularVelocityY * coordSourceXlocal) -
+ angularVelocityZ * angularVelocityZ * coordSourceXlocal +
+ angularVelocityX * angularVelocityY * coordSourceYlocal -
+ 2 * angularVelocityX * angularVelocityX * angularVelocityZ * coordSourceYlocal -
+ 2 * angularVelocityY * angularVelocityY * angularVelocityZ * coordSourceYlocal -
+ 2 * angularVelocityZ * angularVelocityZ * angularVelocityZ * coordSourceYlocal +
+ 2 * angularVelocityX * angularVelocityX * angularVelocityY * coordSourceZlocal +
+ 2 * angularVelocityY * angularVelocityY * angularVelocityY * coordSourceZlocal +
+ angularVelocityX * angularVelocityZ * coordSourceZlocal +
+ 2 * angularVelocityY * angularVelocityZ * angularVelocityZ * coordSourceZlocal +
+ 4 * angularVelocityY * angularVelocityY * vx1 + 4 * angularVelocityZ * angularVelocityZ * vx1 -
+ 4 * angularVelocityX * angularVelocityY * vx2 - 2 * angularVelocityZ * vx2 + 2 * angularVelocityY * vx3 -
+ 4 * angularVelocityX * angularVelocityZ * vx3) /
+ (1 + 4 * angularVelocityX * angularVelocityX + 4 * angularVelocityY * angularVelocityY +
+ 4 * angularVelocityZ * angularVelocityZ));
+
+ forceY = -((angularVelocityX * angularVelocityY * coordSourceXlocal +
+ 2 * angularVelocityX * angularVelocityX * angularVelocityZ * coordSourceXlocal +
+ 2 * angularVelocityY * angularVelocityY * angularVelocityZ * coordSourceXlocal +
+ 2 * angularVelocityZ * angularVelocityZ * angularVelocityZ * coordSourceXlocal -
+ angularVelocityX * angularVelocityX * coordSourceYlocal -
+ angularVelocityZ * angularVelocityZ * coordSourceYlocal -
+ 2 * angularVelocityX * angularVelocityX * angularVelocityX * coordSourceZlocal -
+ 2 * angularVelocityX * angularVelocityY * angularVelocityY * coordSourceZlocal +
+ angularVelocityY * angularVelocityZ * coordSourceZlocal -
+ 2 * angularVelocityX * angularVelocityZ * angularVelocityZ * coordSourceZlocal -
+ 4 * angularVelocityX * angularVelocityY * vx1 + 2 * angularVelocityZ * vx1 +
+ 4 * angularVelocityX * angularVelocityX * vx2 + 4 * angularVelocityZ * angularVelocityZ * vx2 -
+ 2 * angularVelocityX * vx3 - 4 * angularVelocityY * angularVelocityZ * vx3) /
+ (1 + 4 * angularVelocityX * angularVelocityX + 4 * angularVelocityY * angularVelocityY +
+ 4 * angularVelocityZ * angularVelocityZ));
+
+ forceZ = -((-2 * angularVelocityX * angularVelocityX * angularVelocityY * coordSourceXlocal -
+ 2 * angularVelocityY * angularVelocityY * angularVelocityY * coordSourceXlocal +
+ angularVelocityX * angularVelocityZ * coordSourceXlocal -
+ 2 * angularVelocityY * angularVelocityZ * angularVelocityZ * coordSourceXlocal +
+ 2 * angularVelocityX * angularVelocityX * angularVelocityX * coordSourceYlocal +
+ 2 * angularVelocityX * angularVelocityY * angularVelocityY * coordSourceYlocal +
+ angularVelocityY * angularVelocityZ * coordSourceYlocal +
+ 2 * angularVelocityX * angularVelocityZ * angularVelocityZ * coordSourceYlocal -
+ angularVelocityX * angularVelocityX * coordSourceZlocal -
+ angularVelocityY * angularVelocityY * coordSourceZlocal - 2 * angularVelocityY * vx1 -
+ 4 * angularVelocityX * angularVelocityZ * vx1 + 2 * angularVelocityX * vx2 -
+ 4 * angularVelocityY * angularVelocityZ * vx2 + 4 * angularVelocityX * angularVelocityX * vx3 +
+ 4 * angularVelocityY * angularVelocityY * vx3) /
+ (1 + 4 * angularVelocityX * angularVelocityX + 4 * angularVelocityY * angularVelocityY +
+ 4 * angularVelocityZ * angularVelocityZ));
+}
+
+__device__ __inline__ void calculateHalfRotationalForcesByVertex(real *forcesVertexX, real *forcesVertexY, real *forcesVertexZ,
+ uint index, InterpolationVertex direction,
+ const Distributions27 &distRotating, const real *coordSourceX,
+ const real *coordSourceY, const real *coordSourceZ,
+ real angularVelocityX, real angularVelocityY,
+ real angularVelocityZ)
+{
+ real drho, vx1, vx2, vx3;
+ real coordSourceXlocal, coordSourceYlocal, coordSourceZlocal;
+ real forceX, forceY, forceZ;
+
+ getCompressibleMacroscopicValues(distRotating, index, drho, vx1, vx2, vx3);
+
+ coordSourceXlocal = coordSourceX[index];
+ coordSourceYlocal = coordSourceY[index];
+ coordSourceZlocal = coordSourceZ[index];
+
+ calculateRotationalForces(forceX, forceY, forceZ, angularVelocityX, angularVelocityY, angularVelocityZ,
+ coordSourceXlocal, coordSourceYlocal, coordSourceZlocal, vx1, vx2, vx3);
+ forceX *= c1o2;
+ forceY *= c1o2;
+ forceZ *= c1o2;
+
+ forcesVertexX[direction] = forceX;
+ forcesVertexY[direction] = forceY;
+ forcesVertexZ[direction] = forceZ;
+}
+
+__device__ __inline__ void calculateTangentialVelocities(real *tangentialVelocitiesX, real *tangentialVelocitiesY,
+ real *tangentialVelocitiesZ, uint index,
+ InterpolationVertex direction, const real *coordSourceX,
+ const real *coordSourceY, const real *coordSourceZ,
+ real angularVelocityX, real angularVelocityY, real angularVelocityZ, real dx)
+{
+ real coordSourceXlocal = coordSourceX[index] / dx;
+ real coordSourceYlocal = coordSourceY[index] / dx;
+ real coordSourceZlocal = coordSourceZ[index] / dx;
+ tangentialVelocitiesX[direction] = -(angularVelocityZ * coordSourceYlocal) + angularVelocityY * coordSourceZlocal;
+ tangentialVelocitiesY[direction] = angularVelocityZ * coordSourceXlocal - angularVelocityX * coordSourceZlocal;
+ tangentialVelocitiesZ[direction] = -(angularVelocityY * coordSourceXlocal) + angularVelocityX * coordSourceYlocal;
+}
+
+__global__ void interpolateRotatingToStatic(
+ real *distributionsStatic,
+ real *distributionsRotating,
+ unsigned int numberOfLBNodesStatic,
+ unsigned int numberOfLBNodesRotating,
+ unsigned int numberOfInterfaceNodes,
+ const unsigned int *indicesStatic,
+ unsigned int *indicesRotatingCell,
+ const real *coordDestinationX,
+ const real *coordDestinationY,
+ const real *coordDestinationZ,
+ const real *coordSourceX,
+ const real *coordSourceY,
+ const real *coordSourceZ,
+ const uint *neighborXrotating,
+ const uint *neighborYrotating,
+ const uint *neighborZrotating,
+ const uint *neighborMMMrotating,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ,
+ real angularVelocityX,
+ real angularVelocityY,
+ real angularVelocityZ,
+ real omegaRotating,
+ bool isEvenTimestep,
+ real dx)
+{
+ const unsigned listIndex = vf::gpu::getNodeIndex();
+ if (listIndex >= numberOfInterfaceNodes) return;
+ const uint destinationIndex = indicesStatic[listIndex];
+ const uint previousSourceIndex = indicesRotatingCell[listIndex];
+ const uint indexNeighborMMMsource = neighborMMMrotating[previousSourceIndex];
+
+ real rotatedCoordDestinationX;
+ real rotatedCoordDestinationY;
+ real rotatedCoordDestinationZ;
+ transformGlobalToRotating(rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ,
+ coordDestinationX[destinationIndex], coordDestinationY[destinationIndex],
+ coordDestinationZ[destinationIndex], centerCoordX, centerCoordY, centerCoordZ, angleX, angleY,
+ angleZ);
+
+ const uint sourceIndex = traverseSourceCell(rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ,
+ indexNeighborMMMsource, coordSourceX[indexNeighborMMMsource],
+ coordSourceY[indexNeighborMMMsource], coordSourceZ[indexNeighborMMMsource],
+ neighborXrotating, neighborYrotating, neighborZrotating, dx);
+
+ indicesRotatingCell[listIndex] = sourceIndex;
+
+ Distributions27 distRotating;
+ vf::gpu::getPointersToDistributions(distRotating, distributionsRotating, numberOfLBNodesRotating, isEvenTimestep);
+
+ real forcesVertexX [8];
+ real forcesVertexY [8];
+ real forcesVertexZ [8];
+ real tangentialVelocitiesX [8];
+ real tangentialVelocitiesY [8];
+ real tangentialVelocitiesZ [8];
+
+ // forces and tangential velocites
+ // MMM
+ uint indexTemp = sourceIndex;
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::MMM,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::MMM, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // MMP
+ indexTemp = neighborZrotating[indexTemp];
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::MMP,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::MMP, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // PMP
+ indexTemp = neighborXrotating[indexTemp];
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::PMP,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::PMP, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // PMM
+ indexTemp = neighborXrotating[sourceIndex]; // back to the base node of the cell --> sourceIndex
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::PMM,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::PMM, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // PPM
+ indexTemp = neighborYrotating[indexTemp];
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::PPM,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::PPM, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // MPM
+ indexTemp = neighborYrotating[sourceIndex]; // back to the base node of the cell --> sourceIndex
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::MPM,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::MPM, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // MPP
+ indexTemp = neighborZrotating[indexTemp];
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::MPP,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::MPP, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+ // PPP
+ indexTemp = neighborXrotating[indexTemp];
+ calculateHalfRotationalForcesByVertex(forcesVertexX, forcesVertexY, forcesVertexZ, indexTemp, InterpolationVertex::PPP,
+ distRotating, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ);
+ calculateTangentialVelocities(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ, indexTemp,
+ InterpolationVertex::PPP, coordSourceX, coordSourceY, coordSourceZ, angularVelocityX,
+ angularVelocityY, angularVelocityZ, dx);
+
+ // 1. calculate moments for the nodes of the source cell, add half of the rotational force to the velocity during the computation of the moments
+ vf::lbm::MomentsOnSourceNodeSet momentsSet;
+ vf::gpu::calculateMomentSet<false>(momentsSet, listIndex, distributionsRotating, neighborXrotating, neighborYrotating,
+ neighborZrotating, indicesRotatingCell, nullptr, numberOfLBNodesStatic, omegaRotating,
+ isEvenTimestep, forcesVertexX, forcesVertexY, forcesVertexZ);
+
+ // 2. calculate the coefficients for the interpolation
+ // For this the relative coordinates of the destination pooint inside the source cell are needed
+ // this cell coordinate system is centered at the middle of the source cell. The source cell has a side lenght of one
+ // add tangential velocity to velocity for coefficient computaion
+ const real cellCoordDestinationX = (coordSourceX[sourceIndex] + rotatedCoordDestinationX) / dx - c1o2;
+ const real cellCoordDestinationY = (coordSourceY[sourceIndex] + rotatedCoordDestinationY) / dx - c1o2;
+ const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] + rotatedCoordDestinationZ) / dx - c1o2;
+ momentsSet.addToVelocity(tangentialVelocitiesX, tangentialVelocitiesY, tangentialVelocitiesZ);
+ vf::lbm::Coefficients coefficients;
+ momentsSet.calculateCoefficients(coefficients, cellCoordDestinationX, cellCoordDestinationY, cellCoordDestinationZ);
+
+
+}
+
+__global__ void updateGlobalCoordinates(
+ unsigned int numberOfNodes,
+ real *globalX,
+ real *globalY,
+ real *globalZ,
+ const real *localX,
+ const real *localY,
+ const real *localZ,
+ real centerCoordX,
+ real centerCoordY,
+ real centerCoordZ,
+ real angleX,
+ real angleY,
+ real angleZ)
+{
+ const unsigned nodeIndex = vf::gpu::getNodeIndex();
+
+ if (nodeIndex >= numberOfNodes) return;
+
+ real globalXtemp;
+ real globalYtemp;
+ real globalZtemp;
+
+ transformRotatingToGlobal(globalXtemp, globalYtemp, globalZtemp, localX[nodeIndex], localY[nodeIndex], localZ[nodeIndex], centerCoordX,
+ centerCoordY, centerCoordZ, angleX, angleY, angleZ);
+
+ globalX[nodeIndex] = globalXtemp;
+ globalY[nodeIndex] = globalYtemp;
+ globalZ[nodeIndex] = globalZtemp;
+}
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticRotating.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
similarity index 86%
rename from src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticRotating.cu
rename to src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
index 3c077f19b48d1230f66d6dbaf1e82d052cc00144..ee96e6b8bdb0f70c801d1709dddfa5a4f5ccd4f7 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticRotating.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/interpolateStaticToRotating.cu
@@ -113,9 +113,9 @@ __global__ void interpolateStaticToRotating(
// 2. calculate the coefficients for the interpolation
// For this the relative coordinates of the destination pooint inside the source cell are needed
// this cell coordinate system is centered at the middle of the source cell. The source cell has a side lenght of one
- const real cellCoordDestinationX = (coordSourceX[sourceIndex] - globalCoordDestinationX) / dx + c1o2;
- const real cellCoordDestinationY = (coordSourceY[sourceIndex] - globalCoordDestinationY) / dx + c1o2;
- const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] - globalCoordDestinationZ) / dx + c1o2;
+ const real cellCoordDestinationX = (coordSourceX[sourceIndex] + globalCoordDestinationX) / dx - c1o2;
+ const real cellCoordDestinationY = (coordSourceY[sourceIndex] + globalCoordDestinationY) / dx - c1o2;
+ const real cellCoordDestinationZ = (coordSourceZ[sourceIndex] + globalCoordDestinationZ) / dx - c1o2;
vf::lbm::Coefficients coefficients;
momentsSet.calculateCoefficients(coefficients, cellCoordDestinationX, cellCoordDestinationY, cellCoordDestinationZ);
@@ -209,7 +209,6 @@ __global__ void interpolateStaticToRotating(
////////////////////////////////////////////////////////////////////////////////
//! - calculate the forces due to the rotating frame of reference
//!
-
// centrifugal force
real forceX = angularVelocityY * angularVelocityY * coordDestinationXlocal +
angularVelocityZ * angularVelocityZ * coordDestinationXlocal -
@@ -232,10 +231,9 @@ __global__ void interpolateStaticToRotating(
////////////////////////////////////////////////////////////////////////////////
//! - subtract half the force from the velocities
//!
-
- vvx -= c1o2 * forceX;
- vvy -= c1o2 * forceY;
- vvz -= c1o2 * forceZ;
+ m100 -= c1o2 * forceX;
+ m010 -= c1o2 * forceY;
+ m001 -= c1o2 * forceZ;
////////////////////////////////////////////////////////////////////////////////
//! - rotate the velocities
@@ -423,82 +421,3 @@ __global__ void interpolateStaticToRotating(
vf::gpu::ListIndices writeIndicesRotating(destinationIndex, neighborXrotating, neighborYrotating, neighborZrotating);
vf::gpu::write(distRoating, writeIndicesRotating, fRotating);
}
-
-__global__ void interpolateRotatingToStatic(
- unsigned int numberOfInterfaceNodes,
- const unsigned int *indicesStatic,
- unsigned int *indicesRotatingCell,
- const real *coordDestinationX,
- const real *coordDestinationY,
- const real *coordDestinationZ,
- const real *coordSourceX,
- const real *coordSourceY,
- const real *coordSourceZ,
- const uint *neighborXrotating,
- const uint *neighborYrotating,
- const uint *neighborZrotating,
- const uint *neighborMMMrotating,
- real centerCoordX,
- real centerCoordY,
- real centerCoordZ,
- real angleX,
- real angleY,
- real angleZ,
- real angularVelocityX,
- real angularVelocityY,
- real angularVelocityZ,
- real dx)
-{
- const unsigned listIndex = vf::gpu::getNodeIndex();
-
- if (listIndex >= numberOfInterfaceNodes) return;
-
- const uint destinationIndex = indicesStatic[listIndex];
- const uint sourceIndex = indicesRotatingCell[listIndex];
- const uint indexNeighborMMMsource = neighborMMMrotating[sourceIndex];
-
- real rotatedCoordDestinationX;
- real rotatedCoordDestinationY;
- real rotatedCoordDestinationZ;
-
- transformGlobalToRotating(rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ,
- coordDestinationX[destinationIndex], coordDestinationY[destinationIndex],
- coordDestinationZ[destinationIndex], centerCoordX, centerCoordY, centerCoordZ, angleX, angleY,
- angleZ);
-
- indicesRotatingCell[listIndex] = traverseSourceCell(
- rotatedCoordDestinationX, rotatedCoordDestinationY, rotatedCoordDestinationZ, indexNeighborMMMsource,
- coordSourceX[indexNeighborMMMsource], coordSourceY[indexNeighborMMMsource], coordSourceZ[indexNeighborMMMsource],
- neighborXrotating, neighborYrotating, neighborZrotating, dx);
-}
-
-__global__ void updateGlobalCoordinates(
- unsigned int numberOfNodes,
- real *globalX,
- real *globalY,
- real *globalZ,
- const real *localX,
- const real *localY,
- const real *localZ,
- real centerCoordX,
- real centerCoordY,
- real centerCoordZ,
- real angleX,
- real angleY,
- real angleZ)
-{
- const unsigned nodeIndex = vf::gpu::getNodeIndex();
-
- if (nodeIndex >= numberOfNodes) return;
-
- real globalXtemp;
- real globalYtemp;
- real globalZtemp;
-
- transformRotatingToGlobal(globalXtemp, globalYtemp, globalZtemp, localX[nodeIndex], localY[nodeIndex], localZ[nodeIndex], centerCoordX,
- centerCoordY, centerCoordZ, angleX, angleY, angleZ);
-
- globalX[nodeIndex] = globalXtemp;
- globalY[nodeIndex] = globalYtemp;
- globalZ[nodeIndex] = globalZtemp;
-}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
index 668305651b2761cfc1644a7f6e13c8afdd1a0507..54bad32e71ebbd75d55fc5409e5e2f767537be7c 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -4113,6 +4113,10 @@ void InterpolateRotatingToStatic(
dim3 threads(parameterDeviceS->numberofthreads, 1, 1);
interpolateRotatingToStatic<<<grid, threads, 0, CU_STREAM_LEGACY>>>(
+ parameterDeviceS->distributions.f[0],
+ parameterDeviceR->distributions.f[0],
+ parameterDeviceS->numberOfNodes,
+ parameterDeviceR->numberOfNodes,
nestedToBase->numberOfCells,
nestedToBase->coarseCellIndices,
nestedToBase->fineCellIndices,
@@ -4135,6 +4139,8 @@ void InterpolateRotatingToStatic(
paraRotDevice->angularVelocity[0],
paraRotDevice->angularVelocity[1],
paraRotDevice->angularVelocity[2],
+ parameterDeviceR->omega,
+ parameterDeviceS->isEvenTimestep,
parameterDeviceS->gridSpacing
);
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
index 56f8bb0055564bdc5b6875fbbe84a1e95d2edc01..99293080c8c277b0e21c0ebc00eac6d70b568a22 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
@@ -44,7 +44,7 @@
#include <lbm/refinement/Coefficients.h>
using namespace vf::basics::constant;
-using namespace vf::lbm::dir;
+using namespace vf::lbm;
namespace vf::gpu
{
@@ -60,7 +60,10 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
real* turbulentViscosity,
unsigned long long numberOfLBnodes,
const real omega,
- bool isEvenTimestep
+ bool isEvenTimestep,
+ real* forcesForAllNodesX = nullptr,
+ real* forcesForAllNodesY = nullptr,
+ real* forcesForAllNodesZ = nullptr
)
{
real omega_ = omega;
@@ -100,7 +103,13 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
real f_fine[27];
read(f_fine, distFine, indices);
- momentsSet.calculateMMM(f_fine, omega_);
+
+ if (forcesForAllNodesX)
+ momentsSet.calculateMMM(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::MMM],
+ forcesForAllNodesY[InterpolationVertex::MMM],
+ forcesForAllNodesY[InterpolationVertex::MMM]);
+ else
+ momentsSet.calculateMMM(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node TSW = MMP
@@ -118,7 +127,12 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculateMMP(f_fine, omega_);
+
+ if (forcesForAllNodesX)
+ momentsSet.calculateMMP(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::MMP],
+ forcesForAllNodesY[InterpolationVertex::MMP], forcesForAllNodesZ[InterpolationVertex::MMP]);
+ else
+ momentsSet.calculateMMP(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node TSE = PMP
@@ -136,7 +150,11 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculatePMP(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculatePMP(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::PMP],
+ forcesForAllNodesY[InterpolationVertex::PMP], forcesForAllNodesZ[InterpolationVertex::PMP]);
+ else
+ momentsSet.calculatePMP(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node BSE = PMM
@@ -154,7 +172,12 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculatePMM(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculatePMM(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::PMM],
+ forcesForAllNodesY[InterpolationVertex::PMM], forcesForAllNodesZ[InterpolationVertex::PMM]);
+ else
+ momentsSet.calculatePMM(f_fine, omega_);
+
//////////////////////////////////////////////////////////////////////////
// source node BNW = MPM
@@ -182,7 +205,11 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculateMPM(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculateMPM(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::MPM],
+ forcesForAllNodesY[InterpolationVertex::MPM], forcesForAllNodesZ[InterpolationVertex::MPM]);
+ else
+ momentsSet.calculateMPM(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node TNW = MPP
@@ -200,7 +227,11 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculateMPP(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculateMPP(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::MPP],
+ forcesForAllNodesY[InterpolationVertex::MPP], forcesForAllNodesZ[InterpolationVertex::MPP]);
+ else
+ momentsSet.calculateMPP(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node TNE = PPP
@@ -218,7 +249,11 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculatePPP(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculatePPP(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::PPP],
+ forcesForAllNodesY[InterpolationVertex::PPP], forcesForAllNodesZ[InterpolationVertex::PPP]);
+ else
+ momentsSet.calculatePPP(f_fine, omega_);
//////////////////////////////////////////////////////////////////////////
// source node BNE = PPM
@@ -232,11 +267,15 @@ template<bool hasTurbulentViscosity> __device__ void calculateMomentSet(
indices.k_M00 = neighborX[k_base_M00];
indices.k_0M0 = k_base_MM0;
indices.k_MM0 = neighborX[k_base_MM0];
-
+
omega_ = hasTurbulentViscosity ? calculateOmega(omega, turbulentViscosity[indices.k_000]) : omega;
read(f_fine, distFine, indices);
- momentsSet.calculatePPM(f_fine, omega_);
+ if (forcesForAllNodesX)
+ momentsSet.calculatePPM(f_fine, omega_, forcesForAllNodesX[InterpolationVertex::PPM],
+ forcesForAllNodesY[InterpolationVertex::PPM], forcesForAllNodesZ[InterpolationVertex::PPM]);
+ else
+ momentsSet.calculatePPM(f_fine, omega_);
}
}
diff --git a/src/lbm/MacroscopicQuantities.h b/src/lbm/MacroscopicQuantities.h
index e1b066a104874a95701ff3a77f6747884c5a54d7..107b588160bf2f06f2b51ffdaf18853f0c2850fd 100644
--- a/src/lbm/MacroscopicQuantities.h
+++ b/src/lbm/MacroscopicQuantities.h
@@ -8,6 +8,7 @@
#define __device__
#endif
+#include <gpu/VirtualFluids_GPU/LBM/LB.h>
#include <basics/DataTypes.h>
#include <basics/constants/NumericConstants.h>
@@ -101,6 +102,57 @@ inline __host__ __device__ void getCompressibleMacroscopicValues(const real *con
getCompressibleMacroscopicValues(f, drho, oneOverRho, vx1, vx2, vx3);
}
+inline __host__ __device__ void getCompressibleMacroscopicValues(const Distributions27 &dist /*[27] for all nodes*/, uint nodeIndex, real &drho, real &vx1, real &vx2, real &vx3)
+{
+ drho = ((dist.f[dir::DIR_PPP][nodeIndex] + dist.f[dir::DIR_MMM][nodeIndex]) +
+ (dist.f[dir::DIR_PMP][nodeIndex] + dist.f[dir::DIR_MPM][nodeIndex])) +
+ ((dist.f[dir::DIR_PMM][nodeIndex] + dist.f[dir::DIR_MPP][nodeIndex]) +
+ (dist.f[dir::DIR_MMP][nodeIndex] + dist.f[dir::DIR_PPM][nodeIndex])) +
+ (((dist.f[dir::DIR_PP0][nodeIndex] + dist.f[dir::DIR_MM0][nodeIndex]) +
+ (dist.f[dir::DIR_PM0][nodeIndex] + dist.f[dir::DIR_MP0][nodeIndex])) +
+ ((dist.f[dir::DIR_P0P][nodeIndex] + dist.f[dir::DIR_M0M][nodeIndex]) +
+ (dist.f[dir::DIR_P0M][nodeIndex] + dist.f[dir::DIR_M0P][nodeIndex])) +
+ ((dist.f[dir::DIR_0PM][nodeIndex] + dist.f[dir::DIR_0MP][nodeIndex]) +
+ (dist.f[dir::DIR_0PP][nodeIndex] + dist.f[dir::DIR_0MM][nodeIndex]))) +
+ ((dist.f[dir::DIR_P00][nodeIndex] + dist.f[dir::DIR_M00][nodeIndex]) +
+ (dist.f[dir::DIR_0P0][nodeIndex] + dist.f[dir::DIR_0M0][nodeIndex]) +
+ (dist.f[dir::DIR_00P][nodeIndex] + dist.f[dir::DIR_00M][nodeIndex])) +
+ dist.f[dir::DIR_000][nodeIndex];
+ real oneOverRho = vf::basics::constant::c1o1 / (drho + vf::basics::constant::c1o1);
+
+ vx1 = ((((dist.f[dir::DIR_PPP][nodeIndex] - dist.f[dir::DIR_MMM][nodeIndex]) +
+ (dist.f[dir::DIR_PMP][nodeIndex] - dist.f[dir::DIR_MPM][nodeIndex])) +
+ ((dist.f[dir::DIR_PMM][nodeIndex] - dist.f[dir::DIR_MPP][nodeIndex]) +
+ (dist.f[dir::DIR_PPM][nodeIndex] - dist.f[dir::DIR_MMP][nodeIndex]))) +
+ (((dist.f[dir::DIR_P0M][nodeIndex] - dist.f[dir::DIR_M0P][nodeIndex]) +
+ (dist.f[dir::DIR_P0P][nodeIndex] - dist.f[dir::DIR_M0M][nodeIndex])) +
+ ((dist.f[dir::DIR_PM0][nodeIndex] - dist.f[dir::DIR_MP0][nodeIndex]) +
+ (dist.f[dir::DIR_PP0][nodeIndex] - dist.f[dir::DIR_MM0][nodeIndex]))) +
+ (dist.f[dir::DIR_P00][nodeIndex] - dist.f[dir::DIR_M00][nodeIndex]));
+ vx1 *= oneOverRho;
+
+ vx2 = ((((dist.f[dir::DIR_PPP][nodeIndex] - dist.f[dir::DIR_MMM][nodeIndex]) +
+ (dist.f[dir::DIR_MPM][nodeIndex] - dist.f[dir::DIR_PMP][nodeIndex])) +
+ ((dist.f[dir::DIR_MPP][nodeIndex] - dist.f[dir::DIR_PMM][nodeIndex]) +
+ (dist.f[dir::DIR_PPM][nodeIndex] - dist.f[dir::DIR_MMP][nodeIndex]))) +
+ (((dist.f[dir::DIR_0PM][nodeIndex] - dist.f[dir::DIR_0MP][nodeIndex]) +
+ (dist.f[dir::DIR_0PP][nodeIndex] - dist.f[dir::DIR_0MM][nodeIndex])) +
+ ((dist.f[dir::DIR_MP0][nodeIndex] - dist.f[dir::DIR_PM0][nodeIndex]) +
+ (dist.f[dir::DIR_PP0][nodeIndex] - dist.f[dir::DIR_MM0][nodeIndex]))) +
+ (dist.f[dir::DIR_0P0][nodeIndex] - dist.f[dir::DIR_0M0][nodeIndex]));
+ vx2 *= oneOverRho;
+
+ vx3 = ((((dist.f[dir::DIR_PPP][nodeIndex] - dist.f[dir::DIR_MMM][nodeIndex]) +
+ (dist.f[dir::DIR_PMP][nodeIndex] - dist.f[dir::DIR_MPM][nodeIndex])) +
+ ((dist.f[dir::DIR_MPP][nodeIndex] - dist.f[dir::DIR_PMM][nodeIndex]) +
+ (dist.f[dir::DIR_MMP][nodeIndex] - dist.f[dir::DIR_PPM][nodeIndex]))) +
+ (((dist.f[dir::DIR_0MP][nodeIndex] - dist.f[dir::DIR_0PM][nodeIndex]) +
+ (dist.f[dir::DIR_0PP][nodeIndex] - dist.f[dir::DIR_0MM][nodeIndex])) +
+ ((dist.f[dir::DIR_M0P][nodeIndex] - dist.f[dir::DIR_P0M][nodeIndex]) +
+ (dist.f[dir::DIR_P0P][nodeIndex] - dist.f[dir::DIR_M0M][nodeIndex]))) +
+ (dist.f[dir::DIR_00P][nodeIndex] - dist.f[dir::DIR_00M][nodeIndex]));
+ vx3 *= oneOverRho;
+}
/*
* Pressure
diff --git a/src/lbm/refinement/Coefficients.h b/src/lbm/refinement/Coefficients.h
index ad506370398e35566435703b598a69e3c3d4f488..e01e183e5b5f9103402994a7fe87e9c4d09dc41f 100644
--- a/src/lbm/refinement/Coefficients.h
+++ b/src/lbm/refinement/Coefficients.h
@@ -50,6 +50,17 @@ using namespace vf::lbm::dir;
namespace vf::lbm
{
+
+enum InterpolationVertex{
+ PPP,
+ MPP,
+ PMP,
+ PPM,
+ MMP,
+ MPM,
+ PMM,
+ MMM
+};
// The Coefficients struct needs be created like this:
// Coefficients coeffs;
@@ -84,7 +95,7 @@ struct MomentsOnSourceNode
real kxxMyyFromfcNEQ;
real kxxMzzFromfcNEQ;
- __host__ __device__ void calculate(const real* const f, const real omega)
+ __host__ __device__ void calculate(const real* const f, const real omega, real forceX, real forceY, real forceZ)
{
// const real f_000 = f[dir::DIR_000];
const real fP00 = f[dir::DIR_P00];
@@ -117,6 +128,10 @@ struct MomentsOnSourceNode
real oneOverRho;
getCompressibleMacroscopicValues(f, this->drho, oneOverRho, this->velocityX, this->velocityY, this->velocityZ);
+ this->velocityX += forceX;
+ this->velocityY += forceY;
+ this->velocityZ += forceZ;
+
////////////////////////////////////////////////////////////////////////////////////
//! - Calculate second order moments for interpolation
//!
@@ -159,44 +174,72 @@ private:
vf::lbm::MomentsOnSourceNode momentsMMM;
public:
- __host__ __device__ void calculatePPP(const real *const f, const real omega)
+ __host__ __device__ void addToVelocity(const real* deltasVelocitiesX, const real* deltasVelocitiesY, const real* deltasVelocitiesZ)
+ {
+ momentsPPP.velocityX += deltasVelocitiesX[InterpolationVertex::PPP];
+ momentsPPP.velocityY += deltasVelocitiesY[InterpolationVertex::PPP];
+ momentsPPP.velocityZ += deltasVelocitiesZ[InterpolationVertex::PPP];
+ momentsMPP.velocityX += deltasVelocitiesX[InterpolationVertex::MPP];
+ momentsMPP.velocityY += deltasVelocitiesY[InterpolationVertex::MPP];
+ momentsMPP.velocityZ += deltasVelocitiesZ[InterpolationVertex::MPP];
+ momentsPMP.velocityX += deltasVelocitiesX[InterpolationVertex::PMP];
+ momentsPMP.velocityY += deltasVelocitiesY[InterpolationVertex::PMP];
+ momentsPMP.velocityZ += deltasVelocitiesZ[InterpolationVertex::PMP];
+ momentsMMP.velocityX += deltasVelocitiesX[InterpolationVertex::MMP];
+ momentsMMP.velocityY += deltasVelocitiesY[InterpolationVertex::MMP];
+ momentsMMP.velocityZ += deltasVelocitiesZ[InterpolationVertex::MMP];
+ momentsPPM.velocityX += deltasVelocitiesX[InterpolationVertex::PPM];
+ momentsPPM.velocityY += deltasVelocitiesY[InterpolationVertex::PPM];
+ momentsPPM.velocityZ += deltasVelocitiesZ[InterpolationVertex::PPM];
+ momentsMPM.velocityX += deltasVelocitiesX[InterpolationVertex::MPM];
+ momentsMPM.velocityY += deltasVelocitiesY[InterpolationVertex::MPM];
+ momentsMPM.velocityZ += deltasVelocitiesZ[InterpolationVertex::MPM];
+ momentsPMM.velocityX += deltasVelocitiesX[InterpolationVertex::PMM];
+ momentsPMM.velocityY += deltasVelocitiesY[InterpolationVertex::PMM];
+ momentsPMM.velocityZ += deltasVelocitiesZ[InterpolationVertex::PMM];
+ momentsMMM.velocityX += deltasVelocitiesX[InterpolationVertex::MMM];
+ momentsMMM.velocityY += deltasVelocitiesY[InterpolationVertex::MMM];
+ momentsMMM.velocityZ += deltasVelocitiesZ[InterpolationVertex::MMM];
+ }
+
+ __host__ __device__ void calculatePPP(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsPPP.calculate(f, omega);
+ momentsPPP.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculateMPP(const real *const f, const real omega)
+ __host__ __device__ void calculateMPP(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsMPP.calculate(f, omega);
+ momentsMPP.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculatePMP(const real *const f, const real omega)
+ __host__ __device__ void calculatePMP(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsPMP.calculate(f, omega);
+ momentsPMP.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculateMMP(const real *const f, const real omega)
+ __host__ __device__ void calculateMMP(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsMMP.calculate(f, omega);
+ momentsMMP.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculatePPM(const real *const f, const real omega)
+ __host__ __device__ void calculatePPM(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsPPM.calculate(f, omega);
+ momentsPPM.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculateMPM(const real *const f, const real omega)
+ __host__ __device__ void calculateMPM(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsMPM.calculate(f, omega);
+ momentsMPM.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculatePMM(const real *const f, const real omega)
+ __host__ __device__ void calculatePMM(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsPMM.calculate(f, omega);
+ momentsPMM.calculate(f, omega, forceX, forceY, forceZ);
}
- __host__ __device__ void calculateMMM(const real *const f, const real omega)
+ __host__ __device__ void calculateMMM(const real *const f, const real omega, real forceX = 0.0, real forceY = 0.0, real forceZ = 0.0)
{
- momentsMMM.calculate(f, omega);
+ momentsMMM.calculate(f, omega, forceX, forceY, forceZ);
}
__host__ __device__ void calculateCoefficients(Coefficients &coefficients, real xoff, real yoff, real zoff) const