Added prevention of intersections in grid elements after deformation

Common/include/CConfig.hpp

@@ -108,8 +108,14 @@ class CConfig {
   unsigned short FFD_Blending;    /*!< \brief Kind of FFD Blending function. */
   su2double* FFD_BSpline_Order;   /*!< \brief BSpline order in i,j,k direction. */
   su2double FFD_Tol;              /*!< \brief Tolerance in the point inversion problem. */
-  su2double Opt_RelaxFactor;      /*!< \brief Scale factor for the line search. */
-  su2double Opt_LineSearch_Bound; /*!< \brief Bounds for the line search. */
+  bool FFD_IntPrev;                       /*!< \brief Enables self-intersection prevention procedure within the FFD box. */
+  unsigned short FFD_IntPrev_MaxIter;     /*!< \brief Amount of iterations for FFD box self-intersection prevention procedure. */
+  unsigned short FFD_IntPrev_MaxDepth;    /*!< \brief Maximum recursion depth for FFD box self-intersection procedure. */
+  bool ConvexityCheck;                    /*!< \brief Enables convexity check on all mesh elements. */
+  unsigned short ConvexityCheck_MaxIter;  /*!< \brief Amount of iterations for convexity check in deformations. */
+  unsigned short ConvexityCheck_MaxDepth; /*!< \brief Maximum recursion depth for convexity check in deformations.*/
+  su2double Opt_RelaxFactor;              /*!< \brief Scale factor for the line search. */
+  su2double Opt_LineSearch_Bound;         /*!< \brief Bounds for the line search. */
   su2double StartTime;
   bool ContinuousAdjoint,   /*!< \brief Flag to know if the code is solving an adjoint problem. */
   Viscous,                  /*!< \brief Flag to know if the code is solving a viscous problem. */
@@ -8086,6 +8092,27 @@ class CConfig {
    */
   su2double GetFFD_Tol(void) const { return FFD_Tol; }
 
+  /*!
+   * \brief Get information about whether to do a check on self-intersections within 
+      the FFD box based on value on the Jacobian determinant.
+   * \param[out] FFD_IntPrev: <code>TRUE</code> if FFD intersection prevention is active; otherwise <code>FALSE</code>.
+   * \param[out] FFD_IntPrev_MaxIter: Maximum number of iterations in the intersection prevention procedure.
+   * \param[out] FFD_IntPrev_MaxDepth: Maximum recursion depth in the intersection prevention procedure.
+   */
+  tuple<bool, unsigned short, unsigned short> GetFFD_IntPrev(void) {
+    return make_tuple(FFD_IntPrev, FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth);
+  }
+
+  /*!
+   * \brief Get information about whether to do a check on convexity of the mesh elements.
+   * \param[out] ConvexityCheck: <code>TRUE</code> if convexity check is active; otherwise <code>FALSE</code>.
+   * \param[out] ConvexityCheck_MaxIter: Maximum number of iterations in the convexity check.
+   * \param[out] ConvexityCheck_MaxDepth: Maximum recursion depth in the convexity check.
+   */
+  tuple<bool, unsigned short, unsigned short> GetConvexityCheck(void) {
+    return make_tuple(ConvexityCheck, ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth);
+  }
+
   /*!
    * \brief Get the scale factor for the line search.
    * \return Scale factor for the line search.

Common/include/geometry/CGeometry.hpp

@@ -109,7 +109,8 @@ class CGeometry {
   nelem_triangle_bound{0},        /*!< \brief Number of triangles on the mesh boundaries. */
   Global_nelem_triangle_bound{0}, /*!< \brief Total number of triangles on the mesh boundaries across all processors. */
   nelem_quad_bound{0},            /*!< \brief Number of quads on the mesh boundaries. */
-  Global_nelem_quad_bound{0};     /*!< \brief Total number of quads on the mesh boundaries across all processors. */
+  Global_nelem_quad_bound{0},     /*!< \brief Total number of quads on the mesh boundaries across all processors. */
+  nNonconvexElements{0};          /*!< \brief Number of nonconvex elements in the mesh. */
 
   unsigned short nDim{0};         /*!< \brief Number of dimension of the problem. */
   unsigned short nZone{0};        /*!< \brief Number of zones in the problem. */
@@ -1699,5 +1700,16 @@ class CGeometry {
    */
   static void ComputeWallDistance(const CConfig * const *config_container, CGeometry ****geometry_container);
 
+  /*!
+   * \brief Set the amount of nonconvex elements in the mesh.
+   * \param[in] nonconvex_elems - amount of nonconvex elements in the mesh
+   */
+  void SetnNonconvexElements(unsigned long nonconvex_elems) {nNonconvexElements = nonconvex_elems;}
+
+  /*!
+   * \brief Get the amount of nonconvex elements in the mesh.
+   * \param[out] nNonconvexElements- amount of nonconvex elements in the mesh
+   */
+  unsigned long GetnNonconvexElements(void) {return nNonconvexElements;}
 };
 

Common/include/grid_movement/CSurfaceMovement.hpp

@@ -211,7 +211,7 @@ class CSurfaceMovement : public CGridMovement {
    * \param[in] geometry - Geometrical definition of the problem.
    * \param[in] config - Definition of the particular problem.
    */
-  void SetSurface_Deformation(CGeometry *geometry, CConfig *config) override;
+  void SetSurface_Deformation(CGeometry *geometry, CConfig *config, su2double** totaldeformation = nullptr);
 
   /*!
    * \brief Compute the parametric coordinates of a grid point using a point inversion strategy
@@ -268,6 +268,15 @@ class CSurfaceMovement : public CGridMovement {
    */
   void GetCartesianCoordCP(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBoxParent, CFreeFormDefBox *FFDBoxChild);
 
+  /*!
+   * \brief Apply the design variables to the control point position
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] FFDBox - Array with all the free forms FFDBoxes of the computation.
+   * \param[in] iFFDBox - Index of FFD box.
+   */
+  void ApplyDesignVariables(CGeometry *geometry, CConfig *config, CFreeFormDefBox **FFDBox, unsigned short iFFDBox);
+
   /*!
    * \brief Recompute the cartesian coordinates using the control points position.
    * \param[in] geometry - Geometrical definition of the problem.
@@ -484,4 +493,13 @@ class CSurfaceMovement : public CGridMovement {
    * \param[in] config - Definition of the particular problem.
    */
   void SetSurface_Derivative(CGeometry *geometry, CConfig *config);
+
+  /*!
+   * \brief Calculate the determinant of the Jacobian matrix for the FFD problem.
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] FFDBox - Free form deformation box.
+   * \return Number of points with negative Jacobian determinant.
+   */
+  unsigned long calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox);
 };

Common/include/grid_movement/CVolumetricMovement.hpp

@@ -240,6 +240,14 @@ class CVolumetricMovement : public CGridMovement {
    */
   void ComputeDeforming_Element_Volume(CGeometry *geometry, su2double &MinVolume, su2double &MaxVolume, bool Screen_Output);
 
+  /*!
+   * \brief Compute amount of nonconvex elements
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] Screen_Output - determines if text is written to screen
+   */
+  void ComputenNonconvexElements(CGeometry *geometry, bool Screen_Output);
+
+
   /*!
    * \brief Compute the minimum distance to the nearest solid surface.
    * \param[in] geometry - Geometrical definition of the problem.

Common/src/CConfig.cpp

@@ -2564,6 +2564,24 @@ void CConfig::SetConfig_Options() {
   /* DESCRIPTION: Free surface damping coefficient */
   addDoubleOption("FFD_TOLERANCE", FFD_Tol, 1E-10);
 
+  /* DESCRIPTION: Procedure to prevent self-intersections within the FFD box based on Jacobian determinant */
+  addBoolOption("FFD_INTPREV", FFD_IntPrev, NO);
+
+  /* DESCRIPTION: Number of total iterations in the convexity check procedure */
+  addUnsignedShortOption("FFD_INTPREV_ITER", FFD_IntPrev_MaxIter, 10);
+
+  /* DESCRIPTION: Recursion depth in the FFD self-intersection prevention */
+  addUnsignedShortOption("FFD_INTPREV_DEPTH", FFD_IntPrev_MaxDepth, 3);
+
+  /* DESCRIPTION: Convexity check on all mesh elements */
+  addBoolOption("CONVEXITY_CHECK", ConvexityCheck, NO);
+
+  /* DESCRIPTION: Number of total iterations in the convexity check procedure */
+  addUnsignedShortOption("CONVEXITY_CHECK_ITER", ConvexityCheck_MaxIter, 10);
+
+  /* DESCRIPTION: Recursion depth in the FFD self-intersection prevention */
+  addUnsignedShortOption("CONVEXITY_CHECK_DEPTH", ConvexityCheck_MaxDepth, 3);
+
   /* DESCRIPTION: Definition of the FFD boxes */
   addFFDDefOption("FFD_DEFINITION", nFFDBox, CoordFFDBox, TagFFDBox);