Add snap_to_tessellated_mesh #1

masterleinad · 2024-04-19T21:10:16Z

This pull request adds a tool that snaps a tessellated vtk mesh to the geometries surface given by an iges file.
Note that orientation and scaling between the two must match which is not true for the files initially tested on. That's why there are a bunch of comments for curved wall and hourglass.

Rombur

This could use more comments, I don't think I followed everything you were doing. Add the .clang-format from adamantine to the repo.

Rombur · 2024-04-22T14:07:56Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+
+#include <deal.II/opencascade/manifold_lib.h>
+
+dealii::Point<3> my_closest_point(const std::vector<TopoDS_Face> &faces,


Why are you using your own function and not closest_point() directly?

Because we might want to use only part of the Surface to avoid, e.,g., the top and bottom face.

Rombur · 2024-04-22T14:10:27Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+    for (unsigned int i = 0; i < normals.size(); ++i) {
+      for (unsigned int j = i + 1; j < normals.size(); ++j) {
+        auto product = normals[i] * normals[j];
+        minimum_product = std::min(product, minimum_product);


What does this do?

We are trying to figure out if all the normals point roughly in the same direction. Since the input file is tessellated, the products are either 0,1, or -1.
If the normals point in different directions, we identify the vertex as a corner and treat it differently.

Rombur · 2024-04-22T14:13:10Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+    if (minimum_product > .5) {
+      vertex(0) = proj(0);
+      vertex(1) = proj(1);
+      // vertex(2) = proj(2); // hourglass


Why are you getting the z component if you plan to ignore it anyway? Wouldn't it make more sense to search for the closest point in the xy-plane instead?

As the description in the pull request said, none of the iges files I was working with, was entirely conforming with the respective vtk file.
For the hourglass mesh, the iges file is rotated such that the printing direction coincides with the y-axis instead of the z-axis. Thus, we we need to use coordinates 0 and 2 instead of 0 and 1 for that example.

Is this just a case where the iges is aligned with a x/y/z axis or sometimes it's totally random and you need to perform a rotation to match the iges and the vtk file?

For the files provided, I needed to figure out the transformation but ultimately the user should make sure that the files match.

Why did the file not match and would one ensure that they do match?

There is not all that much we can do to detect if the files match apart from checking that the bounding boxes are roughly the same. For the curved wall example, I needed to scale the mesh by a factor of 100 and for the hourglass mesh swap y- and z-component and shift the mesh by 0.07 in x-direction.

All this shouldn't be a problem if STL file (to create the VTK file) and the IGES file are created at the same time.
As I said earlier, we could try to create the VTK file in deal.II from the STL file avoiding OpenFOAM which would also help avoiding these problems (but takes obviously more effort).

Rombur · 2024-04-22T14:14:14Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+
+  bool exclude_z_faces = false;
+
+  // We don't care about faces with z-normal so exclude them here


Is that always true or we just don't care about the bottom and the top face?

The code ignores TopoDS::Faces that are planes with z-normal. This might or might not be a good choice depending on the input files. It should not be detrimental if the face is indeed planar but might be problematic if we wrongly detect it to be planar.

In general, this should help mostly with the edges at the top and bottom.

Rombur · 2024-04-22T14:16:39Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+        auto point_0 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u1, v1));
+        auto point_1 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u1, v2));
+        auto point_2 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u2, v1));
+        auto point_3 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u2, v2));


What are these points?

The four corners of the face.

Rombur · 2024-04-22T14:17:12Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+        double deviation_from_z =
+            normal / normal.norm() * dealii::Tensor<1, 3>({0, 0, 1});
+
+        if (std::abs(deviation) > .1 || std::abs(deviation_from_z) < .9)


How did you pick these numbers?

I just picked something that worked for the meshes provided. I would need more examples to find better numbers.

masterleinad

I'll add more comments.

masterleinad · 2024-04-22T14:44:08Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+
+#include <deal.II/opencascade/manifold_lib.h>
+
+dealii::Point<3> my_closest_point(const std::vector<TopoDS_Face> &faces,


Because we might want to use only part of the Surface to avoid, e.,g., the top and bottom face.

masterleinad · 2024-04-22T14:46:26Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+    for (unsigned int i = 0; i < normals.size(); ++i) {
+      for (unsigned int j = i + 1; j < normals.size(); ++j) {
+        auto product = normals[i] * normals[j];
+        minimum_product = std::min(product, minimum_product);


We are trying to figure out if all the normals point roughly in the same direction. Since the input file is tessellated, the products are either 0,1, or -1.
If the normals point in different directions, we identify the vertex as a corner and treat it differently.

masterleinad · 2024-04-22T14:49:10Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+    if (minimum_product > .5) {
+      vertex(0) = proj(0);
+      vertex(1) = proj(1);
+      // vertex(2) = proj(2); // hourglass


As the description in the pull request said, none of the iges files I was working with, was entirely conforming with the respective vtk file.
For the hourglass mesh, the iges file is rotated such that the printing direction coincides with the y-axis instead of the z-axis. Thus, we we need to use coordinates 0 and 2 instead of 0 and 1 for that example.

masterleinad · 2024-04-22T14:51:40Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+        auto point_0 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u1, v1));
+        auto point_1 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u1, v2));
+        auto point_2 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u2, v1));
+        auto point_3 = dealii::OpenCASCADE::point<3>(SurfToProj->Value(u2, v2));


The four corners of the face.

masterleinad · 2024-04-22T14:53:25Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+
+  bool exclude_z_faces = false;
+
+  // We don't care about faces with z-normal so exclude them here


The code ignores TopoDS::Faces that are planes with z-normal. This might or might not be a good choice depending on the input files. It should not be detrimental if the face is indeed planar but might be problematic if we wrongly detect it to be planar.

In general, this should help mostly with the edges at the top and bottom.

masterleinad · 2024-04-22T14:55:56Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+        double deviation_from_z =
+            normal / normal.norm() * dealii::Tensor<1, 3>({0, 0, 1});
+
+        if (std::abs(deviation) > .1 || std::abs(deviation_from_z) < .9)


I just picked something that worked for the meshes provided. I would need more examples to find better numbers.

Rombur

The code looks ok but we still need:

a README to explain what the code does, how to install it, and how to run it
a test to make sure that we don't break the code when we modify it
a Dockerfile and a GitHub workflow to create a docker image when we push on master because most people are not going to compile deal.II

You can do this in this PR or another one.

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

Rombur · 2024-04-23T12:47:50Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+    if (minimum_product > .5) {
+      vertex(0) = proj(0);
+      vertex(1) = proj(1);
+      // vertex(2) = proj(2); // hourglass


Why did the file not match and would one ensure that they do match?

Rombur · 2024-04-23T12:48:19Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+  B.Get(Xmin, Ymin, Zmin, Xmax, Ymax, Zmax);
+
+  std::cout << "Bounding Box: (" << Xmin << ',' << Ymin << ',' << Zmin << "), ("
+            << Xmax << ',' << Ymax << ',' << Zmax << ")\n";


Any reason to write this on the screen?

It helps a lot when identifying if the IGES file and the VTK file match. Since most viewers can't display IGES files natively.

Rombur · 2024-04-23T12:49:06Z

snap_tessellated_mesh/snap_to_tessellated_mesh.cpp

+  // if the face is indeed planar but might be problematic if we wrongly detect
+  // it to be planar. Excluding those faces should mostly help with the edges at
+  // the top and bottom of the mesh (assuming thet are planar).
+  bool exclude_z_faces = false;


This should be a command line parameter.

Co-authored-by: Bruno Turcksin <[email protected]>

masterleinad · 2024-06-18T18:07:43Z

The only thing missing here is automatic testing. This pull requests includes tests and a Dockerfile that installs all dependencies.

Add snap_to_tessellated_mesh

296c8f0

Rombur reviewed Apr 22, 2024

View reviewed changes

masterleinad commented Apr 22, 2024

View reviewed changes

masterleinad added 2 commits April 22, 2024 11:34

Add more comments

ea32215

Add .clang-format file

cdf52c0

Rombur reviewed Apr 23, 2024

View reviewed changes

masterleinad and others added 7 commits April 23, 2024 09:07

Fix typos

0e43933

Co-authored-by: Bruno Turcksin <[email protected]>

Add README.md

1e90ef7

Output STL file and bounding box fir VTK file

958ee18

Remove special cases for curved wall and hourglass

8b75db6

Debug output for bouding boxes

c466976

Add Dockerfile

bbe0069

Add test files

70fe2d7

masterleinad force-pushed the snap_to_tessellated_mesh branch from 95a9a5f to 70fe2d7 Compare May 24, 2024 15:51

Make exclude_z_faces a runtime parameter

bb294ea

masterleinad changed the title ~~[WIP] Add snap_to_tessellated_mesh~~ Add snap_to_tessellated_mesh Jun 18, 2024


		#include <deal.II/opencascade/manifold_lib.h>

		dealii::Point<3> my_closest_point(const std::vector<TopoDS_Face> &faces,


		bool exclude_z_faces = false;

		// We don't care about faces with z-normal so exclude them here

Add snap_to_tessellated_mesh #1

Are you sure you want to change the base?

Add snap_to_tessellated_mesh #1

Uh oh!

Conversation

masterleinad commented Apr 19, 2024

Uh oh!

Rombur left a comment

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!

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!

masterleinad left a comment

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!

Rombur left a comment

Choose a reason for hiding this comment

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!

masterleinad commented Jun 18, 2024

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

2 participants