Difficulty setting manifold with Opencascade using gmsh mesh + salome STEP (or IGES) - Step-54
202 views
Skip to first unread message
Bruno Blais
Mar 31, 2020, 9:19:23 AM3/31/20
to deal.II User Group
Dear all,
I hope you are well in these uncertain times.
I have been trying to use the OpenCASCADE library to set-up my manifolds. The goal of this is to be able to do CFD simulations in complex geometry, but use an initial very coarse mesh made with GMSH and use the dynamics mesh adaptation to adapt the mesh while having the adaptation be CAD-aware. The geometry is complex, so fixing the manifold analytically appears to be a bad idea.
I started from step-54, which compiles and works well on my machines.
However, I am unable to "make it work" with my own generated Salome STEP and gmsh MSH files.
My CAD is in metric, so I adjusted the scaling to use a scaling of 1. My mesh also contains numerous wires, so I disabled that aspect and only focused on refining with the faces which I label by looping over them.
However, no matter what occurs I get an error thrown:
--------------------------------------------------------
An error occurred in line <114> of file </home/blab/soft/dealii/dealii/dealii/source/opencascade/manifold_lib.cc> in function
dealii::Point<spacedim> dealii::OpenCASCADE::NormalProjectionManifold<dim, spacedim>::project_to_manifold(const dealii::ArrayView<const dealii::Point<spacedim> >&, const dealii::Point<spacedim>&) const [with int dim = 2; int spacedim = 3]
The violated condition was:
closest_point(sh, surrounding_points[i], tolerance) .distance(surrounding_points[i]) < std::max(tolerance * surrounding_points[i].norm(), tolerance)
Additional information:
The point [ 0.1 0 0 ] is not on the manifold.
Stacktrace:
-----------
#0 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::OpenCASCADE::NormalProjectionManifold<2, 3>::project_to_manifold(dealii::ArrayView<dealii::Point<3, double> const, dealii::MemorySpace::Host> const&, dealii::Point<3, double> const&) const
#1 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::FlatManifold<2, 3>::get_new_point(dealii::ArrayView<dealii::Point<3, double> const, dealii::MemorySpace::Host> const&, dealii::ArrayView<double const, dealii::MemorySpace::Host> const&) const
#2 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::Manifold<2, 3>::get_new_point_on_line(dealii::TriaIterator<dealii::TriaAccessor<1, 2, 3> > const&) const
#3 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre:
#4 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre:
#5 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::TriaAccessor<1, 2, 3>::center(bool, bool) const
#6 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::Triangulation<2, 3>::DistortedCellList dealii::internal::TriangulationImplementation::Implementation::execute_refinement<3>(dealii::Triangulation<2, 3>&, bool)
#7 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::Triangulation<2, 3>::execute_refinement()
#8 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::Triangulation<2, 3>::execute_coarsening_and_refinement()
#9 /home/blab/soft/dealii/dealii/build/lib/libdeal_II.g.so.9.2.0-pre: dealii::Triangulation<2, 3>::refine_global(unsigned int)
#10 ./step-54: Step54::TriangulationOnCAD::refine_mesh()
#11 ./step-54: Step54::TriangulationOnCAD::run()
#12 ./step-54: main
--------------------------------------------------------
I have tried some elements :
- Playing with the tolerance of the CAD
- Writing my CAD from the OpenCASCADE STEP object and then opening it again in Salome to see that it is interpreted correctly +
However, I always seem to face an issue where my CAD does not appear to generate a valid manifold. I have joined my code as well as the salome HDF file and gmsh GEO file used to generate the CAD and meshes respectively.
I must admit I am quite at lost here as to why both don't coincide. Any advices on where to proceed from there?
Best :)!
Bruno
Message has been deleted
Message has been deleted
Bruno Blais
Apr 1, 2020, 9:38:27 AM4/1/20
to deal.II User Group
So it seems part of my problem relate to the topology I was trying to adapt OR to the staring complex mesh I was using.
I have managed to make it work starting with a trivial mesh and using a more simple IGES CAD. I will work on complexifying it from there. I will try to keep the information here as it seems that feature has not been used extensively :)!.
What I have working now is this animation.
Bruno Blais
Apr 1, 2020, 11:32:06 AM4/1/20
to deal.II User Group
So my investigation on this issue continues...
It seems the core of my issues stem from using a non-trivial mesh with more than one cell as a starting point.
if I start with a 2 cell spacedim=2 dim=3 mesh such as:
Which is identical to the second step of the adaptation. I get an aberrant result (see attached out_wrong.mp4).
I also get issues when I try to apply the manifold on a starting msh made with GMSH that has a similar topology such as in the following image:
In this case, I get an error about a point not being on the manifold (when it clearly is)...
So I must admit I am at loss here. Right now the code I use to set the manifold is :
for (const auto &cell : tria.active_cell_iterators())
{
cell->set_manifold_id(1);
for (const auto &face : cell->face_iterators())
{
face->set_manifold_id(2);
}
}
If anybody has any experience with using the OpenCascade manifold using initial meshes that have more than 1 cell, I would greatly appreciate advices :)!.
Best
Bruno
luca.heltai
Apr 2, 2020, 11:31:42 AM4/2/20
to Deal.II Users
Bruno, it seems like you are attaching your faces to the *boundary* manifold. Let me try to explain what is happening.
In the code for step-54, there is a wire that is used to describe the manifold of the *boundary* of the surface (a curve of dimension one embedded in dimension three). This is generated by the wire that identifies the boundary of the TopoDS shape of the face. The manifold attached to it is an ArcLengthProjectionManifold, where mid points on the curve are added by computing the distance in arc length, and taking the point in the middle w.r.t. this length.
In your code with more than one cell, you are not assigning *just the boundary faces* to the wire, but *all* faces (including the one in the middle of your surface, delimitating the two cells). Intermediate points on those lines will be projected (correctly) on the boundary wire (that’s what your movie shows).
Of course any point that is in the middle of the surface does *not* belong to the *boundary wire*, and you’ll get an exception when trying to construct the midpoint of an *edge* that belongs to the interior of the domain, with manifold id 2 (which corresponds to a manifold that describes *the boundary* of your topology, not its surface).
With two cells, the internal edge (the only edge that does not belong to the boundary) should *not* have id 2. If you set its id to 2, the result is that a point which is in the middle of the edge, is actually in the midlle of the two points *when running along the curve that describes the boundary*, hence you get the distortion you show in the movie.
Try adding a check to the faces in which you set the manifold id to 2 only if the face is at the boundary.
The principle is:
1. start from the lowest codimension objects, identify how to deform them. In your case, cells of a Tria<2,3> are quads, that should deform as a TopoDS_FACE.
2. Attach your favorite manifold to the TopoDS_Face (I’d personally only use NormalToMeshProjectionManifold) using cell->set_all_manifold_ids(manifold_id) (Notice the use of set_all_manifold_ids, and **not** set_manifold_id: you want all children of these objects to belong to the same manifold, in particular you want all faces that are internal, i.e., that are shared between two obects with the same manifold id, to inherit the same manifold id).
3. Go one codimension lower: in your case, curves (for 3d meshes, surfaces). Follow the same rule as above, setting all_manifold_ids on faces that you know should follow a known codimension one geometry (a known TopoDS_EDGE, or TopoDS_WIRE). In this case I’d only use ArcLengthProjectionManifold objects, attached to the wires that identify your geometry.
Doing things in this order guarantees that internal edges get the correct manifold id, which, in your case, is not happening.
Best,
Luca.
> On 1 Apr 2020, at 17:32, Bruno Blais <blais...@gmail.com> wrote:
>
> So my investigation on this issue continues...
> It seems the core of my issues stem from using a non-trivial mesh with more than one cell as a starting point.
> if I start with a 2 cell spacedim=2 dim=3 mesh such as:
>
>
> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/9decc8c5-0a99-44f1-aaa9-24deb26edc06%40googlegroups.com.
> <out_wrong.mp4><start_point.png><starting_mesh.png>
In the code for step-54, there is a wire that is used to describe the manifold of the *boundary* of the surface (a curve of dimension one embedded in dimension three). This is generated by the wire that identifies the boundary of the TopoDS shape of the face. The manifold attached to it is an ArcLengthProjectionManifold, where mid points on the curve are added by computing the distance in arc length, and taking the point in the middle w.r.t. this length.
In your code with more than one cell, you are not assigning *just the boundary faces* to the wire, but *all* faces (including the one in the middle of your surface, delimitating the two cells). Intermediate points on those lines will be projected (correctly) on the boundary wire (that’s what your movie shows).
Of course any point that is in the middle of the surface does *not* belong to the *boundary wire*, and you’ll get an exception when trying to construct the midpoint of an *edge* that belongs to the interior of the domain, with manifold id 2 (which corresponds to a manifold that describes *the boundary* of your topology, not its surface).
With two cells, the internal edge (the only edge that does not belong to the boundary) should *not* have id 2. If you set its id to 2, the result is that a point which is in the middle of the edge, is actually in the midlle of the two points *when running along the curve that describes the boundary*, hence you get the distortion you show in the movie.
Try adding a check to the faces in which you set the manifold id to 2 only if the face is at the boundary.
The principle is:
1. start from the lowest codimension objects, identify how to deform them. In your case, cells of a Tria<2,3> are quads, that should deform as a TopoDS_FACE.
2. Attach your favorite manifold to the TopoDS_Face (I’d personally only use NormalToMeshProjectionManifold) using cell->set_all_manifold_ids(manifold_id) (Notice the use of set_all_manifold_ids, and **not** set_manifold_id: you want all children of these objects to belong to the same manifold, in particular you want all faces that are internal, i.e., that are shared between two obects with the same manifold id, to inherit the same manifold id).
3. Go one codimension lower: in your case, curves (for 3d meshes, surfaces). Follow the same rule as above, setting all_manifold_ids on faces that you know should follow a known codimension one geometry (a known TopoDS_EDGE, or TopoDS_WIRE). In this case I’d only use ArcLengthProjectionManifold objects, attached to the wires that identify your geometry.
Doing things in this order guarantees that internal edges get the correct manifold id, which, in your case, is not happening.
Best,
Luca.
> On 1 Apr 2020, at 17:32, Bruno Blais <blais...@gmail.com> wrote:
>
> So my investigation on this issue continues...
> It seems the core of my issues stem from using a non-trivial mesh with more than one cell as a starting point.
> if I start with a 2 cell spacedim=2 dim=3 mesh such as:
>
>
> <start_point.png>
>
>
>
> Which is identical to the second step of the adaptation. I get an aberrant result (see attached out_wrong.mp4).
>
>
>
> I also get issues when I try to apply the manifold on a starting msh made with GMSH that has a similar topology such as in the following image:
>
>
>
>
> Which is identical to the second step of the adaptation. I get an aberrant result (see attached out_wrong.mp4).
>
>
>
> I also get issues when I try to apply the manifold on a starting msh made with GMSH that has a similar topology such as in the following image:
>
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/9decc8c5-0a99-44f1-aaa9-24deb26edc06%40googlegroups.com.
> <out_wrong.mp4><start_point.png><starting_mesh.png>
Bruno Blais
Apr 3, 2020, 12:04:00 AM4/3/20
to deal.II User Group
Ciao Luca,
It works perfectly now even with regular GMSH triangulation (see movie). Now I understand the difference and how to set it correctly.
Thank you very much for the help :)!
> To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.
Jean-Paul Pelteret
Apr 3, 2020, 9:43:33 AM4/3/20
to dea...@googlegroups.com
Dear Bruno,
I’m glad that you managed to solve this issue with Luca’s aid. I think that the information and workflow that Luca gave was really valuable. Would you care to add this to the Wiki page?
Best,
Jean-Paul
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7cead34c-5bef-4cfb-84f6-39de20a54936%40googlegroups.com.
<working.mp4>
Bruno Blais
Apr 3, 2020, 9:55:36 AM4/3/20
to deal.II User Group
Dear Jean-Paul,
Do you think it should be added to the wiki page or as additional information to step-54? I have a bit on my hands, but i'd be glad to work to make the information readily available. Clearly I would not have found that out without Luca's help.
Do you think it should be added to the wiki page or as additional information to step-54? I have a bit on my hands, but i'd be glad to work to make the information readily available. Clearly I would not have found that out without Luca's help.
Best
Bruno
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/7cead34c-5bef-4cfb-84f6-39de20a54936%40googlegroups.com.
<working.mp4>
Jean-Paul Pelteret
Apr 3, 2020, 4:11:05 PM4/3/20
to dea...@googlegroups.com
Dear Bruno,
I guess that it would be best to mention this in the tutorial itself, but in the mean time I wrote a segment about it in the FAQ:
Please feel free to edit it as you see fit. Never having used this functionality in the library before, I hope that I have correctly interpreted and massages what Luca had detailed.
Best,
Jean-Paul
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/2867bda9-9820-46e5-966f-580741c817c2%40googlegroups.com.
Bruno Blais
Apr 5, 2020, 12:05:17 PM4/5/20
to deal.II User Group
Dear Jean-Paul, what you wrote in the FAQ is very clear. If I find additional elements as I play with the functionnality I will add it to the FAQ. Thanks for adding this section.
Best
Bruno
To unsubscribe from this group and stop receiving emails from it, send an email to dea...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/2867bda9-9820-46e5-966f-580741c817c2%40googlegroups.com.
Bruno Blais
Apr 8, 2020, 10:30:00 AM4/8/20
to deal.II User Group
Dear all,
Thank you very much for the help and support. I have managed to make it work in 3D within Lethe.
A quick image of the application of CAD mesh adaptation using multiple CADs on a complex mixer :)!
Thank you for your help (especially you Luca :)!)
luca.heltai
Apr 9, 2020, 5:07:14 AM4/9/20
to Deal.II Users
That’s fanstastic!
:)
You sould put this image in the image gallery!
Luca.
> On 8 Apr 2020, at 16:29, Bruno Blais <blais...@gmail.com> wrote:
>
> Dear all,
> Thank you very much for the help and support. I have managed to make it work in 3D within Lethe.
> A quick image of the application of CAD mesh adaptation using multiple CADs on a complex mixer :)!
> Thank you for your help (especially you Luca :)!)
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6548f917-b8b3-407c-9b9a-66f0205d578d%40googlegroups.com.
> <quarter_ribbon.0004.png>
:)
You sould put this image in the image gallery!
Luca.
> On 8 Apr 2020, at 16:29, Bruno Blais <blais...@gmail.com> wrote:
>
> Dear all,
> Thank you very much for the help and support. I have managed to make it work in 3D within Lethe.
> A quick image of the application of CAD mesh adaptation using multiple CADs on a complex mixer :)!
> Thank you for your help (especially you Luca :)!)
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/6548f917-b8b3-407c-9b9a-66f0205d578d%40googlegroups.com.
> <quarter_ribbon.0004.png>
Jean-Paul Pelteret
Apr 9, 2020, 6:37:54 PM4/9/20
to dea...@googlegroups.com
Seconded! It’s a really nice result. I’m glad that you managed to solve your problems, Bruno.
Best,
Jean-Paul
> To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/96520D42-D321-4D23-AE8B-8ACB9612616B%40gmail.com.
Best,
Jean-Paul
Reply all
Reply to author
Forward
0 new messages