Added routines to prevent intersections in grid elements after free form deformation.

After deformation, the surface points are checked for self-intersection within the FFD box. If that is the case, a recursive procedure is done to prevent those. After grid deformation, all elements are checked on convexity. If there are nonconvex elements, a procedure is started that starts from the undeformed mesh and recursively changes the deformation magnitude until no nonconvex elements are present.

Author: tollennaert

Common/include/CConfig.hpp

@@ -108,6 +108,13 @@ 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. */
+  bool FFD_IntPrev;               /*!< \brief Enables self-intersection prevention procedure within the FFD box. */
+  unsigned short FFD_IntPrev_Iter; /*!< \brief Amount of iterations for FFD box self-intersection prevention procedure. */
+  unsigned short FFD_IntPrev_Depth; /*!< \brief Maximum recursion depth for FFD box self-intersection procedure. */
+  bool ConvexityCheck;            /*!< \brief Enables convexity check on all mesh elements. */
+  unsigned short ConvexityCheck_Iter; /*!< \brief Amount of iterations for convexity check in deformations. */
+  unsigned short ConvexityCheck_Depth; /*!< \brief Maximum recursion depth for convexity check in deformations.*/
+  su2double** TotalDeformation;   /*!< \brief Total deformation of the design variables. */
   su2double Opt_RelaxFactor;      /*!< \brief Scale factor for the line search. */
   su2double Opt_LineSearch_Bound; /*!< \brief Bounds for the line search. */
   su2double StartTime;
@@ -8171,6 +8178,72 @@ 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 of the Jacobian determinant.
+   * \return <code>TRUE</code> if FFD intersection prevention is active; otherwise <code>FALSE</code>.
+   */
+  bool GetFFD_IntPrev(void) const { return FFD_IntPrev; }
+
+  /*!
+   * \brief Get the maximum number of iterations in the recursive convexity check procedure in deformation.
+   * \return Maximum number of iterations.
+   */
+  unsigned short GetFFD_IntPrev_Iter(void) const { return FFD_IntPrev_Iter; }
+
+  /*!
+   * \brief Get the maximum recursion depth in the recursive convexity check procedure in deformation.
+   * \return Maximum recursion depth.
+   */
+  unsigned short GetFFD_IntPrev_Depth(void) const { return FFD_IntPrev_Depth; }
+
+  /*!
+   * \brief Get information about whether to do a check on convexity of the mesh elements.
+   * \return <code>TRUE</code> if convexity check is active; otherwise <code>FALSE</code>.
+   */
+  bool GetConvexity_Check(void) const { return ConvexityCheck; }
+
+  /*!
+   * \brief Get the maximum number of iterations in the recursive convexity check procedure in deformation.
+   * \return Maximum number of iterations.
+   */
+  unsigned short GetConvexityCheck_Iter(void) const { return ConvexityCheck_Iter; }
+
+  /*!
+   * \brief Get the maximum recursion depth in the recursive convexity check procedure in deformation.
+   * \return Maximum recursion depth.
+   */
+  unsigned short GetConvexityCheck_Depth(void) const { return ConvexityCheck_Depth; }
+
+  /*!
+   * \brief Set the initial values for the total deformation of the design variables.
+   * \param[in] totalDef - total deformation of all design variables
+   */
+  void Initialize_TotalDeformation(su2double** totalDef) { TotalDeformation = totalDef; }
+
+  /*!
+   * \brief Set the total deformation value of each design variable.
+   * \param[in] val_dv - design variable
+   * \param[in] val_ind - index of design variable
+   * \param[in] val - value of deformation magnitude
+   */
+  void SetTotalDeformation(unsigned short val_dv, unsigned short val_ind, su2double val) {TotalDeformation[val_dv][val_ind] = val; }
+
+  /*!
+   * \brief Add a value to the total deformation value of each design variable.
+   * \param[in] val_dv - design variable
+   * \param[in] val_ind - index of design variable
+   * \param[in] val - value of deformation magnitude that is added
+   */
+  void AddTotalDeformation(unsigned short val_dv, unsigned short val_ind, su2double val) {TotalDeformation[val_dv][val_ind] += val; }
+
+  /*!
+   * \brief Get information about the total deformation of the design variables.
+   * \param[in] val_dv - design variable number
+   * \param[in] val_val - design variable number
+   * \return Value of total deformation for one design variable.
+   */
+  su2double GetTotalDeformation(unsigned short val_dv, unsigned short val_val = 0) const { return TotalDeformation[val_dv][val_val]; }
+
   /*!
    * \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. */
@@ -1713,5 +1714,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

@@ -484,4 +484,22 @@ class CSurfaceMovement : public CGridMovement {
    * \param[in] config - Definition of the particular problem.
    */
   void SetSurface_Derivative(CGeometry *geometry, CConfig *config);
+
+  /*!
+   * \brief Set magnitude of the deformation of all design variables
+   * \param[in] geometry - Geometrical definition of the problem.
+   * \param[in] config - Definition of the particular problem.
+   * \param[in] deformationFactor - Remaining amount of original deformation (between 0 and 1).
+   * \param[in] totalDeformation_initial - Initial deformation magnitude of all design variables.
+   */
+  void SetDeformationMagnitude(CGeometry *geometry, CConfig *config, su2double deformationFactor, su2double** totalDeformation_initial);
+
+  /*!
+   * \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

@@ -2603,6 +2603,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_Iter, 10);
+
+  /* DESCRIPTION: Recursion depth in the FFD self-intersection prevention */
+  addUnsignedShortOption("FFD_INTPREV_DEPTH", FFD_IntPrev_Depth, 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_Iter, 10);
+
+  /* DESCRIPTION: Recursion depth in the FFD self-intersection prevention */
+  addUnsignedShortOption("CONVEXITY_CHECK_DEPTH", ConvexityCheck_Depth, 3);
+
   /* DESCRIPTION: Definition of the FFD boxes */
   addFFDDefOption("FFD_DEFINITION", nFFDBox, CoordFFDBox, TagFFDBox);