diff --git a/Common/include/option_structure.inl b/Common/include/option_structure.inl index fe1b2df385ee..123edcb4ba76 100644 --- a/Common/include/option_structure.inl +++ b/Common/include/option_structure.inl @@ -28,6 +28,7 @@ #include "option_structure.hpp" #include "parallelization/mpi_structure.hpp" +#include "toolboxes/geometry_toolbox.hpp" using namespace std; template @@ -1767,15 +1768,33 @@ class COptionPeriodic : public COptionBase { translation[i][1] = translation[i + nVals / 2][1] = getval(i, 9); translation[i][2] = translation[i + nVals / 2][2] = getval(i, 10); - /*--- Mirror the rotational angles and translation vector (rotational center does not need to move). ---*/ - rot_angles[i + nVals / 2][0] *= -1; - rot_angles[i + nVals / 2][1] *= -1; - rot_angles[i + nVals / 2][2] *= -1; - translation[i + nVals / 2][0] *= -1; - translation[i + nVals / 2][1] *= -1; - translation[i + nVals / 2][2] *= -1; - if (err) return badValue("periodic", name); + + /*--- Invert x' = R (x - c) + c + t. Negating Euler angles in the + * same order is only correct for rotations about a single axis. ---*/ + const auto donor = i + nVals / 2; + su2double rotation[3][3]; + GeometryToolbox::RotationMatrix(rot_angles[i][0], rot_angles[i][1], rot_angles[i][2], rotation); + for (auto iDim = 0u; iDim < 3; ++iDim) { + rot_angles[donor][iDim] = -rot_angles[i][iDim]; + translation[donor][iDim] = 0; + for (auto jDim = 0u; jDim < 3; ++jDim) translation[donor][iDim] -= rotation[jDim][iDim] * translation[i][jDim]; + } + + /*--- Keep single-axis angles unchanged, including their AD recording. ---*/ + const auto nAngles = (rot_angles[i][0] != 0) + (rot_angles[i][1] != 0) + (rot_angles[i][2] != 0); + if (nAngles > 1) { + const auto cosPhi = sqrt(rotation[0][0] * rotation[0][0] + rotation[0][1] * rotation[0][1]); + rot_angles[donor][1] = atan2(-rotation[0][2], cosPhi); + if (cosPhi > 1e-12) { + rot_angles[donor][0] = atan2(rotation[1][2], rotation[2][2]); + rot_angles[donor][2] = atan2(rotation[0][1], rotation[0][0]); + } else { + /*--- At gimbal lock only the combined x/z rotation is determined. ---*/ + rot_angles[donor][0] = atan2(-rotation[2][1], rotation[1][1]); + rot_angles[donor][2] = 0; + } + } } return ""; diff --git a/Common/src/fem/fem_geometry_structure.cpp b/Common/src/fem/fem_geometry_structure.cpp index ae71919060a2..2d097d4e4c01 100644 --- a/Common/src/fem/fem_geometry_structure.cpp +++ b/Common/src/fem/fem_geometry_structure.cpp @@ -34,6 +34,7 @@ #include "../../include/geometry/primal_grid/CPrimalGridBoundFEM.hpp" #include "../../include/adt/CADTElemClass.hpp" #include "../../include/adt/CADTPointsOnlyClass.hpp" +#include "../../include/toolboxes/geometry_toolbox.hpp" /* Prototypes for Lapack functions, if MKL or LAPACK is used. */ #if defined(HAVE_MKL) || defined(HAVE_LAPACK) @@ -1744,31 +1745,15 @@ CMeshFEM::CMeshFEM(CGeometry* geometry, CConfig* config) { transformation from the donor. This is the transpose of the transformation to the donor. ---*/ - /* Store (center-trans) as it is constant and will be added on. */ - su2double translation[] = {center[0] - trans[0], center[1] - trans[1], center[2] - trans[2]}; - - /* Store angles separately for clarity. Compute sines/cosines. */ - su2double theta = angles[0]; - su2double phi = angles[1]; - su2double psi = angles[2]; - - su2double cosTheta = cos(theta), cosPhi = cos(phi), cosPsi = cos(psi); - su2double sinTheta = sin(theta), sinPhi = sin(phi), sinPsi = sin(psi); - - /* Compute the rotation matrix. Note that the implicit - ordering is rotation about the x-axis, y-axis, then z-axis. */ su2double rotMatrix[3][3]; - rotMatrix[0][0] = cosPhi * cosPsi; - rotMatrix[0][1] = cosPhi * sinPsi; - rotMatrix[0][2] = -sinPhi; - - rotMatrix[1][0] = sinTheta * sinPhi * cosPsi - cosTheta * sinPsi; - rotMatrix[1][1] = sinTheta * sinPhi * sinPsi + cosTheta * cosPsi; - rotMatrix[1][2] = sinTheta * cosPhi; - - rotMatrix[2][0] = cosTheta * sinPhi * cosPsi + sinTheta * sinPsi; - rotMatrix[2][1] = cosTheta * sinPhi * sinPsi - sinTheta * cosPsi; - rotMatrix[2][2] = cosTheta * cosPhi; + GeometryToolbox::RotationMatrix(angles[0], angles[1], angles[2], rotMatrix); + for (auto i = 0u; i < 3; ++i) + for (auto j = 0u; j < i; ++j) std::swap(rotMatrix[i][j], rotMatrix[j][i]); + + /*--- Invert the complete affine map: c + R^T (x-c-t). ---*/ + su2double translation[] = {center[0], center[1], center[2]}; + for (auto i = 0u; i < 3; ++i) + for (auto j = 0u; j < 3; ++j) translation[i] -= rotMatrix[i][j] * trans[j]; /* Loop over the halo points for this periodic transformation. */ for (unsigned long i = iLow; i < iUpp; ++i) { diff --git a/TestCases/hom_euler/periodic_affine/periodic_affine.cfg b/TestCases/hom_euler/periodic_affine/periodic_affine.cfg new file mode 100644 index 000000000000..1d32b9f9a6e7 --- /dev/null +++ b/TestCases/hom_euler/periodic_affine/periodic_affine.cfg @@ -0,0 +1,21 @@ +SOLVER= FEM_EULER +MESH_FORMAT= BOX +MESH_BOX_SIZE= (3,3,3) +MESH_BOX_LENGTH= (1,1,1) +MESH_BOX_OFFSET= (0,0,0) +MESH_BOX_POLY_SOL_FEM= 1 +MACH_NUMBER= 0.1 +AOA= 0 +FREESTREAM_PRESSURE= 101325 +FREESTREAM_TEMPERATURE= 300 +MARKER_EULER= (x_plus,y_minus) +MARKER_FAR= (z_minus,z_plus) +NUM_METHOD_FEM_FLOW= DG +RIEMANN_SOLVER_FEM= ROE +TIME_DISCRE_FEM_FLOW= RUNGE-KUTTA_EXPLICIT +CFL_NUMBER= 0.1 +ITER= 2 +OUTPUT_FILES= NONE +SCREEN_OUTPUT= (INNER_ITER,RMS_RES) +HISTORY_OUTPUT= (ITER,RMS_RES) +MARKER_PERIODIC= (x_minus,y_plus, 0,0,0, 0,0,90, 1,1,0) diff --git a/TestCases/serial_regression.py b/TestCases/serial_regression.py index e1a826cbf170..2c8f8ada3807 100755 --- a/TestCases/serial_regression.py +++ b/TestCases/serial_regression.py @@ -554,6 +554,15 @@ def main(): fem_euler_naca0012.test_vals = [-6.519946, -5.976944, 0.255551, 0.000028] test_list.append(fem_euler_naca0012) + # Periodic DG affine map: inverse translation must be rotated. + fem_periodic_affine = TestCase('fem_periodic_affine') + fem_periodic_affine.cfg_dir = "hom_euler/periodic_affine" + fem_periodic_affine.cfg_file = "periodic_affine.cfg" + fem_periodic_affine.test_iter = 1 + fem_periodic_affine.ntest_vals = 5 + fem_periodic_affine.test_vals = [2.083795, 4.631729, 4.386282, 3.676592, 7.559938] + test_list.append(fem_periodic_affine) + ############################ ### DG-FEM Navier-Stokes ### ############################ diff --git a/UnitTests/Common/CConfig_tests.cpp b/UnitTests/Common/CConfig_tests.cpp index 111918785545..1a8a01b1c992 100644 --- a/UnitTests/Common/CConfig_tests.cpp +++ b/UnitTests/Common/CConfig_tests.cpp @@ -29,6 +29,7 @@ #include #include #include "../../Common/include/CConfig.hpp" +#include "../../Common/include/toolboxes/geometry_toolbox.hpp" namespace { @@ -87,3 +88,32 @@ TEST_CASE("INIT_OPTION_INC defaults", "[Config]") { INIT_OPTION_INC::OPERATING_PRESSURE); CHECK(GetInitOptionInc(ideal_gas_options + "INIT_OPTION_INC= DENSITY_INIT\n") == INIT_OPTION_INC::DENSITY_INIT); } + +TEST_CASE("Periodic affine transformation and its inverse", "[Config][Periodic]") { + const auto angles = GENERATE(std::string("0,0,0"), std::string("0,0,30"), std::string("20,30,10"), + std::string("90,0,90"), std::string("-90,0,90")); + std::stringstream options(base_options + "MARKER_PERIODIC= (a,b, 2,-3,1, " + angles + ", 1,2,-3)\n"); + auto* original = std::cout.rdbuf(nullptr); + CConfig config(options, SU2_COMPONENT::SU2_CFD, false); + std::cout.rdbuf(original); + su2double rotations[2][3][3]; + const std::string markers[] = {"a", "b"}; + for (auto side = 0u; side < 2; ++side) { + const auto* rotation = config.GetPeriodicRotAngles(markers[side]); + GeometryToolbox::RotationMatrix(rotation[0], rotation[1], rotation[2], rotations[side]); + } + su2double error = 0; + for (auto iDim = 0u; iDim < 3; ++iDim) { + for (auto jDim = 0u; jDim < 3; ++jDim) { + su2double product = 0; + for (auto kDim = 0u; kDim < 3; ++kDim) product += rotations[1][iDim][kDim] * rotations[0][kDim][jDim]; + error = std::max(error, fabs(product - (iDim == jDim))); + } + const auto* forward = config.GetPeriodicTranslation("a"); + su2double translation = config.GetPeriodicTranslation("b")[iDim]; + for (auto jDim = 0u; jDim < 3; ++jDim) translation += rotations[1][iDim][jDim] * forward[jDim]; + error = std::max(error, fabs(translation)); + } + INFO("Euler angles: " << angles); + CHECK(error < 1e-12); +}