Added prevention of intersections in grid elements after deformation by tollennaert · Pull Request #1076 · su2code/SU2

tollennaert · 2020-09-20T16:20:33Z

Proposed Changes

After surface 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.

Related Work

None that I know of.

PR Checklist

Put an X by all that apply. You can fill this out after submitting the PR. If you have any questions, don't hesitate to ask! We want to help. These are a guide for you to know what the reviewers will be looking for in your contribution.

I am submitting my contribution to the develop branch.
My contribution generates no new compiler warnings (try with the '-Wall -Wextra -Wno-unused-parameter -Wno-empty-body' compiler flags).
My contribution is commented and consistent with SU2 style.
I have added a test case that demonstrates my contribution, if necessary.
I have updated appropriate documentation (Tutorials, Docs Page, config_template.cpp) , if necessary.
With respect to the test case: this can be used with virtually any case that has deformation. If the scale in OPT_OBJECTIVE is set to a high enough value, the deformation will be so big that intersections will occur.

…orm 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.

pcarruscag

Thank you for opening the PR Lennaertt, some comments below regarding making use of a few things we developed recently, and regarding our modernization efforts.

Common/include/CConfig.hpp

SU2_DEF/src/SU2_DEF.cpp

Common/src/grid_movement/CVolumetricMovement.cpp

pcarruscag · 2020-09-20T16:43:32Z

Common/src/grid_movement/CVolumetricMovement.cpp

+          for (iDim = 0; iDim < nDim; iDim++) {
+            edgeVector_1[iDim] = Coords_k[iDim] - Coords_i[iDim];
+            edgeVector_2[iDim] = Coords_j[iDim] - Coords_i[iDim];
+          }
+
+          /*--- Calculate cross product of edge vectors and its length---*/
+          ev1 = edgeVector_1;
+          ev2 = edgeVector_2;
+
+          crossProduct[0] = ev1[1]*ev2[2] - ev1[2]*ev2[1];
+          crossProduct[1] = ev1[2]*ev2[0] - ev1[0]*ev2[2];
+          crossProduct[2] = ev1[0]*ev2[1] - ev1[1]*ev2[0];
+
+          for (iDim = 0; iDim < nDim; iDim++) {
+            crossProductLength += crossProduct[iDim]*crossProduct[iDim];
+          }
+
+          crossProductLength = sqrt(crossProductLength);  


There is a geometry toolbox that implements all these, please use it.

Thanks! There a few more things that the toolbox could do, like computing distance between 2 points, but we can do that latter.

pcarruscag · 2020-09-20T16:45:23Z

Common/src/grid_movement/CVolumetricMovement.cpp

+void CVolumetricMovement::ComputenNonconvexElements(CGeometry *geometry, bool Screen_Output) {
+  unsigned long iElem, PointCorners[8];


If this routine is not super expensive we should include it in the mesh quality metrics.

Common/src/grid_movement/CSurfaceMovement.cpp

Common/include/CConfig.hpp

bigfooted · 2020-09-30T17:09:09Z

Hi Lennaert, I've un-commented the items 'added a testcase' and 'added documentation'. We can use the testcase that you've worked on, and we can use one of your presentations as documentation, but it needs to be put somewhere.

bigfooted · 2020-11-25T10:09:33Z

@tollennaert, can you comment on the points raised by @pcarruscag ? I think you tried to address all points in your latest update? That makes it clear to everybody that all points have been addressed. I hope you still have time for this.

tollennaert · 2020-12-09T07:56:40Z

@tollennaert, can you comment on the points raised by @pcarruscag ? I think you tried to address all points in your latest update? That makes it clear to everybody that all points have been addressed. I hope you still have time for this.

I indeed tried to solve all issues that were mentioned earlier. Could you take another look to see whether I have done this well enough?

tollennaert · 2020-12-09T07:57:57Z

Hi Lennaert, I've un-commented the items 'added a testcase' and 'added documentation'. We can use the testcase that you've worked on, and we can use one of your presentations as documentation, but it needs to be put somewhere.

With respect to the test case and documentation: I have tried to find the best location to put this, but I could not find this. How can this be done best?

bigfooted · 2020-12-09T11:30:38Z

Hi Lennaert, I've un-commented the items 'added a testcase' and 'added documentation'. We can use the testcase that you've worked on, and we can use one of your presentations as documentation, but it needs to be put somewhere.

With respect to the test case and documentation: I have tried to find the best location to put this, but I could not find this. How can this be done best?

should we add a testcase to the SU2/TestCases/incomp_navierstokes repository?

pcarruscag · 2020-12-09T16:22:38Z

Common/include/CConfig.hpp

+  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) {


Thanks for changing the return type 👍 these 2 methods can be const

pcarruscag · 2020-12-09T16:23:09Z

Common/include/geometry/CGeometry.hpp

+   * \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;}


This one const too.

pcarruscag · 2020-12-09T16:24:47Z

Common/include/grid_movement/CSurfaceMovement.hpp

   * \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);


The compiler is probably giving a warning about this function no? (because it has the same signature but it no longer overrides, instead it "hides" the original virtual SetSurface_Deformation)

pcarruscag · 2020-12-09T16:28:34Z

SU2_DEF/src/SU2_DEF.cpp

+      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;


You can use the following syntax to initialize to zero:
TotalDeformation[iDV] = new su2double[config_container[iZone]->GetnDV_Value(iDV)] ();
(note the () at the end).

pcarruscag · 2020-12-09T16:33:00Z

SU2_DEF/src/SU2_DEF.cpp

+      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);
+        }
+      }      


There is an easier way to do this.
auto OriginalCoordinates = geometry_container[iZone]->nodes->GetCoord();
or su2activematrix instead of auto.
This: su2double OriginalCoordinates[nPoint][nDim]; is actually not standard c++

pcarruscag · 2020-12-09T16:47:26Z

Thank you @tollennaert, I added a few notes but if you don't have time we can take care of it after this is merged (I understand that your project is finished).
Regarding the testcase, ideally add it to TestCases/deformation (making it a regression, preferably in parallel_regression.py) and if possible use one of the existing grids (this should be easier than having to add a mesh file to the TestCases repo).
Finally it would be great if you could add a little more info on the options and how to use them, either in config_template or in the new testcase you add.

TobiKattmann

Hi @tollennaert, the code looks generally really tidy and thanks a lot for the contribution 💐
I have just a very few minor things, some of which are also up to debate.

With respect to the test case and documentation: I have tried to find the best location to put this, but I could not find this. How can this be done best?

For a start if you have a pdf, pptx you can share you can just dump them in here for now, so we can look at it together. And in terms of regression test I think just like pedro mentioned to take an existing testcase and modify it so it showcases your contributions. Or of course if you have a nice case at hand you can just add it as a new one.

As you have already finished your project it is ok if you dont want to spend any more time on it. But if you want to finish it: perfect, and disregard the following paragraph, otherwise:
Please allow us to push commits by checking one box on the side if it isn't enabled yet (not sure though who qualifies as 'maintainer' here). If you don't want us to push on your forks branch we could merge this and open a follow-up PR directly if changes are necessary.

Cheers, Tobi

Edit: Oh right and I didn't really understand what the depth-variables are for exactly tbh. I guess you half-en the DV_VALUE and check whether that solved the problem and try that a max of iter times. So what is the depth for 🤔

TobiKattmann · 2020-12-16T17:35:13Z

Common/src/CConfig.cpp

+  /* DESCRIPTION: Procedure to prevent self-intersections within the FFD box based on Jacobian determinant */
+  addBoolOption("FFD_INTPREV", FFD_IntPrev, NO);


I am in favor of FFD_INTERSECTION_PREVENTION just for the cfg-name not the c++ variable name (same for the other two) to be a bit more intuitive. No need to save a few chars here imo

TobiKattmann · 2020-12-16T17:36:13Z

Common/src/CConfig.cpp

+  /* DESCRIPTION: Convexity check on all mesh elements */
+  addBoolOption("CONVEXITY_CHECK", ConvexityCheck, NO);


And here I would go maybe for MESH_CONVEXITY_CHECK just to add a little more intuition

TobiKattmann · 2020-12-16T17:48:43Z

Common/src/grid_movement/CSurfaceMovement.cpp

+
+            }
+
+            /*--- Set total deformation values in config ---*/


well, not anymore in config right?

TobiKattmann · 2020-12-16T17:51:13Z

Common/src/grid_movement/CSurfaceMovement.cpp

+              /*--- 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 ---*/


👍 for including the resource

TobiKattmann · 2020-12-16T18:18:21Z

Common/src/grid_movement/CSurfaceMovement.cpp

+  unsigned long iSurfacePoints;
+  unsigned short iMarker;
+  unsigned long negative_determinants = 0;
+
+  /*--- Loop over the surface points ---*/
+  for (iSurfacePoints = 0; iSurfacePoints < FFDBox->GetnSurfacePoint(); iSurfacePoints++) {


Here one could use:

for(auto iSurfacePoints = 0ul; ...)

for unsigned long also for the i,j,kdegree loops below i guess the with 0us... but that's nothing important

TobiKattmann · 2020-12-16T18:25:57Z

Common/src/grid_movement/CSurfaceMovement.cpp

+    }
+  }
+
+  unsigned long negative_determinants_local = negative_determinants; negative_determinants = 0;


I would spend 2 lines here.

TobiKattmann · 2020-12-17T00:13:00Z

SU2_DEF/src/SU2_DEF.cpp

  for (iZone = 0; iZone < nZone; iZone++){

+    /*--- Initialize total deformation of design variables to zero ---*/
+    su2double **TotalDeformation;


Use su2activematrix instead of su2double**. from @pcarruscag prob still applies here in the new home of TotalDeformation?
You then would most likely

su2activematrix TotalDeformation; TotalDeformation.resize(nDV,nDV_Value) = su2double(0.0);

... that resizing in the loop is not that nice, so I would like to hear pedros opinion before anything is changed here

TobiKattmann · 2020-12-17T00:18:24Z

config_template.cfg

+% Parameters for prevention of self-intersections within FFD box
+FFD_INTPREV = YES
+FFD_INTPREV_ITER = 10
+FFD_INTPREV_DEPTH = 3
+


Here in the config template I would prefer to have an explanation for each option. Imo the best way is to have a suitable explanation in the CConfig.cpp and just pasting it here.

tollennaert · 2020-12-22T08:02:32Z

Hello all,

My apologies for taking so much time to respond. As much as I would have liked to finish this part of the project nicely, I am afraid I must come to the conclusion that I have taken on too many responsibilities to finish the code andfinalize the pull request.

I have attached a pptx to this comment that explains a little more about how the procedures work, the "allow edits by maintainers" box is checked. If you have any questions for me I will still be available to answer them, but I am not able to actively work on this anymore. The test cases should be available for @bigfooted, is it possible for you to place them where they are required?

As for the depth that was mentioned in an earlier comment of @TobiKattmann : the depth is an option that can be used as a stopping criterion. For every iteration without problems the recursion depth is increased. This way, it can be used as a measure for how close to the limit you are going to stop, as opposed to stopping after the maximum amount of iterations.

invalidGridPrevention.pptx

pcarruscag · 2020-12-23T22:54:13Z

Merged because it is a pain to push to branches in forks.

pcarruscag · 2020-12-23T22:54:55Z

@bigfooted please open a PR to add the testcases @tollennaert mentioned you have.

Address comments from #1076...

pr-triage bot added the PR: draft label Sep 20, 2020

pcarruscag reviewed Sep 20, 2020

View reviewed changes

lennaerttol added 2 commits November 5, 2020 11:19

Fixed comments after initial pull request

d576687

Removed tabs for spaces

743d0d9

Merge branch 'develop' into feature_invalidGridPrevention

21a7224

pcarruscag added the changelog:feature label Dec 1, 2020

pcarruscag marked this pull request as ready for review December 1, 2020 17:26

pr-triage bot added PR: unreviewed and removed PR: draft labels Dec 1, 2020

Merge branch 'develop' into feature_invalidGridPrevention

6bbd4a6

pcarruscag reviewed Dec 9, 2020

View reviewed changes

Merge branch 'develop' into feature_invalidGridPrevention

770dc94

Merge branch 'develop' into feature_invalidGridPrevention

8e0085b

TobiKattmann reviewed Dec 17, 2020

View reviewed changes

Merge branch 'develop' into feature_invalidGridPrevention

bced0cd

pcarruscag merged commit 90d156a into su2code:develop Dec 23, 2020

pr-triage bot added PR: merged and removed PR: unreviewed labels Dec 23, 2020

TobiKattmann mentioned this pull request Dec 24, 2020

Address comments from #1076... #1142

Merged

pcarruscag added a commit that referenced this pull request Dec 28, 2020

Merge pull request #1142 from su2code/address_review_1076

dd4acca

Address comments from #1076...

This was referenced Jan 5, 2021

missing testcase and doc for implementation of 'prevention of mesh intersection' #1154

Closed

add testcase for intersection prevention #1158

Merged

pcarruscag mentioned this pull request Jan 20, 2021

Bugs in CVolumetricMovement::ComputenNonconvexElements #1171

Closed

		void CVolumetricMovement::ComputenNonconvexElements(CGeometry *geometry, bool Screen_Output) {
		unsigned long iElem, PointCorners[8];

		/* DESCRIPTION: Procedure to prevent self-intersections within the FFD box based on Jacobian determinant */
		addBoolOption("FFD_INTPREV", FFD_IntPrev, NO);

		/* DESCRIPTION: Convexity check on all mesh elements */
		addBoolOption("CONVEXITY_CHECK", ConvexityCheck, NO);

Conversation

tollennaert commented Sep 20, 2020 • edited by bigfooted Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Proposed Changes

Related Work

PR Checklist

Uh oh!

pcarruscag left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

bigfooted commented Sep 30, 2020

Uh oh!

bigfooted commented Nov 25, 2020

Uh oh!

tollennaert commented Dec 9, 2020

Uh oh!

tollennaert commented Dec 9, 2020

Uh oh!

bigfooted commented Dec 9, 2020

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

pcarruscag commented Dec 9, 2020

Uh oh!

TobiKattmann left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

tollennaert commented Dec 22, 2020

Uh oh!

pcarruscag commented Dec 23, 2020

Uh oh!

pcarruscag commented Dec 23, 2020

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

tollennaert commented Sep 20, 2020 •

edited by bigfooted

Loading

TobiKattmann left a comment •

edited

Loading