diff --git a/Common/include/CConfig.hpp b/Common/include/CConfig.hpp
index 493feb41d3fd..f7ebd7b4ec33 100644
--- a/Common/include/CConfig.hpp
+++ b/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: TRUE if FFD intersection prevention is active; otherwise FALSE.
+ * \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 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: TRUE if convexity check is active; otherwise FALSE.
+ * \param[out] ConvexityCheck_MaxIter: Maximum number of iterations in the convexity check.
+ * \param[out] ConvexityCheck_MaxDepth: Maximum recursion depth in the convexity check.
+ */
+ tuple 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.
diff --git a/Common/include/geometry/CGeometry.hpp b/Common/include/geometry/CGeometry.hpp
index 583c20b963c9..ab3797ac4db4 100644
--- a/Common/include/geometry/CGeometry.hpp
+++ b/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;}
};
diff --git a/Common/include/grid_movement/CSurfaceMovement.hpp b/Common/include/grid_movement/CSurfaceMovement.hpp
index ac901f0379e6..8e0ce8e35878 100644
--- a/Common/include/grid_movement/CSurfaceMovement.hpp
+++ b/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);
};
diff --git a/Common/include/grid_movement/CVolumetricMovement.hpp b/Common/include/grid_movement/CVolumetricMovement.hpp
index fd9874991e2e..79293aabc109 100644
--- a/Common/include/grid_movement/CVolumetricMovement.hpp
+++ b/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.
diff --git a/Common/src/CConfig.cpp b/Common/src/CConfig.cpp
index fab0df9e5dca..c290363577c4 100644
--- a/Common/src/CConfig.cpp
+++ b/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);
diff --git a/Common/src/grid_movement/CSurfaceMovement.cpp b/Common/src/grid_movement/CSurfaceMovement.cpp
index 0ab70a39d1be..1702c71b929c 100644
--- a/Common/src/grid_movement/CSurfaceMovement.cpp
+++ b/Common/src/grid_movement/CSurfaceMovement.cpp
@@ -39,7 +39,7 @@ CSurfaceMovement::CSurfaceMovement(void) : CGridMovement() {
CSurfaceMovement::~CSurfaceMovement(void) {}
-void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *config) {
+void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *config, su2double** totaldeformation) {
unsigned short iFFDBox, iDV, iLevel, iChild, iParent, jFFDBox, iMarker;
unsigned short Degree_Unitary [] = {1,1,1}, BSpline_Unitary [] = {2,2,2};
@@ -295,24 +295,7 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
if (iLevel > 0) UpdateParametricCoord(geometry, config, FFDBox[iFFDBox], iFFDBox);
/*--- Apply the design variables to the control point position ---*/
-
- for (iDV = 0; iDV < config->GetnDV(); iDV++) {
- switch ( config->GetDesign_Variable(iDV) ) {
- case FFD_CONTROL_POINT_2D : SetFFDCPChange_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CAMBER_2D : SetFFDCamber_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_THICKNESS_2D : SetFFDThickness_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_TWIST_2D : SetFFDTwist_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CONTROL_POINT : SetFFDCPChange(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_NACELLE : SetFFDNacelle(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_GULL : SetFFDGull(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_TWIST : SetFFDTwist(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_ROTATION : SetFFDRotation(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CONTROL_SURFACE : SetFFDControl_Surface(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CAMBER : SetFFDCamber(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_THICKNESS : SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_ANGLE_OF_ATTACK : SetFFDAngleOfAttack(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- }
- }
+ ApplyDesignVariables(geometry, config, FFDBox, iFFDBox);
/*--- Recompute cartesian coordinates using the new control point location ---*/
@@ -328,29 +311,149 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
config->SetOpt_RelaxFactor(New_Scale);
/*--- Apply the design variables to the control point position ---*/
+ ApplyDesignVariables(geometry, config, FFDBox, iFFDBox);
+
+ /*--- Recompute cartesian coordinates using the new control point location ---*/
+ MaxDiff = SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, false);
+
+ }
+
+ /*--- Set total deformation values in config ---*/
+ if (config->GetKind_SU2() == SU2_DEF){
+ unsigned short iDV_Value;
+ su2double dv_value;
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
- switch ( config->GetDesign_Variable(iDV) ) {
- case FFD_CONTROL_POINT_2D : SetFFDCPChange_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CAMBER_2D : SetFFDCamber_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_THICKNESS_2D : SetFFDThickness_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_TWIST_2D : SetFFDTwist_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CONTROL_POINT : SetFFDCPChange(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_NACELLE : SetFFDNacelle(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_GULL : SetFFDGull(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_TWIST : SetFFDTwist(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_ROTATION : SetFFDRotation(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CONTROL_SURFACE : SetFFDControl_Surface(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_CAMBER : SetFFDCamber(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_THICKNESS : SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
- case FFD_ANGLE_OF_ATTACK : SetFFDAngleOfAttack(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
+ dv_value = config->GetDV_Value(iDV, iDV_Value);
+ totaldeformation[iDV][iDV_Value] = dv_value;
}
}
- /*--- Recompute cartesian coordinates using the new control point location ---*/
+ /*--- Calculate Jacobian determinant for FFD-deformation to check for self-intersection in FFD box.
+ Procedure from J.E. Gain & N.A. Dodgson, "Preventing Self-Intersection under Free Form Deformation".
+ IEEE transactions on Visualization and Computer Graphics, vol. 7 no. 4. October-December 2001 ---*/
+ unsigned long nNegativeDeterminants, nNegativeDeterminants_previous;
+ su2double DeformationDifference = 1.0, DeformationFactor = 1.0;
- MaxDiff = SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, false);
+ nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
+
+ /*--- Get parameters for FFD self-intersection prevention ---*/
+ bool FFD_IntPrev;
+ unsigned short FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth;
+
+ tie(FFD_IntPrev, FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth) = config->GetFFD_IntPrev();
+
+ /*--- If enabled: start recursive procedure to decrease deformation magnitude
+ to remove self-intersections in FFD box ---*/
+ if (FFD_IntPrev && nNegativeDeterminants > 0) {
+
+ if (rank == MASTER_NODE) {
+ cout << "Self-intersections within FFD box present. ";
+ cout << "Performing iterative deformation reduction procedure." << endl;
+ }
+
+ /*--- Lower the deformation magnitude and add this to the total deformation value in config ---*/
+ for (iDV = 0; iDV < config->GetnDV(); iDV++) {
+ for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
+ dv_value = config->GetDV_Value(iDV, iDV_Value);
+ config->SetDV_Value(iDV, iDV_Value, -dv_value/2);
+ totaldeformation[iDV][iDV_Value] -= dv_value/2;
+ }
+ }
+
+ DeformationDifference /= 2.0;
+ DeformationFactor -= DeformationDifference;
+ nNegativeDeterminants_previous = nNegativeDeterminants;
+
+ /*--- Recursively check for self-intersections. ---*/
+ unsigned short FFD_IntPrev_Iter, FFD_IntPrev_Depth = 0;
+ for (FFD_IntPrev_Iter = 1; FFD_IntPrev_Iter <= FFD_IntPrev_MaxIter; FFD_IntPrev_Iter++){
+
+ if (rank == MASTER_NODE) cout << "Checking FFD box intersections with the solid surfaces." << endl;
+ CheckFFDIntersections(geometry, config, FFDBox[iFFDBox], iFFDBox);
+
+ /*--- Compute the parametric coordinates of the child box
+ control points (using the parent FFDBox) ---*/
+ for (iChild = 0; iChild < FFDBox[iFFDBox]->GetnChildFFDBox(); iChild++) {
+ FFDBoxTag = FFDBox[iFFDBox]->GetChildFFDBoxTag(iChild);
+ for (jFFDBox = 0; jFFDBox < GetnFFDBox(); jFFDBox++)
+ if (FFDBoxTag == FFDBox[jFFDBox]->GetTag()) break;
+ SetParametricCoordCP(geometry, config, FFDBox[iFFDBox], FFDBox[jFFDBox]);
+ }
+
+ /*--- Update the parametric coordinates if it is a child FFDBox ---*/
+ iLevel = FFDBox[iFFDBox]->GetLevel();
+ if (iLevel > 0) UpdateParametricCoord(geometry, config, FFDBox[iFFDBox], iFFDBox);
+
+ /*--- Apply the design variables to the control point position ---*/
+ ApplyDesignVariables(geometry, config, FFDBox, iFFDBox);
+
+ /*--- Recompute cartesian coordinates using the new control point location ---*/
+ MaxDiff = SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, false);
+
+ /*--- Check for self-intersections in FFD box ---*/
+ nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
+
+ if (rank == MASTER_NODE) {
+ cout << "Amount of points with negative Jacobian determinant for iteration ";
+ cout << FFD_IntPrev_Iter << ": " << nNegativeDeterminants << endl;
+ cout << "Remaining amount of original deformation: " << DeformationFactor*100.0 << " percent." << endl;
+ }
+
+ /*--- Recursively change deformation magnitude.
+ Increase if there are no points with negative determinants, decrease otherwise. ---*/
+ if (nNegativeDeterminants == 0){
+ DeformationDifference = abs(DeformationDifference/2.0);
+
+ /*--- Update recursion depth if there are no points with negative determinant.
+ Quit if maximum depth is reached. ---*/
+ FFD_IntPrev_Depth++;
+
+ if (FFD_IntPrev_Depth == FFD_IntPrev_MaxDepth) {
+ if (rank == MASTER_NODE) {
+ cout << "Maximum recursion depth reached." << endl;
+ cout << "Remaining amount of original deformation: " << endl;
+ cout << DeformationFactor*100.0 << " percent." << endl;
+ }
+ break;
+ }
+ } else {
+ DeformationDifference = -abs(DeformationDifference/2.0);
+ }
+
+ if (FFD_IntPrev_Iter < FFD_IntPrev_MaxIter) {
+ DeformationFactor += DeformationDifference;
+ }
+
+ /*--- Set DV values for next iteration. Always decrease absolute value,
+ since starting point of iteration is previously deformed value. ---*/
+ if (FFD_IntPrev_Iter < FFD_IntPrev_MaxIter) {
+ if ((nNegativeDeterminants_previous > 0 && nNegativeDeterminants > 0) ||
+ (nNegativeDeterminants_previous == 0 && nNegativeDeterminants == 0)){
+ for (iDV = 0; iDV < config->GetnDV(); iDV++) {
+ for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
+ dv_value = config->GetDV_Value(iDV, iDV_Value);
+ config->SetDV_Value(iDV, iDV_Value, dv_value/2.0);
+ totaldeformation[iDV][iDV_Value] += dv_value/2.0;
+ }
+ }
+ } else {
+ for (iDV = 0; iDV < config->GetnDV(); iDV++) {
+ for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
+ dv_value = config->GetDV_Value(iDV, iDV_Value);
+ config->SetDV_Value(iDV, iDV_Value, -dv_value/2.0);
+ totaldeformation[iDV][iDV_Value] -= dv_value/2.0;
+ }
+ }
+ }
+ }
+
+ nNegativeDeterminants_previous = nNegativeDeterminants;
+ }
+
+ }
}
/*--- Reparametrization of the parent FFD box ---*/
@@ -1344,6 +1447,29 @@ void CSurfaceMovement::UpdateParametricCoord(CGeometry *geometry, CConfig *confi
}
+void CSurfaceMovement::ApplyDesignVariables(CGeometry *geometry, CConfig *config, CFreeFormDefBox **FFDBox, unsigned short iFFDBox) {
+
+ unsigned short iDV;
+
+ for (iDV = 0; iDV < config->GetnDV(); iDV++) {
+ switch ( config->GetDesign_Variable(iDV) ) {
+ case FFD_CONTROL_POINT_2D : SetFFDCPChange_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_CAMBER_2D : SetFFDCamber_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_THICKNESS_2D : SetFFDThickness_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_TWIST_2D : SetFFDTwist_2D(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_CONTROL_POINT : SetFFDCPChange(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_NACELLE : SetFFDNacelle(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_GULL : SetFFDGull(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_TWIST : SetFFDTwist(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_ROTATION : SetFFDRotation(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_CONTROL_SURFACE : SetFFDControl_Surface(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_CAMBER : SetFFDCamber(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_THICKNESS : SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ case FFD_ANGLE_OF_ATTACK : SetFFDAngleOfAttack(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
+ }
+ }
+}
+
su2double CSurfaceMovement::SetCartesianCoord(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox, unsigned short iFFDBox, bool ResetDef) {
su2double *CartCoordNew, Diff, my_MaxDiff = 0.0, MaxDiff,
@@ -4785,5 +4911,70 @@ void CSurfaceMovement::WriteFFDInfo(CSurfaceMovement** surface_movement, CGeomet
output_file.close();
}
+}
+
+unsigned long CSurfaceMovement::calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox){
+
+ su2double *ParamCoord;
+ unsigned long iSurfacePoints;
+ unsigned short iMarker;
+ unsigned long negative_determinants = 0;
+
+ /*--- Loop over the surface points ---*/
+ for (iSurfacePoints = 0; iSurfacePoints < FFDBox->GetnSurfacePoint(); iSurfacePoints++) {
+
+ /*--- Get the marker of the surface point ---*/
+ iMarker = FFDBox->Get_MarkerIndex(iSurfacePoints);
+
+ if (config->GetMarker_All_DV(iMarker) == YES) {
+
+ ParamCoord = FFDBox->Get_ParametricCoord(iSurfacePoints);
+
+ /*--- Calculate partial derivatives ---*/
+ unsigned short iDegree, jDegree, kDegree, lDegree, mDegree, nDegree;
+ su2double Ba, Bb, Bc, Ba_der, Bb_der, Bc_der;
+ su2double determinant, d_du[3] = {0.0, 0.0, 0.0}, d_dv[3] = {0.0, 0.0, 0.0}, d_dw[3] = {0.0, 0.0, 0.0};
+
+ for (iDegree = 0; iDegree <= FFDBox->lDegree; iDegree++){
+ Ba = FFDBox->BlendingFunction[0]->GetBasis(iDegree, ParamCoord[0]);
+ Ba_der = FFDBox->BlendingFunction[0]->GetDerivative(iDegree, ParamCoord[0], 1);
+
+ for (jDegree = 0; jDegree <= FFDBox->mDegree; jDegree++){
+ Bb = FFDBox->BlendingFunction[1]->GetBasis(jDegree, ParamCoord[1]);
+ Bb_der = FFDBox->BlendingFunction[1]->GetDerivative(jDegree, ParamCoord[1], 1);
+
+ for (kDegree = 0; kDegree <= FFDBox->nDegree; kDegree++){
+
+ Bc = FFDBox->BlendingFunction[2]->GetBasis(kDegree, ParamCoord[2]);
+ Bc_der = FFDBox->BlendingFunction[2]->GetDerivative(kDegree, ParamCoord[2], 1);
+
+ d_du[0] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
+ d_du[1] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
+ d_du[2] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
+
+ d_dv[0] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
+ d_dv[1] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
+ d_dv[2] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
+
+ d_dw[0] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
+ d_dw[1] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
+ d_dw[2] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
+
+ }
+ }
+ }
+
+ /*--- Calculate determinant ---*/
+ determinant = d_du[0]*(d_dv[1]*d_dw[2] - d_dv[2]*d_dw[1]) - d_dv[0]*(d_du[1]*d_dw[2] - d_du[2]*d_dw[1]) + d_dw[0]*(d_du[1]*d_dv[2] - d_du[2]*d_dv[1]);
+
+ if (determinant < 0){
+ negative_determinants++;
+ }
+ }
+ }
+
+ unsigned long negative_determinants_local = negative_determinants; negative_determinants = 0;
+ SU2_MPI::Allreduce(&negative_determinants_local, &negative_determinants, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
+ return negative_determinants;
}
diff --git a/Common/src/grid_movement/CVolumetricMovement.cpp b/Common/src/grid_movement/CVolumetricMovement.cpp
index c4cefce292b7..3ee6914f460b 100644
--- a/Common/src/grid_movement/CVolumetricMovement.cpp
+++ b/Common/src/grid_movement/CVolumetricMovement.cpp
@@ -28,6 +28,7 @@
#include "../../include/grid_movement/CVolumetricMovement.hpp"
#include "../../include/adt/CADTPointsOnlyClass.hpp"
+#include "../../include/toolboxes/geometry_toolbox.hpp"
CVolumetricMovement::CVolumetricMovement(void) : CGridMovement(), System(true) {
@@ -200,6 +201,9 @@ void CVolumetricMovement::SetVolume_Deformation(CGeometry *geometry, CConfig *co
ComputeDeforming_Element_Volume(geometry, MinVolume, MaxVolume, Screen_Output);
+ /*--- Calculate amount of nonconvex elements ---*/
+ ComputenNonconvexElements(geometry, Screen_Output);
+
/*--- Set number of iterations in the mesh update. ---*/
Set_nIterMesh(Tot_Iter);
@@ -293,6 +297,135 @@ void CVolumetricMovement::ComputeDeforming_Element_Volume(CGeometry *geometry, s
}
+void CVolumetricMovement::ComputenNonconvexElements(CGeometry *geometry, bool Screen_Output) {
+ unsigned long iElem;
+ unsigned short iDim, nNodes;
+ unsigned long nNonconvexElements = 0;
+
+ /*--- Load up each tetrahedron to check for convex properties. ---*/
+ if (nDim == 2){
+ for (iElem = 0; iElem < geometry->GetnElem(); iElem++) {
+ su2double minCrossProduct = 1.e6, maxCrossProduct = -1.e6;
+
+ if (geometry->elem[iElem]->GetVTK_Type() == TRIANGLE) nNodes = 3;
+ if (geometry->elem[iElem]->GetVTK_Type() == QUADRILATERAL) nNodes = 4;
+
+ /*--- Get coordinates of corner points ---*/
+ unsigned short iNodes;
+ unsigned long PointCorners[8];
+ su2double CoordCorners[8][3];
+
+ for (iNodes = 0; iNodes < nNodes; iNodes++) {
+ PointCorners[iNodes] = geometry->elem[iElem]->GetNode(iNodes);
+ for (iDim = 0; iDim < nDim; iDim++) {
+ CoordCorners[iNodes][iDim] = geometry->nodes->GetCoord(PointCorners[iNodes], iDim);
+ }
+ }
+
+ /*--- Determine whether element is convex ---*/
+ for (iNodes = 0; iNodes < nNodes; iNodes ++) {
+
+ /*--- Calculate minimum and maximum angle between edge vectors adjacent to each node ---*/
+ su2double edgeVector_i[nDim], edgeVector_j[nDim];
+
+ for (iDim = 0; iDim < nDim; iDim ++) {
+ if (iNodes == 0) {
+ edgeVector_i[iDim] = CoordCorners[nNodes-1][iDim] - CoordCorners[iNodes][iDim];
+ } else {
+ edgeVector_i[iDim] = CoordCorners[iNodes-1][iDim] - CoordCorners[iNodes][iDim];
+ }
+
+ if (iNodes == nNodes-1) {
+ edgeVector_j[iDim] = CoordCorners[0][iDim] - CoordCorners[iNodes][iDim];
+ } else {
+ edgeVector_j[iDim] = CoordCorners[iNodes+1][iDim] - CoordCorners[iNodes][iDim];
+ }
+ }
+
+ /*--- Calculate cross product of edge vectors ---*/
+ su2double crossProduct;
+ crossProduct = edgeVector_i[1]*edgeVector_j[0] - edgeVector_i[0]*edgeVector_j[1];
+
+ if (crossProduct < minCrossProduct) minCrossProduct = crossProduct;
+ if (crossProduct > maxCrossProduct) maxCrossProduct = crossProduct;
+ }
+
+ /*--- Element is nonconvex if cross product of at least one set of adjacent edges is negative ---*/
+ if (minCrossProduct < 0 && maxCrossProduct > 0){
+ nNonconvexElements++;
+ }
+ }
+ } else {
+
+ /*--- 3D elements ---*/
+ unsigned short iNode, iFace, nFaceNodes;
+
+ for (iElem = 0; iElem < geometry->GetnElem(); iElem++){
+ for (iFace = 0; iFace < geometry->elem[iElem]->GetnFaces(); iFace++){
+ nFaceNodes = geometry->elem[iElem]->GetnNodesFace(iFace);
+
+ su2double minCrossProductLength = 1.e6, maxCrossProductLength = -1.e6;
+ for (iNode = 0; iNode < nFaceNodes; iNode++) {
+ su2double crossProductLength = 0.0;
+
+ /*--- Get coords of node face_point_i and its adjacent nodes in the face ---*/
+ unsigned long face_point_i, face_point_j, face_point_k;
+
+ face_point_i = geometry->elem[iElem]->GetNode(geometry->elem[iElem]->GetFaces(iFace, iNode));
+
+ if (iNode == 0) {
+ face_point_j = geometry->elem[iElem]->GetNode(geometry->elem[iElem]->GetFaces(iFace, nFaceNodes-1));
+ face_point_k = geometry->elem[iElem]->GetNode(geometry->elem[iElem]->GetFaces(iFace, iNode+1));
+ } else if (iNode == nFaceNodes-1) {
+ face_point_j = geometry->elem[iElem]->GetNode(geometry->elem[iElem]->GetFaces(iFace, iNode-1));
+ face_point_k = geometry->elem[iElem]->GetNode(geometry->elem[iElem]->GetFaces(iFace, 0));
+ }
+
+ su2double *Coords_i, *Coords_j, *Coords_k;
+
+ Coords_i = geometry->nodes->GetCoord(face_point_i);
+ Coords_j = geometry->nodes->GetCoord(face_point_j);
+ Coords_k = geometry->nodes->GetCoord(face_point_k);
+
+ /*--- Get edge vectors from point k to i and point j to i ---*/
+ su2double edgeVector_i[nDim], edgeVector_j[nDim], normalVector[nDim];
+
+ for (iDim = 0; iDim < nDim; iDim++) {
+ edgeVector_i[iDim] = Coords_k[iDim] - Coords_i[iDim];
+ edgeVector_j[iDim] = Coords_j[iDim] - Coords_i[iDim];
+ }
+
+ /*--- Calculate cross product of edge vectors and its length---*/
+ su2double crossProduct[nDim];
+ GeometryToolbox::CrossProduct(edgeVector_i, edgeVector_j, crossProduct);
+ crossProductLength = GeometryToolbox::Norm(nDim, crossProduct);
+
+ /*--- Check if length is minimum or maximum ---*/
+ if (crossProductLength < minCrossProductLength) {
+ minCrossProductLength = crossProductLength;
+ } else if (crossProductLength > maxCrossProductLength) {
+ maxCrossProductLength = crossProductLength;
+ }
+ }
+
+ /*--- If minimum cross product length is smaller than 0,
+ face (and therefore element) is not convex ---*/
+ if (minCrossProductLength < 0) {
+ nNonconvexElements++;
+ break;
+ }
+
+ }
+ }
+ }
+
+ unsigned long nNonconvexElements_Local = nNonconvexElements; nNonconvexElements = 0;
+ SU2_MPI::Allreduce(&nNonconvexElements_Local, &nNonconvexElements, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
+
+ /*--- Set number of nonconvex elements in geometry ---*/
+ geometry->SetnNonconvexElements(nNonconvexElements);
+}
+
void CVolumetricMovement::ComputeSolid_Wall_Distance(CGeometry *geometry, CConfig *config, su2double &MinDistance, su2double &MaxDistance) const {
diff --git a/SU2_DEF/src/SU2_DEF.cpp b/SU2_DEF/src/SU2_DEF.cpp
index 96f967c692c2..1a2837100b4a 100644
--- a/SU2_DEF/src/SU2_DEF.cpp
+++ b/SU2_DEF/src/SU2_DEF.cpp
@@ -230,11 +230,34 @@ int main(int argc, char *argv[]) {
for (iZone = 0; iZone < nZone; iZone++){
+ /*--- Initialize total deformation of design variables to zero ---*/
+ su2double **TotalDeformation;
+ unsigned short iDV, iDV_Value;
+ TotalDeformation = new su2double*[config_container[iZone]->GetnDV()];
+
+ for (iDV = 0; iDV < config_container[iZone]->GetnDV(); iDV++){
+ TotalDeformation[iDV] = new su2double[config_container[iZone]->GetnDV_Value(iDV)];
+ for (iDV_Value = 0; iDV_Value < config_container[iZone]->GetnDV_Value(iDV); iDV_Value++){
+ TotalDeformation[iDV][iDV_Value] = 0.0;
+ }
+ }
+
if (config_container[iZone]->GetDesign_Variable(0) != NO_DEFORMATION) {
/*--- Definition of the Class for grid movement ---*/
grid_movement[iZone] = new CVolumetricMovement(geometry_container[iZone], config_container[iZone]);
+ /*--- Save original coordinates to be reused in convexity checking procedure ---*/
+ unsigned long iPoint, nPoint = geometry_container[iZone]->GetnPoint();
+ unsigned short iDim, nDim = geometry_container[iZone]->GetnDim();
+ su2double OriginalCoordinates[nPoint][nDim];
+
+ for (iPoint = 0; iPoint < nPoint; iPoint++) {
+ for (iDim = 0; iDim < nDim; iDim++) {
+ OriginalCoordinates[iPoint][iDim] = geometry_container[iZone]->nodes->GetCoord(iPoint, iDim);
+ }
+ }
+
/*--- First check for volumetric grid deformation/transformations ---*/
if (config_container[iZone]->GetDesign_Variable(0) == SCALE_GRID) {
@@ -274,7 +297,7 @@ int main(int argc, char *argv[]) {
/*--- Surface grid deformation ---*/
if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl;
- surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone], config_container[iZone]);
+ surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone], config_container[iZone], TotalDeformation);
if (config_container[iZone]->GetDesign_Variable(0) != FFD_SETTING) {
@@ -285,10 +308,107 @@ int main(int argc, char *argv[]) {
cout << "Performing the deformation of the volumetric grid." << endl;
grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone], config_container[iZone], false);
+ /*--- Get parameters for convexity check ---*/
+ bool ConvexityCheck;
+ unsigned short ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth;
+
+ tie(ConvexityCheck, ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth) = config_container[iZone]->GetConvexityCheck();
+
+ /*--- Recursively change deformations if there are nonconvex elements. ---*/
+
+ if (ConvexityCheck && geometry_container[iZone]->GetnNonconvexElements() > 0) {
+ if (rank == MASTER_NODE) {
+ cout << "Nonconvex elements present after deformation. " << endl;
+ cout << "Recursively lowering deformation magnitude." << endl;
+ }
+
+ /*--- Load initial deformation values ---*/
+ su2double** InitialDeformation;
+
+ InitialDeformation = new su2double*[config_container[iZone]->GetnDV()];
+ for (iDV = 0; iDV < config_container[iZone]->GetnDV(); iDV++) {
+ InitialDeformation[iDV] = new su2double[config_container[iZone]->GetnDV_Value(iDV)];
+ for (iDV_Value = 0; iDV_Value < config_container[iZone]->GetnDV_Value(iDV); iDV_Value++) {
+ InitialDeformation[iDV][iDV_Value] = TotalDeformation[iDV][iDV_Value];
+ }
+ }
+
+ unsigned short ConvexityCheckIter, RecursionDepth = 0;
+ su2double DeformationFactor = 1.0, DeformationDifference = 1.0;
+ for (ConvexityCheckIter = 1; ConvexityCheckIter <= ConvexityCheck_MaxIter; ConvexityCheckIter++) {
+
+ /*--- Recursively change deformation magnitude:
+ decrease if there are nonconvex elements, increase otherwise ---*/
+ DeformationDifference /= 2.0;
+
+ if (geometry_container[iZone]->GetnNonconvexElements() > 0) {
+ DeformationFactor -= DeformationDifference;
+ } else {
+ RecursionDepth += 1;
+
+ if (RecursionDepth == ConvexityCheck_MaxDepth) {
+ if (rank == MASTER_NODE) {
+ cout << "Maximum recursion depth reached." << endl;
+ cout << "Remaining amount of original deformation: ";
+ cout << DeformationFactor*100.0 << " percent. " << endl;
+ }
+ break;
+ }
+
+ DeformationFactor += DeformationDifference;
+ }
+
+ /*--- Load mesh to start every iteration with an undeformed grid ---*/
+ for (iPoint = 0; iPoint < nPoint; iPoint++) {
+ for (iDim = 0; iDim < nDim; iDim++) {
+ geometry_container[iZone]->nodes->SetCoord(iPoint, iDim, OriginalCoordinates[iPoint][iDim]);
+ }
+ }
+
+ /*--- Set deformation magnitude as percentage of initial deformation ---*/
+ for (iDV = 0; iDV < config->GetnDV(); iDV++) {
+ for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
+ config_container[iZone]->SetDV_Value(iDV, iDV_Value, InitialDeformation[iDV][iDV_Value]*DeformationFactor);
+ }
+ }
+
+ /*--- Surface grid deformation ---*/
+ if (rank == MASTER_NODE) cout << "Performing the deformation of the surface grid." << endl;
+
+ surface_movement[iZone]->SetSurface_Deformation(geometry_container[iZone], config_container[iZone], TotalDeformation);
+
+ if (rank == MASTER_NODE)
+ cout << endl << "------------------- Volumetric grid deformation (ZONE " << iZone <<") ----------------" << endl;
+
+ if (rank == MASTER_NODE)
+ cout << "Performing the deformation of the volumetric grid." << endl;
+ grid_movement[iZone]->SetVolume_Deformation(geometry_container[iZone], config_container[iZone], false);
+
+ if (rank == MASTER_NODE) {
+ cout << "Number of nonconvex elements for iteration " << ConvexityCheckIter << ": ";
+ cout << geometry_container[iZone]->GetnNonconvexElements() << endl;
+ cout << "Remaining amount of original deformation: ";
+ cout << DeformationFactor*100.0 << " percent. " << endl;
+ }
+
+ }
+
+ for (iDV = 0; iDV < config_container[iZone]->GetnDV(); iDV++) {
+ delete InitialDeformation[iDV];
+ }
+ delete [] InitialDeformation;
+
+ }
+
}
}
+ for (iDV = 0; iDV < config_container[iZone]->GetnDV(); iDV++){
+ delete TotalDeformation[iDV];
+ }
+ delete [] TotalDeformation;
+
}
}
diff --git a/config_template.cfg b/config_template.cfg
index 6958a53ecfdb..2a05e54e25f2 100644
--- a/config_template.cfg
+++ b/config_template.cfg
@@ -1352,6 +1352,17 @@ FFD_TOLERANCE= 1E-10
%
% Maximum number of iterations in the Free-Form Deformation point inversion
FFD_ITERATIONS= 500
+
+% Parameters for prevention of self-intersections within FFD box
+FFD_INTPREV = YES
+FFD_INTPREV_ITER = 10
+FFD_INTPREV_DEPTH = 3
+
+% Parameters for prevention of nonconvex elements in mesh after deformation
+CONVEXITY_CHECK = YES
+CONVEXITY_CHECK_ITER = 10
+CONVEXITY_CHECK_DEPTH = 3
+
%
% FFD box definition: 3D case (FFD_BoxTag, X1, Y1, Z1, X2, Y2, Z2, X3, Y3, Z3, X4, Y4, Z4,
% X5, Y5, Z5, X6, Y6, Z6, X7, Y7, Z7, X8, Y8, Z8)