gmsh v4

[ryancoe]

It seems that meshmagick cannot read .msh files created with with the latest version of gmsh (4.*). I'm aware that you can use gmsh's command line option -format msh22 to output a v2.2 .msh file, but I haven't found a way to do this from the gmsh Python API. It would be great to be able to use gmsh and meshmagick python APIs together, but I'm not sure how to do this with the current version incongruence. Any thoughts?

[mancellin]

Thanks, I wasn't aware that there is a new version of gmsh files. After a quick look, it does not seem very different from the previous one, so it should be possible to update the code in meshmagick.

As for the gmsh python API, I'm not sure what you have in mind. Would you mind giving an example of how you'd like to use it?

[ryancoe]

Matthieu - Thanks for your response. I was essentially hoping to use gmsh to generate parametric geometries and meshes, and then use meshmagick to prep/check those meshes before analyzing them in Nemoh.

Ryan

--

Ryan Coe, PhD Water Power Technologies (Org. 8822) | Sandia National Laboratories Email: Ryan.Coe@sandia.govmailto:Ryan.Coe@sandia.gov | Office: (505)845-9064<tel:(505)845-9064> | Mobile: (513)910-8809<tel:(513)910-8809>

On Mar 20, 2019, at 8:44 AM, Matthieu Ancellin notifications@github.com<mailto:notifications@github.com> wrote:

Thanks, I wasn't aware that there is a new version of gmsh files. After a quick look, it does not seem very different from the previous one, so it should be possible to update the code in meshmagick.

As for the gmsh python API, I'm not sure what you have in mind. Would you mind giving an example of how you'd like to use it?

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/LHEEA/meshmagick/issues/9#issuecomment-474862725, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AFVLf7VLYrCc5m7H_YZ-IGVy-KKfbIRpks5vYklegaJpZM4bsWka.

[mancellin]

I took a look at pygmsh. It seems fairly easy to interface it with meshmagick:

# [...] define a mesh with pygmsh 

points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

def all_faces_as_tetra(cells):
    triangles_as_tetra = np.empty((cells['triangle'].shape[0], 4), dtype=np.int)
    triangles_as_tetra[:, :3] = cells['triangle'][:, :]
    triangles_as_tetra[:, 3] = cells['triangle'][:, 2]  # Repeat one node to make a tetra
    return np.concatenate([cells['tetra'], triangles_as_tetra])

from meshmagick.mesh import Mesh
mesh = Mesh(vertices=points, faces=all_faces_as_tetra(cells))
mesh.show()

pygmsh seems to be compatible with Python 2.7, so you should be able to use it with meshmagick. I have also an experimental fork of meshmagick for Python 3 available at conda install -c mancellin meshmagick. Or was it something else that you meant by "version incongruence"?

For the interface with Nemoh, you might be interested in this project for Python 3.

[ryancoe]

Matthieu - Wow, neat! Thanks very much.

I’m running OSX, so I’m having a bit of trouble getting Capytaine working… I will see if I can’t get this working.

Ryan

--

Ryan Coe, PhD Water Power Technologies (Org. 8822) | Sandia National Laboratories Email: Ryan.Coe@sandia.govmailto:Ryan.Coe@sandia.gov | Office: (505)845-9064 | Mobile: (513)910-8809

On Mar 20, 2019, at 10:37 AM, Matthieu Ancellin notifications@github.com<mailto:notifications@github.com> wrote:

I took a look at pygmshhttps://github.com/nschloe/pygmsh. It seems fairly easy to interface it with meshmagick:

[...] define a mesh with pygmsh

points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

def all_faces_as_tetra(cells): triangles_as_tetra = np.empty((cells['triangle'].shape[0], 4), dtype=np.int) triangles_as_tetra[:, :3] = cells['triangle'][:, :] triangles_as_tetra[:, 3] = cells['triangle'][:, 2] # Repeat one node to make a tetra return np.concatenate([cells['tetra'], triangles_as_tetra])

from meshmagick.mesh import Mesh mesh = Mesh(vertices=points, faces=all_faces_as_tetra(cells)) mesh.show()

pygmsh seems to be compatible with Python 2.7, so you should be able to use it with meshmagick. I have also an experimental fork of meshmagick for Python 3 available at conda install -c mancellin meshmagick.

For the interface with Nemoh, you might be interested in this projecthttps://github.com/mancellin/capytaine for Python 3.

— You are receiving this because you authored the thread. Reply to this email directly, view it on GitHubhttps://github.com/LHEEA/meshmagick/issues/9#issuecomment-474919543, or mute the threadhttps://github.com/notifications/unsubscribe-auth/AFVLfwmb9QEWacHwNDyng3ku0kjH7XTIks5vYmO8gaJpZM4bsWka.

[mancellin]

There is no Capyitaine package for OSX because I don't have a device to test it. But I don't see why it wouldn't work. Feel free to open an issue in Capytaine repository if you'd like some help. \end{off-topic}

Although it's possible to update the code of Meshmagick to read Gmsh v4 files, it might make more sense to put the effort on a binding with meshio instead.

[ryancoe]

FYI - I played a bit more with pygmsh. It is neat, but it seems that you cannot control/refine the mesh: https://github.com/nschloe/pygmsh/issues/149. This is somewhat of a showstopper for using it with BEM...

[ryancoe]

I returned this after a long hiatus and realized the pygmsh can indeed do refinement via the extra_gmsh_arguments arg for the method generate_mesh. However, I don't know (yet) how to do this without creating an STL file and reading it in (i.e., should be able to just pass the faces and vertices).

import pygmsh
import numpy as np
import capytaine as cpt
import logging
import matplotlib.pyplot as plt

logging.basicConfig(level=logging.INFO,
                    format="%(levelname)s:\t%(message)s")

geom = pygmsh.opencascade.Geometry()

fbname = 'test'
ofst = 0.1
geom.add_cylinder(x0=[0.0, 0.0, 0+ofst],
                         axis=[0.0, 0.0, -2],
                         radius=1,
                         angle=2 * np.pi,
                         char_length=1)

mshRefFactor = 0.25
mshArgs = ['-clscale', str(mshRefFactor),
           '-clcurv', str(360/50)]
mesh = pygmsh.generate_mesh(geom,
                            dim=2,
                            extra_gmsh_arguments=mshArgs,
                            remove_lower_dim_cells=True,
                            geo_filename=fbname + '.geo',
                            msh_filename=fbname + '.stl',
                            mesh_file_type='stl')

fb = cpt.FloatingBody.from_file(fbname + '.stl',
                                file_format='stl',
                                name=fbname)
fb.keep_immersed_part()
fb.show()

f_range = np.logspace(-1, 0.4, 50)
omega_range = 2*np.pi*f_range

fb.add_all_rigid_body_dofs()
problems = [
    cpt.RadiationProblem(body=fb, radiating_dof=dof, omega=omega)
    for dof in fb.dofs
    for omega in omega_range
]
solver = cpt.BEMSolver()
results = [solver.solve(pb) for pb in sorted(problems)]
data = cpt.assemble_dataset(results)

plt.figure()
for dof in fb.dofs:
    plt.plot(
        f_range,
        data['added_mass'].sel(radiating_dof=dof, influenced_dof=dof),
        label=dof,
        marker='o',
    )
plt.xlabel('Frequencey [Hz]')
plt.ylabel('Added mass')
plt.grid(True)
plt.legend()
plt.tight_layout()

plt.figure()
for dof in fb.dofs:
    plt.plot(
        f_range,
        data['radiation_damping'].sel(radiating_dof=dof, influenced_dof=dof),
        label=dof,
        marker='o',
    )
plt.xlabel('Frequencey [Hz]')
plt.ylabel('Radiation damping')
plt.grid(True)
plt.legend()
plt.tight_layout()

plt.show()

[mancellin]

I did not test it, but the example I posted above should also work for Capytaine:

# [...] define a mesh with pygmsh 

points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

def all_faces_as_tetra(cells):
    triangles_as_tetra = np.empty((cells['triangle'].shape[0], 4), dtype=np.int)
    triangles_as_tetra[:, :3] = cells['triangle'][:, :]
    triangles_as_tetra[:, 3] = cells['triangle'][:, 2]  # Repeat one node to make a tetra
    return np.concatenate([cells['tetra'], triangles_as_tetra])

mesh = cpt.Mesh(vertices=points, faces=all_faces_as_tetra(cells))
fb = cpt.FloatingBody(mesh=mesh, name="my_body")

[ryancoe]

This probably can work with some tweaking, but doesn't quite work yet. Unfortunately, I'm afraid this is far enough outside of my current knowledge base that it'd be another year before I can realistically contribute/address this error ;)

(using your snipped and the geometry from my post above)

Traceback (most recent call last):

  File "/Users/rcoe/mwePygmsh.py", line 41, in <module>
    mesh = cpt.Mesh(vertices=points, faces=all_faces_as_tetra(cells))

  File "/Users/rcoe/mwePygmsh.py", line 36, in all_faces_as_tetra
    triangles_as_tetra = np.empty((cells['triangle'].shape[0], 4), dtype=np.int)

TypeError: list indices must be integers or slices, not str

[mancellin]

Indeed, it seems that Pygmsh has changed since last year and returns a slightly different representation of the mesh. Also, my previous snippet might fail when there are only triangles or only quadrilaterals.

Below is an updated version that works in your case:

# mesh = pygmsh.generate_mesh(geom, ...)

def all_faces_as_tetra(cells):
    all_faces = []
    if 'tetra' in cells:
        all_faces.append(cells['tetra'])
    if 'triangle' in cells:
        triangles_as_tetra = np.empty((cells['triangle'].shape[0], 4), dtype=np.int)
        triangles_as_tetra[:, :3] = cells['triangle'][:, :]
        triangles_as_tetra[:, 3] = cells['triangle'][:, 2]  # Repeat one node to make a tetra
        all_faces.append(triangles_as_tetra)
    return np.concatenate(all_faces)

cpt_mesh = cpt.Mesh(vertices=mesh.points, faces=all_faces_as_tetra(mesh.cells_dict))
fb = cpt.FloatingBody(mesh=cpt_mesh, name="my_body")

And the same should work for meshmagick's Mesh object.

[ryancoe]

Very cool. Perhaps we can incorporate your all_faces_as_tetra method above into capytaine (perhaps as part of heal_mesh)?

[mancellin]

I just realized that what pygmsh calls a 'tetra' is probably a tetrahedron for 3D meshes and not a 2D quadrilateral as I'm used to. So the code snippet still need some polishing before we incorporate it.

Could you help me by doing some tests with several options of pygmsh to check if the snippet always works and if pygmsh can actually handle quadrilateral faces?

[ryancoe]

Sounds good. I messed with this a bit and wasn't able to get pygmsh to create quads in the same way that basic gmsh does... raised issue here https://github.com/nschloe/pygmsh/issues/331

[ryancoe]

@mancellin - Any preference on whether I put together tests for this via meshmagick or capytaine??

[mancellin]

As you wish. It should be fairly easy to convert from one to the other and I plan to add the pygmsh importer into both in the long run.

[saltynexus]

This is with regard to the original post. I'm using GMSH 4.6 and meshmagick 1.0.6 and it appears that meshmagick can not read ".msh" files. Here's my traceback

$ meshmagick t2.msh --show

============================================= meshmagick - version 1.0.6 Copyright 2014-2020, Ecole Centrale de Nantes

Traceback (most recent call last): File "anaconda3/bin/meshmagick", line 11, in sys.exit(main()) File "anaconda3/lib/python3.7/site-packages/meshmagick/init.py", line 5, in main mm.main() File "anaconda3/lib/python3.7/site-packages/meshmagick/meshmagick.py", line 752, in main V, F = mmio.load_mesh(args.infilename, format) File "anaconda3/lib/python3.7/site-packages/meshmagick/mmio.py", line 46, in load_mesh vertices, faces = loader(filename) File "anaconda3/lib/python3.7/site-packages/meshmagick/mmio.py", line 869, in load_MSH nb_nodes, nodes_data = re.search(r'\$Nodes\n(\d+)\n(.+)\$EndNodes', data, re.DOTALL).groups() AttributeError: 'NoneType' object has no attribute 'groups'

[mancellin]

@saltynexus Yes, the original issue of this thread is still there. Did you try to convert the mesh file to an older gmsh format using the -format option of gmsh? This is the only solution at the moment.