Fix inverse periodic affine transforms in the parser and FEM

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 <class Tenum, class TField>
@@ -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 "";

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) {

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)

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 ###
     ############################

UnitTests/Common/CConfig_tests.cpp

@@ -29,6 +29,7 @@
 #include <sstream>
 #include <string>
 #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);
+}