Merge pull request #2406 from su2code/better_fea_symmetry

pcarruscag · web-flow · commit 99e6e5484dcd · 2024-12-28T14:08:25.000-08:00
Remove axis-perpendicular restriction on FEA symmetry planes
diff --git a/Common/include/linear_algebra/CSysMatrix.hpp b/Common/include/linear_algebra/CSysMatrix.hpp
@@ -805,10 +805,10 @@ class CSysMatrix {
   void EnforceSolutionAtNode(unsigned long node_i, const OtherType* x_i, CSysVector<OtherType>& b);
 
   /*!
-   * \brief Version of EnforceSolutionAtNode for a single degree of freedom.
+   * \brief Similar to EnforceSolutionAtNode, but for 0 projection in a given direction.
    */
   template <class OtherType>
-  void EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i, CSysVector<OtherType>& b);
+  void EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector<OtherType>& b);
 
   /*!
    * \brief Sets the diagonal entries of the matrix as the sum of the blocks in the corresponding column.
diff --git a/Common/src/linear_algebra/CSysMatrix.cpp b/Common/src/linear_algebra/CSysMatrix.cpp
@@ -1027,36 +1027,59 @@ void CSysMatrix<ScalarType>::EnforceSolutionAtNode(const unsigned long node_i, c
 
 template <class ScalarType>
 template <class OtherType>
-void CSysMatrix<ScalarType>::EnforceSolutionAtDOF(unsigned long node_i, unsigned long iVar, OtherType x_i,
-                                                  CSysVector<OtherType>& b) {
+void CSysMatrix<ScalarType>::EnforceZeroProjection(unsigned long node_i, const OtherType* n, CSysVector<OtherType>& b) {
   for (auto index = row_ptr[node_i]; index < row_ptr[node_i + 1]; ++index) {
     const auto node_j = col_ind[index];
 
-    /*--- Delete row iVar of block j on row i (bij) and ATTEMPT
-     *    to delete column iVar block i on row j (bji). ---*/
+    /*--- Remove product components of block j on row i (bij) and ATTEMPT
+     *    to remove solution components of block i on row j (bji).
+     *    This is identical to symmetry correction applied to gradients
+     *    but extended to the entire matrix. ---*/
 
     auto bij = &matrix[index * nVar * nVar];
     auto bji = GetBlock(node_j, node_i);
 
-    /*--- The "attempt" part. ---*/
+    /*--- Attempt to remove solution components. ---*/
+    ScalarType nbn{};
     if (bji != nullptr) {
-      for (auto jVar = 0ul; jVar < nVar; ++jVar) {
-        /*--- Column product. ---*/
-        b[node_j * nVar + jVar] -= bji[jVar * nVar + iVar] * x_i;
-        /*--- Delete entries. ---*/
-        bji[jVar * nVar + iVar] = 0.0;
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        ScalarType proj{};
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          proj += bji[iVar * nVar + jVar] * PassiveAssign(n[jVar]);
+        }
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          bji[iVar * nVar + jVar] -= proj * PassiveAssign(n[jVar]);
+        }
+        nbn += proj * PassiveAssign(n[iVar]);
       }
     }
 
-    /*--- Delete row. ---*/
-    for (auto jVar = 0ul; jVar < nVar; ++jVar) bij[iVar * nVar + jVar] = 0.0;
+    /*--- Product components. ---*/
+    for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+      ScalarType proj{};
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        proj += bij[iVar * nVar + jVar] * PassiveAssign(n[iVar]);
+      }
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        bij[iVar * nVar + jVar] -= proj * PassiveAssign(n[iVar]);
+      }
+    }
 
-    /*--- Set the diagonal entry of the block to 1. ---*/
-    if (node_j == node_i) bij[iVar * (nVar + 1)] = 1.0;
+    /*--- This part doesn't have the "*2" factor because the product components
+     *    were removed from the result of removing the solution components
+     *    instead of from the original block (bji == bij). ---*/
+    if (node_i == node_j) {
+      for (auto iVar = 0ul; iVar < nVar; ++iVar) {
+        for (auto jVar = 0ul; jVar < nVar; ++jVar) {
+          bij[iVar * nVar + jVar] += PassiveAssign(n[iVar]) * nbn * PassiveAssign(n[jVar]);
+        }
+      }
+    }
   }
 
-  /*--- Set known solution in rhs vector. ---*/
-  b(node_i, iVar) = x_i;
+  OtherType proj{};
+  for (auto iVar = 0ul; iVar < nVar; ++iVar) proj += b(node_i, iVar) * n[iVar];
+  for (auto iVar = 0ul; iVar < nVar; ++iVar) b(node_i, iVar) -= proj * n[iVar];
 }
 
 template <class ScalarType>
@@ -1203,8 +1226,7 @@ void CSysMatrix<ScalarType>::ComputePastixPreconditioner(const CSysVector<Scalar
 #define INSTANTIATE_MATRIX(TYPE)                                                                                  \
   template class CSysMatrix<TYPE>;                                                                                \
   template void CSysMatrix<TYPE>::EnforceSolutionAtNode(unsigned long, const su2double*, CSysVector<su2double>&); \
-  template void CSysMatrix<TYPE>::EnforceSolutionAtDOF(unsigned long, unsigned long, su2double,                   \
-                                                       CSysVector<su2double>&);                                   \
+  template void CSysMatrix<TYPE>::EnforceZeroProjection(unsigned long, const su2double*, CSysVector<su2double>&); \
   INSTANTIATE_COMMS(TYPE)
 
 #ifdef CODI_FORWARD_TYPE
diff --git a/SU2_CFD/src/solvers/CFEASolver.cpp b/SU2_CFD/src/solvers/CFEASolver.cpp
@@ -1588,6 +1588,20 @@ void CFEASolver::BC_Clamped_Post(CGeometry *geometry, const CConfig *config, uns
 
 }
 
+namespace {
+/*--- Helper for BC_Sym_Plane ---*/
+template <typename Read, typename Write>
+void SubtractProjection(unsigned short nDim, const su2double* n, const Read& read, const Write& write) {
+  su2double tmp[3] = {};
+  for (auto iDim = 0u; iDim < nDim; ++iDim) tmp[iDim] = read(iDim);
+  const su2double proj = DotProduct(nDim, tmp, n);
+  for (auto iDim = 0u; iDim < nDim; ++iDim) {
+    tmp[iDim] -= proj * n[iDim];
+    write(iDim, tmp[iDim]);
+  }
+}
+}
+
 void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsigned short val_marker) {
 
   if (geometry->GetnElem_Bound(val_marker) == 0) return;
@@ -1610,14 +1624,8 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsign
     case 3: TriangleNormal(nodeCoord, normal); break;
     case 4: QuadrilateralNormal(nodeCoord, normal); break;
   }
-
-  auto axis = 0u;
-  for (auto iDim = 1u; iDim < MAXNDIM; ++iDim)
-    axis = (fabs(normal[iDim]) > fabs(normal[axis]))? iDim : axis;
-
-  if (fabs(normal[axis]) < 0.99*Norm(int(MAXNDIM),normal)) {
-    SU2_MPI::Error("The structural solver only supports axis-aligned symmetry planes.",CURRENT_FUNCTION);
-  }
+  const su2double area = Norm(int(MAXNDIM), normal);
+  for (auto iDim = 0u; iDim < MAXNDIM; ++iDim) normal[iDim] /= area;
 
   /*--- Impose zero displacement perpendicular to the symmetry plane. ---*/
 
@@ -1627,20 +1635,29 @@ void CFEASolver::BC_Sym_Plane(CGeometry *geometry, const CConfig *config, unsign
     const auto iPoint = geometry->vertex[val_marker][iVertex]->GetNode();
 
     /*--- Set and enforce solution at current and previous time-step ---*/
-    nodes->SetSolution(iPoint, axis, 0.0);
+#define SUBTRACT_PROJECTION(READ, WRITE)                                                                \
+    SubtractProjection(nDim, normal, [&](unsigned short iDim) { return nodes->READ(iPoint, iDim); },    \
+                       [&](unsigned short iDim, const su2double& x) { nodes->WRITE(iPoint, iDim, x); });
+    SUBTRACT_PROJECTION(GetSolution, SetSolution)
     if (dynamic) {
-      nodes->SetSolution_Vel(iPoint, axis, 0.0);
-      nodes->SetSolution_Accel(iPoint, axis, 0.0);
-      nodes->Set_Solution_time_n(iPoint, axis, 0.0);
-      nodes->SetSolution_Vel_time_n(iPoint, axis, 0.0);
-      nodes->SetSolution_Accel_time_n(iPoint, axis, 0.0);
+      SUBTRACT_PROJECTION(GetSolution_Vel, SetSolution_Vel)
+      SUBTRACT_PROJECTION(GetSolution_Accel, SetSolution_Accel)
+      SUBTRACT_PROJECTION(GetSolution_time_n, Set_Solution_time_n)
+      SUBTRACT_PROJECTION(GetSolution_Vel_time_n, SetSolution_Vel_time_n)
+      SUBTRACT_PROJECTION(GetSolution_Accel_time_n, SetSolution_Accel_time_n)
     }
 
     /*--- Set and enforce 0 solution for mesh deformation ---*/
-    nodes->SetBound_Disp(iPoint, axis, 0.0);
-    LinSysSol(iPoint, axis) = 0.0;
-    if (LinSysReact.GetLocSize() > 0) LinSysReact(iPoint, axis) = 0.0;
-    Jacobian.EnforceSolutionAtDOF(iPoint, axis, su2double(0.0), LinSysRes);
+    SUBTRACT_PROJECTION(GetBound_Disp, SetBound_Disp)
+#undef SUBTRACT_PROJECTION
+
+    SubtractProjection(nDim, normal, [&](unsigned short iDim) { return LinSysSol(iPoint, iDim); },
+                       [&](unsigned short iDim, const su2double& x) { LinSysSol(iPoint, iDim) = x; });
+    if (LinSysReact.GetLocSize() > 0) {
+      SubtractProjection(nDim, normal, [&](unsigned short iDim) { return LinSysReact(iPoint, iDim); },
+                         [&](unsigned short iDim, const su2double& x) { LinSysReact(iPoint, iDim) = x; });
+    }
+    Jacobian.EnforceZeroProjection(iPoint, normal, LinSysRes);
 
   }
 
diff --git a/TestCases/parallel_regression.py b/TestCases/parallel_regression.py
@@ -1233,7 +1233,7 @@ def main():
     # For a thin disk with the inner and outer radius of this geometry, from
     # "Formulas for Stress, Strain, and Structural Matrices", 2nd Edition, figure 19-4,
     # the maximum stress is 165.6MPa, we get a Von Misses stress very close to that.
-    rotating_cylinder_fea.test_vals = [-6.005467, -5.615543, -5.615527, 38, -8.126591, 1.6457e8]
+    rotating_cylinder_fea.test_vals = [-6.861940, -6.835545, -6.895500, 22, -8.313847, 1.6502e+08]
     test_list.append(rotating_cylinder_fea)
 
     # Dynamic beam, 2d
