fprotais / hexsmoothing

Hexahedral smoothing codes done during my Ph.D.
GNU General Public License v3.0
5 stars 2 forks source link

issue of internal smoothing with concave angle in geometry #2

Closed otaolafranc closed 6 months ago

otaolafranc commented 10 months ago

Hello @fprotais , sorry to bother again, I succesfully smooth the surface of my mesh which increased a lot the quality of the surface, nevertheless, I am having issues with the innerSmoother, the smoother is generating a bad cell in a concave angle which blocks the smoother (and well the quality it is quite bad) I have left the hole night and the smoother is 'blocked' it did already 11.6 k iterations but there are this bad cells.

iteration #11655 E: 17.5806 eps: 4.55821 detmin: -22.7911 ninv: 328421 ||g||/max(1,||x||) <= 1e-14 E : 17.5806 ---> 17.5806

here is an image of the surface mesh (white), mesh before smoothing (red) and after the smoothing (blue) for only 20 interations, and zoom in of what I imagine is the problematic cell.....

image1-26

furthemore, here is the original mesh if you want to play with it: testMesh.vtk.tar.gz

furthemore, I just noticed that it is waaaay worst than I thought here is a video of passing throught the mesh in the two cases, the red original one and the green/blueish the one outputed after several interations (after 5 if I remember corrrectly the output log does not evolve any more.

https://github.com/fprotais/hexsmoothing/assets/44319721/c4b11482-ce3c-4454-82be-b0048ac82952

fprotais commented 10 months ago

Hey there, Apologies for the slight delay, but happy new year!

E: 17.5806 eps: 4.55821 detmin: -22.7911 ninv: 328421

Indeed something went very wrong, and the culprit to look at is this very large epsilon (it is supposed to stay low!). Could you share some of the starting iteration? They could help me see what's going on. I cannot run the code right now, but could you try the following: Edit the following line: https://github.com/fprotais/hexsmoothing/blob/e31c4d862f31e1099238472504777bf50836a18b/lib/standard_smoother.cpp#L96 By replacing false by true. Then recompile using make, and try running the code again. The code might be slower in general, but significantly more robust, well at least I hope.

otaolafranc commented 10 months ago

Hello françois, no apologies needed, you are taking the time to maintain an old code for me I am more than grateful and happy new year for you too. well.... changing this flag impacted directly to the quality of the mesh actually, and it was not too much slower I would say even faster (but I think that this is due to the smoother completly breaking the mesh in the firsts iterations and then trying to recover from it takes more time...) One last thing, I saw that in standard_smoother.cpp it has also tetrahedral smoothing? Is it possible to use this code for hybrid meshes? (hexa+tetra+pyramids in the boundaries between tetra and hexa!?)

fprotais commented 10 months ago

Actually, the way I smooth hexes, is by considering them to be a set of tetrahedra! Meaning that it is almost free to do, I would just have to merge some code. Pyramids are a bit of a problem, but nothing out of reach it seems. It looks like I have code to read them, and splitting in 4 (overlapping) tets should do the tricks for smoothing.

One issue is the way I read the model, I may lose the inside/outside which is a bit annoying. Would you have a file that I could try on? I'll try to find a spot this weekend to look at it.
I won't have a way to check the quality though, I don't have a tool to display pyramids, so I'll leave that to you!

otaolafranc commented 10 months ago

thanks a lot man! here you have the most fully featured cell type that I can actually right now output. this means also prisms. FYI: my code what actually does is to: take as input a geometry block divided in a CAD format (B-rep or STEP) and accepts 3 different types of 'blocks', bricks, prisms and polyhedron blocks, where bricks are the 'classic' blocks of 6 faces, the prisms n polygon prisms, and polyhedron whatever else. the idea is that it meshes each type of block using different algorithms and different cell types, this mean structured hexas for bricks, triangular prisms for prism blocks, and tetrahedrons (and sometimes pyramids) for polyhedrons (pyramids in case that the polyhedron shares a face with a brick of one of the extruded face of a prism. here you have the output in vtk of a mesh with every type of element inside of it, the quality is already quite good so smoothing would not modify a lot but for reading file purpose it might be everything you need: mixedElementMesh.vtk.tar.gz the mesh looks like this: image with the blue elements being prisms, the green hexas, and the red section, has tetras and pyramids (as it shares one face with the green section) PS. in case that you are intereseted, what I use to visualize the meshes is paraview :) which support any kind of cell type

fprotais commented 10 months ago

Hey there, I have not forgotten you. Code is almost there, I am having a small issue in element definition.

And Paraview allows to read the file, but it is not great at displaying volumetric meshes... Or at least I don't know where to click. But it should do for my purpose.

otaolafranc commented 9 months ago

Hey there, I have not forgotten you. Code is almost there, I am having a small issue in element definition.

And Paraview allows to read the file, but it is not great at displaying volumetric meshes... Or at least I don't know where to click. But it should do for my purpose.

Hi francois, no worries from my side I finished of integrating it to my workflow and when the update to the smoother will come I will fit right into it :). here a small photo of what I was looking to achieve, smoothing the blocks of a CAD blocked geometry. still strugling with some small details but at 90% there. this allows me to easily generate meshes that are 'pre' smoothed and therefore economize the smoother step for other meshes. not sure but might help you, in paraview for better visualization, you can select the mesh object, and use the clip filter (also ctrl+space bar and write clip) and in the clip filter options click in the small gear for advanced options and turn on the crinkle cells) image

otaolafranc commented 7 months ago

Hey there, I have not forgotten you. Code is almost there, I am having a small issue in element definition.

And Paraview allows to read the file, but it is not great at displaying volumetric meshes... Or at least I don't know where to click. But it should do for my purpose.

hello @fprotais by any chance have you managed to make it work? thanks in advance!

fprotais commented 7 months ago

Hey @otaolafranc! I was pretty busy in the past months, so I was not able to find time. This should be better now, but not immediately. I should have time in the next weeks. If you have no news from me in 2 weeks, you can ping me again, I'll get you that code.

fprotais commented 7 months ago

Hi @otaolafranc, taking a look at the mixedElementsMesh.vtk you sent me, all elements seems to be not numbered correctly: loading the mesh in hexalab.net shows all hexes with a -1 scaled jacobian. Sadly, I am outside of my usual dev environment, and it is difficult for me to check why. Could you take a look at it? My parser, that follows the vtk standard, finds all the elements inverted (tets, hexes, pyramids and wedges). Otherwise my code seems ready to do mixed-smoothing!

otaolafranc commented 7 months ago

Hello @fprotais sorry for the delay, i have been having a lot of work.... well.... looks like my vtk writer still have some issues... here is a mesh that was writed in gmsh directly so it should do the work normally. sorry for the delay and the wrong mesh. example.tar.gz

fprotais commented 6 months ago

Hi @otaolafranc, A new commit is there! I only had your test case, so there might still be issues, but it smoothed. After updating the repo, you can run the following line from the build directory: cmake .. && make -j && ./mixedSmoothing ../mixedElementsMesh.vtk result.vtk And it should smooth you mixed mesh! I struggled with orientations, but I've set-up the code to easily swap between configuration. if you have any issue, it should be easy to troubleshoot.

otaolafranc commented 6 months ago

hello @fprotais suuuuper nice :)! thanks a lot for this addition! this is awsome, from my knowledge is the only smoother with mesquite that is open source and smoothes in 3D, would be nice to have a way to cite this github (maybe something to add to the github readme?) for giving the corresponding credits when I am going to publish with. also, small question, what do you think in smoothing polyhedral meshes? up to what you have created up to now it completely fully feature for my meshes, maybe it would miss the possibility of fixing nodes (and obviously surface smooth but as you mention this is super difficult). also small question, in comparison to hexSmoothing, before the command to its utilisation was: OMP_NUM_THREADS=1 ./innerSmoother mesh.vtk output.vtk 50 so one could set the number of threads and the number of iterations, is this possible now with mixedSmoothing?

but another utilisation I can see of your smoother, would be for smoothing meshes generated by a quite common automatic mesher, snappyHexMesh (https://www.openfoam.com/documentation/user-guide/4-mesh-generation-and-conversion/4.4-mesh-generation-with-the-snappyhexmesh-utility) an octree mesher, but this I imagine that would be radically difficult to do as the decomposition to tets would be super hard (as you can get meshes with awful polyhedrons (ie., cells of 21 faces! this happens during the snapping to a surface stage) it is just a question.....

thanks a lot and an awsome work :)

fprotais commented 6 months ago

I've added bibtexs to the ReadMe :) Not in order, but hopefully responding to everything:

otaolafranc commented 6 months ago

I actually forgot about OMP_NUM_THREADS, technically the optimization code I am using should be compatible with it, but every time I tried it, it ended being slower that using only 1 thread. So I would recommend actually adding OMP_NUM_THREADS=1 before the exe.

nice to know, thanks for the clarification!

Smoothing codes are indeed difficult to find, or at least the smoothing that we actually want. That is one the reason why this code is there, I had one on my side, and people asked me about it. Another "major" player is MFEM, that could be seen as the successor of Mesquite with their TMOP method. I never tried running it though.

I will give a look to TMOP, it looks awfully complex in comparation ^^

Regarding Polyhedra, my code relies on tetrahedra. My first guess would be to "subdivide" each Polyhedron in a set of tetrahedra, so that I can smooth them (what I did with hexes, pyramids and wedges). If the Polyhedra are very complex, it could be a very poor way of doing it. I would say that the best is to introduce the smoothing energy directly in the metric used to produce de Polyhedra, as it would directly solve the right problem. The framework proposed by TMOP is more general than mine, it is possible that they have tools to smooth polyhedra (I know that they have for high order meshes).

yeah, that is the thing, transforming to tetra with such complex cells I would imagine that would be quite complex to do, at least automatically....

Regarding fixing point, it rather trivial to do so using the cpp interface, but would take a bit of work to set-up using input files. If you have a precise use case we can discuss that

maybe a flag field inside the vtk file as in the case of mesquite?

And lastly: I tried playing a bit with smoothing meshes from other people. The main issue I ran into is that problems usually are on the boundary, where my code doesn't work the best... Hopefully I will have time in the future to iron that thing out. And a word of caution about Octree meshes: the large difference in scale of elements introduces bad things in my untangling approach. If you start with an invalid mesh, it may struggle a lot to find a good solution. A bit more research should be done there.

well, if you are interested in developing this I think I have an idea on how to do it without a lot or required files but it would be a lot of developpment, but in my head, having a scalar field inside the vtk file which would be equal to: 0 if it is frelly to move, -1 if it is fixed, and then a scalar [1,inf) where the node would be associated with a bunch of step files which would be called 1.step, 2.step etc, where for example 1, is an edge/polyline (saved in .step/.iges or a highly refined surface mesh) and therefore it could move over this edge/polyline, then another vertex has 2, and it could move over a surface... even blocking nodes inside a specific volume. Then create a folder associated with the vtk file or something in that regard.... and if there isnt a folder associated, all nodes with a index different to 0 would be interpreted as -1 (so this would 'fix' the nodes when the folder is not detected,and therefore would work as current smoother ie., internal smoothing without touching the external nodes). but yeah, this is my current issue, I have succesfully done some surface smoothing pre smoothing the internal cells (in the input vtk of hexsmoothing) but it would be better to have at the same time as the smoother I am using gives some issues, and also the results of the internal and external mesh would be different if they are coupled

fprotais commented 6 months ago

maybe a flag field inside the vtk file as in the case of mesquite?

My parser can't parse flags currently, but it could be done. You can raise a new issue, and when I have a strike of motivation, I can do that. It shouldn't be too problematic if I focus on only 1 flag. (And if you could provide me with an input file for Mesquite it could be really useful, I don't even know how to generate those)

Surface stuff

There is two major issues in surface meshing:

  1. matching of features
  2. Smoothing on the surface and in the volume at the same time

For issue 1, when you have already nicely labelled data, you indeed provide a good solution. But without it, it is a very complex poorly defined problem. SnappyHexMesh probably has heuristic for that, but none are perfect, as there is no reason for a given mesh to match with the complexity of the initial features. Adding to that, some matching may actually be impossible to smooth because there exists no valid geometrical solution. That is the major issue when trying to smooth meshes from other methods: they may be loosely related to original feature set.

For issue 2, there is two sides of it. I could provide an iterative process (boundary then volume then boundary), but this requires a lot of fine tuning and study of different cases. Those methods have proven to be decent, but it is a lot of work to get right. The second possibility is too actually smooth both at the same time. Now this gets tricky: You have to parametrize points on surface and somehow smoothly optimize your energy on those optimized point. I had people tell me it shouldn't be too hard to do, but I never saw a working code :)

All in all, as this was research code, the best would have be to try to go with the second solution, but it is a. mathematically a complex thing, and b. would time consuming to achieve (a lot of testing involved to get the thing right). And most importantly, while all have done here is just re-using my existing available code, this would involve a new important coding project, which I am currently not in a position to do! (or at least open-source) If ever you want to get your hand deep in smoothing codes, I can give some pointer on what I think can be good directions.

otaolafranc commented 6 months ago

Hello @fprotais I might have missexplained, when I was talking about the surface smoothing, I was imagining not meshes generated with something like SHM which is a an automatic mesher that begins from a background mesh, but the meshes I compose myself in salome (or in gmsh), where you begin from a CAD geometry, so what I was saying is reading a vtk file (as what I gave you) with a dataset attribute (https://docs.vtk.org/en/latest/design_documents/VTKFileFormats.html) so in a vtk file you define the POINTS (nodes of mesh) then CELLS (the index of the nodes used for each element of the mesh) then CELL TYPES (where you indicate the type of element), but one can have a also a scalar which would be a table of scalar that gives a value to each element of POINTS (to 'attach' the node to a geometry element (I mean a STEP file of a node/edge/surface/volume) ) and the smoother would check if each element of the POINTS is on (or really near ~error distance) of the respective STEP file. (this could be done automatically in salome or gmsh) this would be a lot of development.... but if you are interested I could prepare an example of a file type. for something like the geometry I uploaded. I wait for your answer if you are interested in something like this and if yes I will create a new issue, but understand if you don't want to develop further as it would be a lot of work, in any case I really appreciate! I will not doubt in citing the corresponding credits :) regards!

fprotais commented 6 months ago

What you propose sadly only cover my point 1, not point 2. Currently, I am not in a position to do any open-source development for this second one (in a way that would compete with my current position in Siemens). If that code is of real interest to you, I would recommend reaching out to the TMOP people and see if they have anything matching your need! Unlike my small code that I made alone, they have actual aim to propose functional codes to users in the MFEM project, and have resources allocated for that. Cheers :)