Skip to content

Commit dd4acca

Browse files
pcarruscagpcarruscag
authored
Merge pull request #1142 from su2code/address_review_1076
Address comments from #1076...
2 parents 90d156a + 00c2ec4 commit dd4acca

7 files changed

Lines changed: 112 additions & 161 deletions

File tree

‎Common/include/CConfig.hpp‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -8093,13 +8093,13 @@ class CConfig {
80938093
su2double GetFFD_Tol(void) const { return FFD_Tol; }
80948094

80958095
/*!
8096-
* \brief Get information about whether to do a check on self-intersections within
8096+
* \brief Get information about whether to do a check on self-intersections within
80978097
the FFD box based on value on the Jacobian determinant.
80988098
* \param[out] FFD_IntPrev: <code>TRUE</code> if FFD intersection prevention is active; otherwise <code>FALSE</code>.
80998099
* \param[out] FFD_IntPrev_MaxIter: Maximum number of iterations in the intersection prevention procedure.
81008100
* \param[out] FFD_IntPrev_MaxDepth: Maximum recursion depth in the intersection prevention procedure.
81018101
*/
8102-
tuple<bool, unsigned short, unsigned short> GetFFD_IntPrev(void) {
8102+
tuple<bool, unsigned short, unsigned short> GetFFD_IntPrev(void) const {
81038103
return make_tuple(FFD_IntPrev, FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth);
81048104
}
81058105

@@ -8109,7 +8109,7 @@ class CConfig {
81098109
* \param[out] ConvexityCheck_MaxIter: Maximum number of iterations in the convexity check.
81108110
* \param[out] ConvexityCheck_MaxDepth: Maximum recursion depth in the convexity check.
81118111
*/
8112-
tuple<bool, unsigned short, unsigned short> GetConvexityCheck(void) {
8112+
tuple<bool, unsigned short, unsigned short> GetConvexityCheck(void) const {
81138113
return make_tuple(ConvexityCheck, ConvexityCheck_MaxIter, ConvexityCheck_MaxDepth);
81148114
}
81158115

‎Common/include/geometry/CGeometry.hpp‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1710,6 +1710,6 @@ class CGeometry {
17101710
* \brief Get the amount of nonconvex elements in the mesh.
17111711
* \param[out] nNonconvexElements- amount of nonconvex elements in the mesh
17121712
*/
1713-
unsigned long GetnNonconvexElements(void) {return nNonconvexElements;}
1713+
unsigned long GetnNonconvexElements() const {return nNonconvexElements;}
17141714
};
17151715

‎Common/include/grid_movement/CGridMovement.hpp‎

Lines changed: 5 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -55,10 +55,12 @@ class CGridMovement {
5555
virtual ~CGridMovement(void);
5656

5757
/*!
58-
* \brief A pure virtual member.
58+
* \brief Set the surface/boundary deformation.
5959
* \param[in] geometry - Geometrical definition of the problem.
6060
* \param[in] config - Definition of the particular problem.
61+
* \return Total deformation applied, which may be less than target if intersection prevention is used.
6162
*/
62-
inline virtual void SetSurface_Deformation(CGeometry *geometry, CConfig *config) { }
63-
63+
inline virtual vector<vector<su2double> > SetSurface_Deformation(CGeometry *geometry, CConfig *config) {
64+
return vector<vector<su2double> >();
65+
}
6466
};

‎Common/include/grid_movement/CSurfaceMovement.hpp‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -210,8 +210,9 @@ class CSurfaceMovement : public CGridMovement {
210210
* \brief Set the surface/boundary deformation.
211211
* \param[in] geometry - Geometrical definition of the problem.
212212
* \param[in] config - Definition of the particular problem.
213+
* \return Total deformation applied, which may be less than target if intersection prevention is used.
213214
*/
214-
void SetSurface_Deformation(CGeometry *geometry, CConfig *config, su2double** totaldeformation = nullptr);
215+
vector<vector<su2double> > SetSurface_Deformation(CGeometry *geometry, CConfig *config) override;
215216

216217
/*!
217218
* \brief Compute the parametric coordinates of a grid point using a point inversion strategy
@@ -501,5 +502,5 @@ class CSurfaceMovement : public CGridMovement {
501502
* \param[in] FFDBox - Free form deformation box.
502503
* \return Number of points with negative Jacobian determinant.
503504
*/
504-
unsigned long calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox);
505+
unsigned long calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox) const;
505506
};

‎Common/src/grid_movement/CSurfaceMovement.cpp‎

Lines changed: 68 additions & 76 deletions
Original file line numberDiff line numberDiff line change
@@ -39,19 +39,21 @@ CSurfaceMovement::CSurfaceMovement(void) : CGridMovement() {
3939

4040
CSurfaceMovement::~CSurfaceMovement(void) {}
4141

42-
void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *config, su2double** totaldeformation) {
42+
vector<vector<su2double> > CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *config) {
4343

4444
unsigned short iFFDBox, iDV, iLevel, iChild, iParent, jFFDBox, iMarker;
4545
unsigned short Degree_Unitary [] = {1,1,1}, BSpline_Unitary [] = {2,2,2};
4646
su2double MaxDiff, Current_Scale, Ratio, New_Scale;
4747
string FFDBoxTag;
48-
bool allmoving;
48+
bool allmoving;
4949

50-
bool cylindrical = (config->GetFFD_CoordSystem() == CYLINDRICAL);
51-
bool spherical = (config->GetFFD_CoordSystem() == SPHERICAL);
52-
bool polar = (config->GetFFD_CoordSystem() == POLAR);
53-
bool cartesian = (config->GetFFD_CoordSystem() == CARTESIAN);
54-
su2double BoundLimit = config->GetOpt_LineSearch_Bound();
50+
const bool cylindrical = (config->GetFFD_CoordSystem() == CYLINDRICAL);
51+
const bool spherical = (config->GetFFD_CoordSystem() == SPHERICAL);
52+
const bool polar = (config->GetFFD_CoordSystem() == POLAR);
53+
const bool cartesian = (config->GetFFD_CoordSystem() == CARTESIAN);
54+
const su2double BoundLimit = config->GetOpt_LineSearch_Bound();
55+
56+
vector<vector<su2double> > totaldeformation;
5557

5658
/*--- Setting the Free Form Deformation ---*/
5759

@@ -154,7 +156,7 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
154156
}
155157

156158
else {
157-
SU2_MPI::Error("There are not FFD boxes in the mesh file!!", CURRENT_FUNCTION);
159+
SU2_MPI::Error("There are no FFD boxes in the mesh file!!", CURRENT_FUNCTION);
158160
}
159161

160162
}
@@ -320,33 +322,35 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
320322
}
321323

322324
/*--- Set total deformation values in config ---*/
323-
if (config->GetKind_SU2() == SU2_DEF){
324-
unsigned short iDV_Value;
325-
su2double dv_value;
325+
if (config->GetKind_SU2() == SU2_DEF) {
326+
327+
totaldeformation.resize(config->GetnDV());
326328
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
327-
for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
328-
dv_value = config->GetDV_Value(iDV, iDV_Value);
329-
totaldeformation[iDV][iDV_Value] = dv_value;
329+
totaldeformation[iDV].resize(config->GetnDV_Value(iDV));
330+
for (auto iDV_Value = 0u; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
331+
totaldeformation[iDV][iDV_Value] = config->GetDV_Value(iDV, iDV_Value);
330332
}
331333
}
332334

335+
/*--- Get parameters for FFD self-intersection prevention ---*/
336+
bool FFD_IntPrev;
337+
unsigned short FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth;
338+
339+
tie(FFD_IntPrev, FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth) = config->GetFFD_IntPrev();
340+
333341
/*--- Calculate Jacobian determinant for FFD-deformation to check for self-intersection in FFD box.
334342
Procedure from J.E. Gain & N.A. Dodgson, "Preventing Self-Intersection under Free Form Deformation".
335343
IEEE transactions on Visualization and Computer Graphics, vol. 7 no. 4. October-December 2001 ---*/
336-
unsigned long nNegativeDeterminants, nNegativeDeterminants_previous;
344+
unsigned long nNegativeDeterminants = 0, nNegativeDeterminants_previous;
337345
su2double DeformationDifference = 1.0, DeformationFactor = 1.0;
338346

339-
nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
340-
341-
/*--- Get parameters for FFD self-intersection prevention ---*/
342-
bool FFD_IntPrev;
343-
unsigned short FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth;
344-
345-
tie(FFD_IntPrev, FFD_IntPrev_MaxIter, FFD_IntPrev_MaxDepth) = config->GetFFD_IntPrev();
347+
if (FFD_IntPrev) {
348+
nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
349+
}
346350

347-
/*--- If enabled: start recursive procedure to decrease deformation magnitude
351+
/*--- If enabled: start recursive procedure to decrease deformation magnitude
348352
to remove self-intersections in FFD box ---*/
349-
if (FFD_IntPrev && nNegativeDeterminants > 0) {
353+
if (nNegativeDeterminants > 0) {
350354

351355
if (rank == MASTER_NODE) {
352356
cout << "Self-intersections within FFD box present. ";
@@ -355,8 +359,8 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
355359

356360
/*--- Lower the deformation magnitude and add this to the total deformation value in config ---*/
357361
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
358-
for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
359-
dv_value = config->GetDV_Value(iDV, iDV_Value);
362+
for (auto iDV_Value = 0u; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
363+
auto dv_value = config->GetDV_Value(iDV, iDV_Value);
360364
config->SetDV_Value(iDV, iDV_Value, -dv_value/2);
361365
totaldeformation[iDV][iDV_Value] -= dv_value/2;
362366
}
@@ -379,8 +383,9 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
379383
for (iChild = 0; iChild < FFDBox[iFFDBox]->GetnChildFFDBox(); iChild++) {
380384
FFDBoxTag = FFDBox[iFFDBox]->GetChildFFDBoxTag(iChild);
381385
for (jFFDBox = 0; jFFDBox < GetnFFDBox(); jFFDBox++)
382-
if (FFDBoxTag == FFDBox[jFFDBox]->GetTag()) break;
383-
SetParametricCoordCP(geometry, config, FFDBox[iFFDBox], FFDBox[jFFDBox]);
386+
if (FFDBoxTag == FFDBox[jFFDBox]->GetTag())
387+
break;
388+
SetParametricCoordCP(geometry, config, FFDBox[iFFDBox], FFDBox[jFFDBox]);
384389
}
385390

386391
/*--- Update the parametric coordinates if it is a child FFDBox ---*/
@@ -394,27 +399,27 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
394399
MaxDiff = SetCartesianCoord(geometry, config, FFDBox[iFFDBox], iFFDBox, false);
395400

396401
/*--- Check for self-intersections in FFD box ---*/
397-
nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
402+
nNegativeDeterminants = calculateJacobianDeterminant(geometry, config, FFDBox[iFFDBox]);
398403

399404
if (rank == MASTER_NODE) {
400405
cout << "Amount of points with negative Jacobian determinant for iteration ";
401406
cout << FFD_IntPrev_Iter << ": " << nNegativeDeterminants << endl;
402407
cout << "Remaining amount of original deformation: " << DeformationFactor*100.0 << " percent." << endl;
403408
}
404409

405-
/*--- Recursively change deformation magnitude.
410+
/*--- Recursively change deformation magnitude.
406411
Increase if there are no points with negative determinants, decrease otherwise. ---*/
407412
if (nNegativeDeterminants == 0){
408413
DeformationDifference = abs(DeformationDifference/2.0);
409414

410-
/*--- Update recursion depth if there are no points with negative determinant.
415+
/*--- Update recursion depth if there are no points with negative determinant.
411416
Quit if maximum depth is reached. ---*/
412417
FFD_IntPrev_Depth++;
413418

414419
if (FFD_IntPrev_Depth == FFD_IntPrev_MaxDepth) {
415420
if (rank == MASTER_NODE) {
416421
cout << "Maximum recursion depth reached." << endl;
417-
cout << "Remaining amount of original deformation: " << endl;
422+
cout << "Remaining amount of original deformation: " << endl;
418423
cout << DeformationFactor*100.0 << " percent." << endl;
419424
}
420425
break;
@@ -427,34 +432,28 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
427432
DeformationFactor += DeformationDifference;
428433
}
429434

430-
/*--- Set DV values for next iteration. Always decrease absolute value,
435+
/*--- Set DV values for next iteration. Always decrease absolute value,
431436
since starting point of iteration is previously deformed value. ---*/
432437
if (FFD_IntPrev_Iter < FFD_IntPrev_MaxIter) {
433-
if ((nNegativeDeterminants_previous > 0 && nNegativeDeterminants > 0) ||
438+
su2double sign = -1.0;
439+
if ((nNegativeDeterminants_previous > 0 && nNegativeDeterminants > 0) ||
434440
(nNegativeDeterminants_previous == 0 && nNegativeDeterminants == 0)){
435-
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
436-
for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
437-
dv_value = config->GetDV_Value(iDV, iDV_Value);
438-
config->SetDV_Value(iDV, iDV_Value, dv_value/2.0);
439-
totaldeformation[iDV][iDV_Value] += dv_value/2.0;
440-
}
441-
}
442-
} else {
443-
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
444-
for (iDV_Value = 0; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
445-
dv_value = config->GetDV_Value(iDV, iDV_Value);
446-
config->SetDV_Value(iDV, iDV_Value, -dv_value/2.0);
447-
totaldeformation[iDV][iDV_Value] -= dv_value/2.0;
448-
}
441+
sign = 1.0;
442+
}
443+
for (iDV = 0; iDV < config->GetnDV(); iDV++) {
444+
for (auto iDV_Value = 0u; iDV_Value < config->GetnDV_Value(iDV); iDV_Value++) {
445+
auto dv_value = sign * config->GetDV_Value(iDV, iDV_Value);
446+
config->SetDV_Value(iDV, iDV_Value, dv_value/2.0);
447+
totaldeformation[iDV][iDV_Value] += dv_value/2.0;
449448
}
450449
}
451450
}
452-
453451
nNegativeDeterminants_previous = nNegativeDeterminants;
454452
}
455453

456454
}
457-
}
455+
456+
} // end SU2_DEF
458457

459458
/*--- Reparametrization of the parent FFD box ---*/
460459

@@ -679,6 +678,7 @@ void CSurfaceMovement::SetSurface_Deformation(CGeometry *geometry, CConfig *conf
679678
cout << "Design Variable not implemented yet" << endl;
680679
}
681680

681+
return totaldeformation;
682682
}
683683

684684

@@ -694,7 +694,7 @@ void CSurfaceMovement::SetSurface_Derivative(CGeometry *geometry, CConfig *confi
694694
DV_Value = config->GetDV_Value(iDV, iDV_Value);
695695

696696
/*--- If value of the design variable is not 0.0 we apply the differentation.
697-
* Note if multiple variables are non-zero, we end up with the sum of all the derivatives. ---*/
697+
* Note if multiple variables are non-zero, we end up with the sum of all the derivatives. ---*/
698698

699699
if (DV_Value != 0.0) {
700700

@@ -1467,7 +1467,7 @@ void CSurfaceMovement::ApplyDesignVariables(CGeometry *geometry, CConfig *config
14671467
case FFD_THICKNESS : SetFFDThickness(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
14681468
case FFD_ANGLE_OF_ATTACK : SetFFDAngleOfAttack(geometry, config, FFDBox[iFFDBox], FFDBox, iDV, false); break;
14691469
}
1470-
}
1470+
}
14711471
}
14721472

14731473
su2double CSurfaceMovement::SetCartesianCoord(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox, unsigned short iFFDBox, bool ResetDef) {
@@ -4913,9 +4913,8 @@ void CSurfaceMovement::WriteFFDInfo(CSurfaceMovement** surface_movement, CGeomet
49134913
}
49144914
}
49154915

4916-
unsigned long CSurfaceMovement::calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox){
4916+
unsigned long CSurfaceMovement::calculateJacobianDeterminant(CGeometry *geometry, CConfig *config, CFreeFormDefBox *FFDBox) const {
49174917

4918-
su2double *ParamCoord;
49194918
unsigned long iSurfacePoints;
49204919
unsigned short iMarker;
49214920
unsigned long negative_determinants = 0;
@@ -4925,41 +4924,34 @@ unsigned long CSurfaceMovement::calculateJacobianDeterminant(CGeometry *geometry
49254924

49264925
/*--- Get the marker of the surface point ---*/
49274926
iMarker = FFDBox->Get_MarkerIndex(iSurfacePoints);
4928-
4927+
49294928
if (config->GetMarker_All_DV(iMarker) == YES) {
49304929

4931-
ParamCoord = FFDBox->Get_ParametricCoord(iSurfacePoints);
4930+
const auto ParamCoord = FFDBox->Get_ParametricCoord(iSurfacePoints);
49324931

49334932
/*--- Calculate partial derivatives ---*/
4934-
unsigned short iDegree, jDegree, kDegree, lDegree, mDegree, nDegree;
4933+
unsigned short iDegree, jDegree, kDegree;
49354934
su2double Ba, Bb, Bc, Ba_der, Bb_der, Bc_der;
4936-
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};
4935+
su2double determinant, d_du[3] = {0.0}, d_dv[3] = {0.0}, d_dw[3] = {0.0};
49374936

4938-
for (iDegree = 0; iDegree <= FFDBox->lDegree; iDegree++){
4937+
for (iDegree = 0; iDegree <= FFDBox->lDegree; iDegree++){
49394938
Ba = FFDBox->BlendingFunction[0]->GetBasis(iDegree, ParamCoord[0]);
49404939
Ba_der = FFDBox->BlendingFunction[0]->GetDerivative(iDegree, ParamCoord[0], 1);
49414940

4942-
for (jDegree = 0; jDegree <= FFDBox->mDegree; jDegree++){
4941+
for (jDegree = 0; jDegree <= FFDBox->mDegree; jDegree++){
49434942
Bb = FFDBox->BlendingFunction[1]->GetBasis(jDegree, ParamCoord[1]);
49444943
Bb_der = FFDBox->BlendingFunction[1]->GetDerivative(jDegree, ParamCoord[1], 1);
49454944

4946-
for (kDegree = 0; kDegree <= FFDBox->nDegree; kDegree++){
4947-
4945+
for (kDegree = 0; kDegree <= FFDBox->nDegree; kDegree++){
4946+
49484947
Bc = FFDBox->BlendingFunction[2]->GetBasis(kDegree, ParamCoord[2]);
49494948
Bc_der = FFDBox->BlendingFunction[2]->GetDerivative(kDegree, ParamCoord[2], 1);
49504949

4951-
d_du[0] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
4952-
d_du[1] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
4953-
d_du[2] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
4954-
4955-
d_dv[0] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
4956-
d_dv[1] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
4957-
d_dv[2] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
4958-
4959-
d_dw[0] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][0];
4960-
d_dw[1] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][1];
4961-
d_dw[2] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][2];
4962-
4950+
for (int i=0; i<3; ++i) {
4951+
d_du[i] += Ba_der*Bb*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][i];
4952+
d_dv[i] += Ba*Bb_der*Bc*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][i];
4953+
d_dw[i] += Ba*Bb*Bc_der*FFDBox->Coord_Control_Points[iDegree][jDegree][kDegree][i];
4954+
}
49634955
}
49644956
}
49654957
}
@@ -4973,8 +4965,8 @@ unsigned long CSurfaceMovement::calculateJacobianDeterminant(CGeometry *geometry
49734965
}
49744966
}
49754967

4976-
unsigned long negative_determinants_local = negative_determinants; negative_determinants = 0;
4977-
SU2_MPI::Allreduce(&negative_determinants_local, &negative_determinants, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
4968+
unsigned long tmp = negative_determinants;
4969+
SU2_MPI::Allreduce(&tmp, &negative_determinants, 1, MPI_UNSIGNED_LONG, MPI_SUM, MPI_COMM_WORLD);
49784970

49794971
return negative_determinants;
49804972
}

0 commit comments

Comments
 (0)