nschloe/pygmsh

:spider_web: Gmsh for Python

Gmsh for Python.

pygmsh combines the power of Gmsh with the versatility of Python. It provides useful abstractions from Gmsh's own Python interface so you can create complex geometries more easily.

To use, install Gmsh itself and pygmsh from pypi:

[sudo] apt install python3-gmsh
pip install pygmsh

This document and the tests/ directory contain many small examples. See here for the full documentation.

Flat shapes

Polygon	Circle	(B-)Splines

Codes:

import pygmsh

with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(
        [
            [0.0, 0.0],
            [1.0, -0.2],
            [1.1, 1.2],
            [0.1, 0.7],
        ],
        mesh_size=0.1,
    )
    mesh = geom.generate_mesh()

# mesh.points, mesh.cells, ...
# mesh.write("out.vtk")

import pygmsh

with pygmsh.geo.Geometry() as geom:
    geom.add_circle([0.0, 0.0], 1.0, mesh_size=0.2)
    mesh = geom.generate_mesh()

import pygmsh

with pygmsh.geo.Geometry() as geom:
    lcar = 0.1
    p1 = geom.add_point([0.0, 0.0], lcar)
    p2 = geom.add_point([1.0, 0.0], lcar)
    p3 = geom.add_point([1.0, 0.5], lcar)
    p4 = geom.add_point([1.0, 1.0], lcar)
    s1 = geom.add_bspline([p1, p2, p3, p4])

    p2 = geom.add_point([0.0, 1.0], lcar)
    p3 = geom.add_point([0.5, 1.0], lcar)
    s2 = geom.add_spline([p4, p3, p2, p1])

    ll = geom.add_curve_loop([s1, s2])
    pl = geom.add_plane_surface(ll)

    mesh = geom.generate_mesh()

The return value is always a meshio mesh, so to store it to a file you can

mesh.write("test.vtk")

The output file can be visualized with various tools, e.g., ParaView.

With

pygmsh.write("test.msh")

you can access Gmsh's native file writer.

Extrusions

extrude	revolve	twist

import pygmsh

with pygmsh.geo.Geometry() as geom:
    poly = geom.add_polygon(
        [
            [0.0, 0.0],
            [1.0, -0.2],
            [1.1, 1.2],
            [0.1, 0.7],
        ],
        mesh_size=0.1,
    )
    geom.extrude(poly, [0.0, 0.3, 1.0], num_layers=5)
    mesh = geom.generate_mesh()

from math import pi
import pygmsh

with pygmsh.geo.Geometry() as geom:
    poly = geom.add_polygon(
        [
            [0.0, 0.2, 0.0],
            [0.0, 1.2, 0.0],
            [0.0, 1.2, 1.0],
        ],
        mesh_size=0.1,
    )
    geom.revolve(poly, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0], 0.8 * pi)
    mesh = geom.generate_mesh()

from math import pi
import pygmsh

with pygmsh.geo.Geometry() as geom:
    poly = geom.add_polygon(
        [
            [+0.0, +0.5],
            [-0.1, +0.1],
            [-0.5, +0.0],
            [-0.1, -0.1],
            [+0.0, -0.5],
            [+0.1, -0.1],
            [+0.5, +0.0],
            [+0.1, +0.1],
        ],
        mesh_size=0.05,
    )

    geom.twist(
        poly,
        translation_axis=[0, 0, 1],
        rotation_axis=[0, 0, 1],
        point_on_axis=[0, 0, 0],
        angle=pi / 3,
    )

    mesh = geom.generate_mesh()

OpenCASCADE

Gmsh also supports OpenCASCADE (occ), allowing for a CAD-style geometry specification.

from math import pi, cos
import pygmsh

with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    r = 0.5
    disks = [
        geom.add_disk([-0.5 * cos(7 / 6 * pi), -0.25], 1.0),
        geom.add_disk([+0.5 * cos(7 / 6 * pi), -0.25], 1.0),
        geom.add_disk([0.0, 0.5], 1.0),
    ]
    geom.boolean_intersection(disks)

    mesh = geom.generate_mesh()

# ellpsoid with holes
import pygmsh

with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_max = 0.1
    ellipsoid = geom.add_ellipsoid([0.0, 0.0, 0.0], [1.0, 0.7, 0.5])

    cylinders = [
        geom.add_cylinder([-1.0, 0.0, 0.0], [2.0, 0.0, 0.0], 0.3),
        geom.add_cylinder([0.0, -1.0, 0.0], [0.0, 2.0, 0.0], 0.3),
        geom.add_cylinder([0.0, 0.0, -1.0], [0.0, 0.0, 2.0], 0.3),
    ]
    geom.boolean_difference(ellipsoid, geom.boolean_union(cylinders))

    mesh = geom.generate_mesh()

# puzzle piece
import pygmsh

with pygmsh.occ.Geometry() as geom:
    geom.characteristic_length_min = 0.1
    geom.characteristic_length_max = 0.1

    rectangle = geom.add_rectangle([-1.0, -1.0, 0.0], 2.0, 2.0)
    disk1 = geom.add_disk([-1.2, 0.0, 0.0], 0.5)
    disk2 = geom.add_disk([+1.2, 0.0, 0.0], 0.5)

    disk3 = geom.add_disk([0.0, -0.9, 0.0], 0.5)
    disk4 = geom.add_disk([0.0, +0.9, 0.0], 0.5)
    flat = geom.boolean_difference(
        geom.boolean_union([rectangle, disk1, disk2]),
        geom.boolean_union([disk3, disk4]),
    )

    geom.extrude(flat, [0, 0, 0.3])

    mesh = geom.generate_mesh()

Mesh refinement/boundary layers

# boundary refinement
import pygmsh

with pygmsh.geo.Geometry() as geom:
    poly = geom.add_polygon(
        [
            [0.0, 0.0],
            [2.0, 0.0],
            [3.0, 1.0],
            [1.0, 2.0],
            [0.0, 1.0],
        ],
        mesh_size=0.3,
    )

    field0 = geom.add_boundary_layer(
        edges_list=[poly.curves[0]],
        lcmin=0.05,
        lcmax=0.2,
        distmin=0.0,
        distmax=0.2,
    )
    field1 = geom.add_boundary_layer(
        nodes_list=[poly.points[2]],
        lcmin=0.05,
        lcmax=0.2,
        distmin=0.1,
        distmax=0.4,
    )
    geom.set_background_mesh([field0, field1], operator="Min")

    mesh = geom.generate_mesh()

# mesh refinement with callback
import pygmsh

with pygmsh.geo.Geometry() as geom:
    geom.add_polygon(
        [
            [-1.0, -1.0],
            [+1.0, -1.0],
            [+1.0, +1.0],
            [-1.0, +1.0],
        ]
    )
    geom.set_mesh_size_callback(
        lambda dim, tag, x, y, z: 6.0e-2 + 2.0e-1 * (x**2 + y**2)
    )

    mesh = geom.generate_mesh()

# ball with mesh refinement
from math import sqrt
import pygmsh


with pygmsh.occ.Geometry() as geom:
    geom.add_ball([0.0, 0.0, 0.0], 1.0)

    geom.set_mesh_size_callback(
        lambda dim, tag, x, y, z: abs(sqrt(x**2 + y**2 + z**2) - 0.5) + 0.1
    )
    mesh = geom.generate_mesh()

Optimization

pygmsh can optimize existing meshes, too.

import meshio

mesh = meshio.read("mymesh.vtk")
optimized_mesh = pygmsh.optimize(mesh, method="")

You can also use the command-line utility

pygmsh-optimize input.vtk output.xdmf

where input and output can be any format supported by meshio.

Testing

To run the pygmsh unit tests, check out this repository and type

pytest

Building Documentation

Docs are built using Sphinx.

To build, run

sphinx-build -b html doc doc/_build

License

This software is published under the GPLv3 license.

pygmsh's Issues

3D boolean operation problem

It seems 3D boolean operations don't work properly. Try meshing the geom created by:

def generate():
    geom = pg.Geometry()
    geom.set_factory('OpenCASCADE')
    lcar = 0.1
    ball=geom.add_ball([0.0, 0.0, 0.0], 0.5, lcar)
    box=geom.add_box(-1, 1, -1, 1, -1, 1, lcar)
    geom.boolean_difference([box.volume], [ball.volume])
    return geom

Or have I done something wrong?

How to menage a mesh?

Dear pygmsh users, can you please give me examples (links) how to manage cells inside of mesh - with their size, form,dimensions,surface and amount ? I want to set size, intervals and other things in python code by myself

generate_mesh raises two RuntimeWarnings for one-dimensional geometries

The function helper._sit_in_plane added in #50 raises two RuntimeWarnings, first on trying to divide by the zero value of the squared-magnitude of the cross-product of two collinear vectors, second trying to compare the resulting quotient which is NaN to tol.
To reproduce:

g = Geometry()
g.add_line(*(g.add_point((x, 0, 0), 1) for x in [0, 2]))
print(generate_mesh(g))

The code runs fine and does generate correct one-dimensional meshes, but I wonder whether there isn't a yet more elegant test for planarity that would avoid these warnings.

Ruled Surfaces

The lack of a facility for emitting Gmsh.geo code for Ruled Surfaces arose in #121.

accept expression-lists for Physical entities

the expression-list on the right hand side should contain the identification numbers of all the elementary points that need to be grouped inside the physical point. [http://gmsh.info/doc/texinfo/gmsh.html#index-Physical-points]

the expression-list on the right hand side should contain the identification numbers of all the elementary lines that need to be grouped inside the physical line [http://gmsh.info/doc/texinfo/gmsh.html#index-Physical-lines]

Usage is widespread; e.g. https://onelab.info/svn/gmsh/trunk/tutorial/t1.geo

Physical Point(1) = {1,2} ;

I hope to propose a pull-request shortly.

Examples - passing tuples where tuples aren't expected

In the test/examples/ directory the files with lines out = pg.generate_mesh(generate()) should to be changed to out = generate_mesh(generate()[0]) or they don't run

installation

I just want to verify, but do you need to have gmsh installed first. The documentation makes no mention of installing gmsh.

Dimension header when saving a 2D mesh

Hi Nico,

Thanks for creating such useful tools as pygmsh and meshio. I have been using them for a few days now, and it has been a big help. I did run into a small problem when I tried to create a 2D mesh. I understand that gmsh requires 3D points, and the "dim" variable in the generate_mesh function here is passed correctly to gmsh. However, if I save the output geometry as a .mesh file, it always seems to list "Dimension 3" regardless of "dim".

If I first manually adapt the array containing the points and then save the mesh meshio does correctly store it as "Dimension 2":

out = pg.generate_mesh(geom,dim=2)

points=out[0][:,0:2]
cells=out[1]

meshio.write(
    'test_forced_2D.mesh',
    points,
    cells,
    )

However, this doesn't seem very elegant. Am I doing something wrong, or is there maybe a better way of ensuring that a 2D mesh is saved with the correct value for the Dimension?

Best regards,
Ruben

Structured mesh of a spherical lune

Is it possible to create a structured mesh of a spherical lune using pygmsh ?

triangular/quadrilateral mesh in pygmsh

Hi Nico,

I believe that you are the only person who can help me with my master dissertation at this point .We have had conversation on github and you have helped me a lot. I am sending you my python code and I have 2 questions:

1. How can I make a triangular mesh with my code(or any code where I have to read data from file and where I only know which points are creating splines,line loops and surface)? Is it possible to make the quadrilateral mesh in pygmsh? I've manually put all data from geo file in geometry like you said I needed to do (#121). Do you have some links or examples where you did it manually and made some mesh with that data?

2)How can I manipulate and remember only points on the surface of geometric character and the points which create triangular mesh? I will appreciate all links and examples you have on this too.

Thank you so much for your work and help

reconcile quad meshes with optimize-lloyd

See http://onelab.info/pipermail/gmsh/2017/011814.html.

import stl

Hi, I am trying to import an STL file and I get an error:

% python
Python 3.5.3 |Intel Corporation| (default, Apr 27 2017, 18:08:47)
[GCC 4.8.2 20140120 (Red Hat 4.8.2-15)] on linux
Type "help", "copyright", "credits" or "license" for more information.
Intel(R) Distribution for Python is brought to you by Intel Corporation.
Please check out: https://software.intel.com/en-us/python-distribution
>>> import meshio as io
>>> io.read('./data/03_FULL_DELIVERY_22/VHP22_STL_FILES/VHPC_Heart.stl',file_format='stl')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/opt/intel/intelpython3/lib/python3.5/site-packages/meshio/helpers.py", line 121, in read
    return format_to_reader[file_format].read(filename)
KeyError: 'stl'

Any ideas?

Adding Polygons or Rectangles that share common points and lines

Hi Nico

The add_polygon or add_rectangle commands are very useful in that points, lines, and line loops are all created for you and so do not complicate the code. However, if you add polygons or rectangles that share common points and lines new points and lines will still be created with the same coordinates, and Gmsh does not seem to like that. It would be useful I think if an option could be added to these commands in which a list of existing polygons/rectangles could be added as an argument so that the command could then identify existing points and lines (within a tolerance bound) and not generate new ones.

Imagine how much simpler the following code would look without the add_point, add_line, add_line_loop commands

for i in range(numfins):
    points[8*i] = geom.add_point([i*(fingap+finbase), -base, 0.0],edge)
    points[8*i+1] = geom.add_point([i*(fingap+finbase)+finbase/2.0-fintip/2.0, -base-finlength, 0.0],edge)
    points[8*i+2] = geom.add_point([i*(fingap+finbase)+finbase/2.0-fintip/2.0, -base-finlength, -length],edge)
    points[8*i+3] = geom.add_point([i*(fingap+finbase), -base, -length],edge)
    lines[8*i] = geom.add_line(points[8*i],points[8*i+1])
    lines[8*i+1] = geom.add_line(points[8*i+1],points[8*i+2])
    lines[8*i+2] = geom.add_line(points[8*i+2],points[8*i+3])
    lines[8*i+3] = geom.add_line(points[8*i+3],points[8*i])
    loops[2*i] = geom.add_line_loop({lines[8*i],lines[8*i+1],lines[8*i+2],lines[8*i+3]})
    surfaces[2*i] = geom.add_plane_surface(loops[2*i])

accept expressions for Physical entity identifiers

This picks up a tangent from #19 . That was originally about the expression-list in the right-hand side of the constructor for a Physical entity but @beckerrh also raised there a related issue about the left-hand side only accepting a char-expression whereas according to e.g. §5.5.1 Points the syntax is

Physical Point ( expression | char-expression <, expression> ) <+|->= { expression-list };

i.e. it should also accept an expression, where 'The expression inside the parentheses is the physical point’s identification number'.

This is trivially exemplified in the same snippet from t1.geo that I used to open #19:

Physical Point(1) = {1,2} ;

So the answer to the query of @nschloe

@beckerrh Is this really supported? I only ever see strings as Physical labels.

is yes, numerically valued non-char-expression expressions are supported for Physical labels in Gmsh and maybe should be in pygmsh to.

-optimize -optimize_lloyd 10 and gmsh 3.0.1

import pygmsh as pg
import numpy as np

geom = pg.Geometry()

# Draw a cross.
poly = geom.add_polygon([
    [0.0,   0.5, 0.0],
    [-0.1,  0.1, 0.0],
    [-0.5,  0.0, 0.0],
    [-0.1, -0.1, 0.0],
    [0.0,  -0.5, 0.0],
    [0.1,  -0.1, 0.0],
    [0.5,   0.0, 0.0],
    [0.1,   0.1, 0.0]
    ],
    lcar=0.05
    )

axis = [0, 0, 1]

geom.extrude(
    poly,
    translation_axis=axis,
    rotation_axis=axis,
    point_on_axis=[0, 0, 0],
    angle=2.0 / 6.0 * np.pi
    )

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

Warning : The '-optimize' option is now obsolete: Gmsh optimizes tetrahedral meshes by default
Warning : Use '-optimize_threshold threshold' to control which elements are optimized
Warning : Option '-optimize_threshold 0' leads to no optimization
RuntimeError: Gmsh exited with error (return code -4).

Generating a Square Directly Using geom.add_polygon

Hi there, and thanks for making Gmsh available to us through Python!

I am new to both Gmsh and PyGmsh, and am trying to use these to create a simple mesh for use with libMesh.

I'll begin with the simplest case: while similar code is working to generate both a disk and the extruded star (example) the following:

Question 1
"Square Example"

geom = pygmsh.built_in.Geometry()

pts = [[x,y,z] for x in [-1.,1.] for y in [-1.,1.] for z in [0.]];

Draw a cross.

poly = geom.add_polygon(pts,lcar=0.05)
points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom,optimize=False, num_quad_lloyd_steps=10, num_lloyd_steps=10, verbose=True, dim=3, prune_vertices=True, gmsh_path=None)

import meshio
meshio.write('cube.vtu', points, cells, cell_data=cell_data)

results in an error:

Info    : Running 'gmsh -3 -bin /tmp/tmprgpnttr6.geo -o /tmp/tmpz12edeie.msh -optimize_lloyd 10' [Gmsh 2.10.1, 1 node, max. 1 thread]
Info    : Started on Wed Nov 22 18:48:18 2017
Info    : Reading '/tmp/tmprgpnttr6.geo'...
Info    : Done reading '/tmp/tmprgpnttr6.geo'
Info    : Meshing 1D...
Info    : Meshing curve 1 (Line)
Info    : Meshing curve 2 (Line)
Info    : Meshing curve 3 (Line)
Info    : Meshing curve 4 (Line)
Info    : Done meshing 1D (0 s)
Info    : Meshing 2D...
Info    : Meshing surface 6 (Plane, Delaunay)
Warning : :-( There are 2 intersections in the 1D mesh (curves 2 4)
Warning : 8-| Gmsh splits those edges and tries again
Error   : Unable to recover an edge 0.0175439 0.0175439 && 2.72926e-12 2.72926e-12 (28/57)
Info    : Done meshing 2D (0.00740314 s)
Info    : Meshing 3D...
Info    : Done meshing 3D (0 s)
Info    : 196 vertices 200 elements
Error   : ------------------------------
Error   : Mesh generation error summary
Error   :     2 warnings
Error   :     1 error
Error   : Check the full log for details
Error   : ------------------------------
Info    : Writing '/tmp/tmpz12edeie.msh'...
Info    : Done writing '/tmp/tmpz12edeie.msh'
Info    : Stopped on Wed Nov 22 18:48:18 2017
---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-33-5cdd6c88155a> in <module>()
     16 #     )
     17 
---> 18 points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom,optimize=False, num_quad_lloyd_steps=10, num_lloyd_steps=10, verbose=True, dim=3, prune_vertices=True, gmsh_path=None)
     19 
     20 import meshio

~/anaconda3/lib/python3.6/site-packages/pygmsh/helpers.py in generate_mesh(geo_object, optimize, num_quad_lloyd_steps, num_lloyd_steps, verbose, dim, prune_vertices, gmsh_path)
    129     p.communicate()
    130     assert p.returncode == 0, \
--> 131         'Gmsh exited with error (return code {}).'.format(p.returncode)
    132 
    133     X, cells, pt_data, cell_data, field_data = meshio.read(msh_filename)

AssertionError: Gmsh exited with error (return code 1).

Changing the number of lloyd steps and toggling optimization and prune vertices does not seem to resolve this.

Question 2
Also, as a simpler question: do I need OpenCascade in order to use this in a CAD-like style?

Thanks, and Happy Thanksgiving
Charlie

Cannot open stp file in python

Dear pygmsh and python users, I have connected pygmsh, vtk files,paraview , python and gmsh. My question is how to connect .stp file with python. Inside of that file I have geometry which I have to use to create mesh. I do not know what all text in .stp file mean so I have not tried to read line by line from .stp file yet.

Voropy throws AssertionError

Hi pygmsh folks,

I want to create a 2D mesh from raw gmsh code which works fine, if i do it directly in gmsh. Since i want to create it dynamically, i'd like to use pygmsh for that. Unfortunately I always get an AssertionError from Voropy:

assert all(abs(_row_dot(e2, e2) - _row_dot(e1, e1)) < 1.0e-14)

I think it would be nice if there was a switch in the pygmsh.generate_mesh function to enable or disable Lloyd smoothing.

But nevertheless: great work! ;)

Adding Physical Groups

Hi

I am not sure if I am applying the add_physical_volume command correctly. For instance in the given extruded and twisted 4 point star example the extruded volume vol gets assigned the string ex1 when I run the python script. But when I add

geom.add_physical_volume(vol, label='twisted_star')

no physical group volume in Gmsh appears because the volume is assigned in Gmsh the label name (p1) and not ex1. When I do

geom.add_physical_volume('p1', label='twisted_star')

that works.

I also note that when I inspect the Elementary Entities in Gmsh some lines are assigned names starting with p. I started writing a python script of my own using pygmsh and saw later in Gmsh some line names starting with 'l' and some with 'p' . I though 'l' was for lines and 'p' for points. In the twisted star example Line1 (first assigned line of the polygon) also has the name (p1) in Gmsh. In the .geo text file produced by Python it is 'l1' as it might be expected. Not sure if this is a bug in Gmsh, or my understanding is flawed. I am running Gmsh 2.16.0 on OS X 10.11.6. In any case consistent naming is important because, for example, the geom.add_polygon command creates a host of points, lines and surfaces that one might want to refer to later by 'p1', 'l1', 'surf1' etc., as more polygons get added.

Screenshot:

Could you possibly add to the twisted star Python script example say a physical group name to the bottom and/or top star surfaces and also to the volume, and check in Gmsh that the physical groups appear as expected? I am using Elmer as the solver and assignment of the boundary conditions is made easier when the surfaces are defined as physical groups.

With thanks and appreciation for your efforts

holes no longer respected in extrusion

Geometry.extrude of a surface with a hole, as in test_layers.py or test_screw.py, doesn't quite respect the hole anymore. The volume is toroidal but the far surface is simply connected.

This isn't caught by the test using compute_volume.

Python 3.5 support ?

Hello,

When I try to pip install pygmsh (MacOS 10.10.5) that yields the following error:
ImportError: cannot import name 'HTTPSHandler'

I didn't find information on python version that you support, is it supposed to work on Python 3.x ?

Spline not BSpline for built_in

Often an interpolating spline is desired. In Gmsh, with the default built-in factor, this is provided by Spline rather than BSpline; however, pygmsh/built_in currently only provides the latter.

As the .geo syntax is identical, this should be a straightforward addition.

This is a counterpart to #91 ‘Spline (not BSpline) for OpenCascade’.

Gmsh_executable under MacOs

Gmsh could not be installed in /Applications, for example for brew-installed Gmsh.
Maybe we could have a check in 'helper.py':

if sys.platform == 'darwin':
    # likely there.
    gmsh_executable = '/Applications/Gmsh.app/Contents/MacOS/gmsh'
else:
    gmsh_executable = 'gmsh'

Issues when running using MPI

I use pygmsh at the beginning of some finite element code which is run in parallel using MPI. Since upgrading from 3.0.15 to the current version, the code aborts as soon as pygmsh.Geometry() is called.
MPI Errorcode is 59, unfortunately that doesn't help me very much.

Any idea what this might be? I'm running Ubuntu 17.04 with OpenMPI if that helps.

Thanks!

Set path to gmsh executable

The example didn't run right away as here the gmsh executable path is hardcoded for MacOS. Maybe add a seperate function that allows you to set the correct path manually?

pygmsh/pygmsh/helper.py

Line 70 in fb10b7b

macos_gmsh_location = '/Applications/Gmsh.app/Contents/MacOS/gmsh'

use pygmsh for mesh

Does pygmsh support creation of geometry only or does it support mesh commands of gmsh as well? I'm looking for something like the attached geo file.

Merge "tetra-02.brep";


///////////////////////////////////////////////////////
//Geometry Options, "= 0"  is False!
Geometry.OCCFixDegenerated = 0;
// Fix degenerated edges/faces in STEP, IGES and BRep models

Geometry.OCCFixSmallEdges = 0;
// Fix small edges in STEP, IGES and BRep models

Geometry.OCCFixSmallFaces = 0;
// Fix small faces in STEP, IGES and BRep models

Geometry.OCCSewFaces = 0;
// connect edges --> for shells 

Geometry.OCCConnectFaces = 0;
// connect faces --> for solids


///////////////////////////////////////////////////////
// elementsize (watch your units!)
//Mesh.CharacteristicLengthMax = 0.5;
Mesh.CharacteristicLengthMax = 1e+22;
//Default value: 1e+22

Mesh.CharacteristicLengthMin = 20;
// Minimum mesh element size



// optimization
//Mesh.Optimize = 0;
// Optimize the mesh to improve the quality of tetrahedral elements
//Mesh.OptimizeNetgen = 0;
// Optimize the mesh using Netgen to improve the quality of tetrahedral elements

//Mesh.HighOrderOptimize = 0;
// Optimize high order meshes?



///////////////////////////////////////////////////////
Mesh.ElementOrder = 2;
// Element order (1=linear elements, N (<6) = elements of higher order)
// Default value: 1

Mesh.Algorithm3D = 5;
// 3D mesh algorithm: 1=Delaunay, 4=Frontal, 5=Frontal Delaunay, 6=Frontal Hex, 7=MMG3D, 9=R-tree
// default = 1


///////////////////////////////////////////////////////
// meshen immediately after loading the file into GUI --> not needed in konsole
//Mesh  2;  // surface
//Mesh  3;  // volume


///////////////////////////////////////////////////////
Mesh.Format = 33;
// 1=msh, 2=unv, 10=automatic, 19=vrml, 27=stl, 30=mesh, 31=bdf, 32=cgns, 33=med, 40=ply2
// default = 10

Mesh.SaveAll = 1;
// Ignore Physical definitions and save all elements

//Save "tetra-03.unv";
//Save "tetra-03.msh";
Save "tetra-03.med";

meshio.write

I tried to store my mesh, using meshio.write, but it print an assertion error:

AssertionError: Only contiguous arrays are supported.

Kindly assist me.
-Thanks

Code used:

import meshio
meshio.write('test.vtu', points, cells, cell_data=cell_data)

Error:

AssertionError                            Traceback (most recent call last)
<ipython-input-1-2defc9f31bea> in <module>()
     12     import meshio
     13     out = pg.generate_mesh(generate()[0])
---> 14     meshio.write('cube.vtu', *out)
     15 

/Users/home/grahake/anaconda2/lib/python2.7/site-packages/meshio/helpers.pyc in write(filename, points, cells, point_data, cell_data, field_data, file_format)
    185             point_data=point_data,
    186             cell_data=cell_data,
--> 187             field_data=field_data
    188             )
    189     elif file_format == 'vtk-ascii':

/Users/home/grahake/anaconda2/lib/python2.7/site-packages/meshio/vtk_io.pyc in write(filetype, filename, points, cells, point_data, cell_data, field_data)
    234                     'Cell data mismatch.'
    235 
--> 236     vtk_mesh = _generate_vtk_mesh(points, cells)
    237 
    238     # add point data

/Users/home/grahake/anaconda2/lib/python2.7/site-packages/meshio/vtk_io.pyc in _generate_vtk_mesh(points, cells)
    333     vtk_points = vtk.vtkPoints()
    334     # Not using a deep copy here results in a segfault.
--> 335     vtk_array = numpy_support.numpy_to_vtk(points, deep=True)
    336     vtk_points.SetData(vtk_array)
    337     mesh.SetPoints(vtk_points)

/Users/home/grahake/anaconda2/lib/python2.7/site-packages/vtk/util/numpy_support.pyc in numpy_to_vtk(num_array, deep, array_type)
    129 
    130     shape = z.shape
--> 131     assert z.flags.contiguous, 'Only contiguous arrays are supported.'
    132     assert len(shape) < 3, \
    133            "Only arrays of dimensionality 2 or lower are allowed!"

AssertionError: Only contiguous arrays are supported.

unicode labels for physical groups broken

My code was broken I think by #47. My code attempts to give unicode labels to physical entities. Since pipupdating today, geom.add_physical_point now raises an AssertionError at assert isinstance(label, str).

Short example:

geom = Geometry()
point = geom.add_point([0, 0, 0], 1)
geom.add_physical_point(point, u'origin')

without the u, it works as before.

I think it should be a matter of accepting unicode or str rather than insisting on the latter?

Is Gmsh version 3 required to use opencascade features?

Hi again,

First, quickly: do we need gmsh version 3 to use pygmsh?

I ran pytest and received:
============================= test session starts ==============================
platform linux -- Python 3.5.2, pytest-3.3.1, py-1.5.2, pluggy-0.6.0
rootdir: /home/chaztikov/anaconda3/pkgs/pygmsh, inifile:
collected 32 items

test/test_airfoil.py F [ 3%]
test/test_boundary_layers.py F [ 6%]
test/test_bsplines.py F [ 9%]
test/test_circle.py F [ 12%]
test/test_cube.py *** Error in `gmsh': free(): invalid pointer: 0x00007fba2f71a348 ***
======= Backtrace: =========

Not sure if this is an issue with my gmsh install or pygmsh, but I was wondering if you might have some insight?

My original question:

As far as I know I have opencascade and gmsh (installed from Synaptic Package Manager) and pygmsh (installed through anaconda I believe) and I get this error not only when trying to use the opencascade feature, but also when I try to create something simple, say, an ellipsoid without opencascade features.

In particular, when I run test_opencascade_ball.py I get

Info : Running 'gmsh -3 -bin /tmp/tmp50u0hdt8.geo -o /tmp/tmp4pnhhetn.msh -optimize -optimize_lloyd 10' [Gmsh 2.10.1, 1 node, max. 1 thread]
Info : Started on Fri Dec 15 20:58:53 2017
Info : Reading '/tmp/tmp50u0hdt8.geo'...
Error : '/tmp/tmp50u0hdt8.geo', line 2 : syntax error (")
Error : '/tmp/tmp50u0hdt8.geo', line 4 : Sphere 1 has to be defined using 2 points (center + any point) and not 7
Warning : '/tmp/tmp50u0hdt8.geo', line 5 : Unknown volume 1
Error : '/tmp/tmp50u0hdt8.geo', line 5 : Unknown command on multiple shapes: 'PointsOf'
Info : Done reading '/tmp/tmp50u0hdt8.geo'
Info : Meshing 1D...
Info : Done meshing 1D (0 s)
Info : Meshing 2D...
Info : Done meshing 2D (5.96046e-06 s)
Info : Meshing 3D...
Info : Done meshing 3D (0 s)
Info : 0 vertices 0 elements
Info : Writing '/tmp/tmp4pnhhetn.msh'...
Info : Done writing '/tmp/tmp4pnhhetn.msh'
Info : Stopped on Fri Dec 15 20:58:53 2017
Traceback (most recent call last):
File "test_opencascade_ball.py", line 31, in
meshio.write('opencascade_ball.vtu', *test())
File "test_opencascade_ball.py", line 24, in test
points, cells, _, _, _ = pygmsh.generate_mesh(geom)
File "/usr/local/lib/python3.5/dist-packages/pygmsh/helpers.py", line 131, in generate_mesh
'Gmsh exited with error (return code {}).'.format(p.returncode)
AssertionError: Gmsh exited with error (return code 1).

Insight is appreciated, thanks!

Update: it seems there is a greater issue at hand,
pytest result :
======================================================================= 18 failed, 1 passed, 13 skipped, 1 warnings in 66.57 seconds =======================================================================

built_in spline faulty assert

In built_in/spline.py:13, a check is made more than 3 points were passed as input. However, gmsh's spline works also with 3 (or even 2) points.

unrotated extrusion-layers no longer respected

I think #126 has broken the case of Geometry.extrude with rotation_axis None when num_layers is not None; see comment on 666b6f8.

This is the case of test_layers.py #109. Running that now gives unstructured extrusion without layers.

Test error

I tried the the test code that you posted. Importing pygmsh is fine but the: geom.extrude raised an error:

geom.extrude(
    poly,
    translation_axis=axis,
    rotation_axis=axis,
    point_on_axis=[0, 0, 0],
    angle=2.0 / 6.0 * np.pi
    )

Out[20]: (<pygmsh.dummy.Dummy at 0xb028080>, <pygmsh.dummy.Dummy at 0xb0284e0>)

I'm using spyder python 3.5, Windows 64bits.
Please help me to solve this issue!

Support of (truncated) cones

https://en.wikipedia.org/wiki/Frustum

Translational extrusion of structured grids

As in the example t3.geo, we would like to be able to have pygmsh created extruded structured grids, eg. assuming you've created a sensible line-loop called ll0:

s0 = news;
Plane Surface(s0) = {ll0};
ex1[] = Extrude{10,0,0}{Surface{s0}; Layers{1}; Recombine;};

This would be useful for CFD codes which require 3d meshes which are one element thick to run 2d simulations (eg OpenFOAM).

This could be done by extending the extrude method:

1. change way function is called:

    def extrude(
            self,
            input_entity,
            translation_axis=None,
            rotation_axis=None,
            point_on_axis=None,
            angle=None,
            layers=None,
            recombine=False,
            ):

2. then in the part of the code where the translation only part is described:

       elif translation_axis is not None:
           # Only translation
           if layers is not None:
               self._GMSH_CODE.append(
                   '{}[] = Extrude{{{}}}{{{}; Layers{{{}}}; {}}};'.format(
                       name,
                       ','.join(repr(x) for x in translation_axis),
                       entity.id,
                       layers,
                       'Recombine;' if recombine else ''
                   ))
           else:
               self._GMSH_CODE.append(
                   '{}[] = Extrude{{{}}}{{{};}};'.format(
                       name,
                       ','.join(repr(x) for x in translation_axis),
                       entity.id
                   ))

3. in your python script:

   surface = geom.add_plane_surface(line_loop)
   axis = [10, 0, 0]
   top, volume, lat = geom.extrude(surface, translation_axis=axis, layers=1, recombine=True)

Ellipsoid

A function for an ellipsoid would be appreciated.

compound surface/volume from union

Wait for gmsh 3.1+ to get a compund surface from unions. See http://onelab.info/pipermail/gmsh/2017/011794.html.

ValueError: invalid literal for int() with base 10: '2-git'

line 93, in generate_mesh
gmsh_version = _get_gmsh_version()
line 71, in _get_gmsh_version
    return [int(x) for x in ex]
ValueError: invalid literal for int() with base 10: '2-git'

gmsh -version
--------------------------------------------------------------------------
Petsc Release Version 3.7.5, Jan, 01, 2017 
       The PETSc Team
    [email protected]
 http://www.mcs.anl.gov/petsc/
See docs/changes/index.html for recent updates.
See docs/faq.html for problems.
See docs/manualpages/index.html for help. 
Libraries linked from /home/geuzaine/src/petsc-3.7.5/complex_mumps_seq/lib
--------------------------------------------------------------------------
3.0.2-git

meshio in README demo doesn't work

Environment
Windows 10 :: Python 3.6.1 :: Anaconda 4.4.0 (64-bit)
pygmsh (3.1.0)
meshio (1.8.9)

Problem Desc

• pg.generate_mesh(geom) works well
• meshio.write('test.vtu', points, cells, cell_data=cell_data) returns with exit code 1

Debug Info

Traceback (most recent call last):
  File "D:/Documents/PycharmProjects/demo_pygmsh/demo_pygmsh.py", line 36, in <module>
    meshio.write('test.vtu', points, cells, cell_data=cell_data)
  File "D:\Documents\PycharmProjects\demo_pygmsh\meshio\helpers.py", line 187, in write
    field_data=field_data
  File "D:\Documents\PycharmProjects\demo_pygmsh\meshio\vtk_io.py", line 236, in write
    vtk_mesh = _generate_vtk_mesh(points, cells)
  File "D:\Documents\PycharmProjects\demo_pygmsh\meshio\vtk_io.py", line 335, in _generate_vtk_mesh
    vtk_array = numpy_support.numpy_to_vtk(points, deep=True)
  File "D:\Anaconda3\Lib\site-packages\vtk\util\numpy_support.py", line 131, in numpy_to_vtk
    assert z.flags.contiguous, 'Only contiguous arrays are supported.'
AssertionError: Only contiguous arrays are supported.

Process finished with exit code 1

Getting structured square cells in mesh

Hi,

this is my python code (the important part):
`with open('proba.txt', 'r') as myfile:
for line in myfile:
if 'Surface' in line:

		Surface.append(find_id(line))
		geom.add_raw_code(line)
	
	else:
		geom.add_raw_code(line)

myfile.close()

geom.add_raw_code('Mesh.Algorithm=8;')

for i in Surface:

geom.add_raw_code('Recombine Surface {%s}=0;' %i)

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

meshio.write('test8.vtu', points, cells, cell_data=cell_data)

`
And the pictures of mesh I get:

The problem is that I have to have all squares or rectangles in my mesh and I am not able to get squares and rectangles at some surfaces(picture above). Do you have any advice how to repare that, which (py)gmsh command to use?

Error with pygmsh library

Good morning,
I have a problem that causes me error with the pygmsh library, at the moment of executing code for that library.
I describe the steps I followed:
Install python 3.6.4
Install pygmsh 4.1.1
Execute Gmsh version 3.0.6
charge in the IDLE 3.6.4 (64 bits) a test code:

# -*- coding: utf-8 -*-
"""
Created on Thu Feb 15 12:11:41 2018
@author: Me
"""
import pygmsh
import numpy as np
import meshio
geom = pygmsh.built_in.Geometry()
ballHole = geom.add_ball(x0=[0.5,0.5,0.5], radius=0.25, lcar=0.10,)
geom.add_box(0, 1, 0, 1, 0, 1, 0.10, holes=[ballHole.surface_loop])
points, cells, point_data, cell_data, field_data = pygmsh.helpers.generate_mesh(geom)
meshio.write('test.msh', points, cells, cell_data=cell_data) # The .msh file in gmsh format
meshio.write('test.vtu', points, cells, cell_data=cell_data) #The output file can be visualized with various tools, e.g., ParaView.

I run the test code
and I get the following error:

>>> 
 RESTART: I:\Users\mramos\OneDrive\Documentos\IMP\Manuel\Python\pygmsh_test.py 
Traceback (most recent call last):
  File "I:\Users\mramos\OneDrive\Documentos\IMP\Manuel\Python\pygmsh_test.py", line 12, in <module>
    points, cells, point_data, cell_data, field_data = pygmsh.helpers.generate_mesh(geom)
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\site-packages\pygmsh\helpers.py", line 119, in generate_mesh
    gmsh_major_version = get_gmsh_major_version(gmsh_executable)
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\site-packages\pygmsh\helpers.py", line 82, in get_gmsh_major_version
    stderr=subprocess.STDOUT
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\subprocess.py", line 336, in check_output
    **kwargs).stdout
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\subprocess.py", line 403, in run
    with Popen(*popenargs, **kwargs) as process:
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\subprocess.py", line 709, in __init__
    restore_signals, start_new_session)
  File "C:\Users\mramos\AppData\Local\Programs\Python\Python36\lib\subprocess.py", line 997, in _execute_child
    startupinfo)
FileNotFoundError: [WinError 2] El sistema no puede encontrar el archivo especificado
>>> 

I appreciate your help in this regard.

best regards

Support of hyperboloids

https://en.wikipedia.org/wiki/Hyperboloid

FileNotFoundError: [WinError 2] The system cannot find the file specified

Running first code block example in ReadMe yields following error:

Traceback (most recent call last):
  File "C:/Users/Sterling Butters/PycharmProjects/Monte Carlo Neutron Transport/PyGMSH", line 31, in <module>
    points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)
  File "C:\ProgramData\Anaconda3\lib\site-packages\pygmsh\helpers.py", line 120, in generate_mesh
    gmsh_major_version = get_gmsh_major_version(gmsh_executable)
  File "C:\ProgramData\Anaconda3\lib\site-packages\pygmsh\helpers.py", line 82, in get_gmsh_major_version
    stderr=subprocess.STDOUT
  File "C:\ProgramData\Anaconda3\lib\subprocess.py", line 336, in check_output
    **kwargs).stdout
  File "C:\ProgramData\Anaconda3\lib\subprocess.py", line 403, in run
    with Popen(*popenargs, **kwargs) as process:
  File "C:\ProgramData\Anaconda3\lib\subprocess.py", line 709, in __init__
    restore_signals, start_new_session)
  File "C:\ProgramData\Anaconda3\lib\subprocess.py", line 997, in _execute_child
    startupinfo)
FileNotFoundError: [WinError 2] The system cannot find the file specified

Process finished with exit code 1

FileNotFoundError: [WinError 2] The system cannot find the file specified

Dear pygmsh users,

I have just tried to run the first example on the pygmsh website, but I am getting the following error:

FileNotFoundError: [WinError 2] The system cannot find the file specified

I am using Anaconda3 and running the code via Spyder.

Could you please advise how I can fix this issue?

The code I am trying to run is as follows:

import pygmsh
import numpy as np

geom = pygmsh.built_in.Geometry()

# Draw a cross.
poly = geom.add_polygon([
    [0.0,   0.5, 0.0],
    [-0.1,  0.1, 0.0],
    [-0.5,  0.0, 0.0],
    [-0.1, -0.1, 0.0],
    [0.0,  -0.5, 0.0],
    [0.1,  -0.1, 0.0],
    [0.5,   0.0, 0.0],
    [0.1,   0.1, 0.0]
    ],
    lcar=0.05
    )

axis = [0, 0, 1]

geom.extrude(
    poly,
    translation_axis=axis,
    rotation_axis=axis,
    point_on_axis=[0, 0, 0],
    angle=2.0 / 6.0 * np.pi
    )

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

Thanks!

gmsh 3 features

Would you accept patches that work only with gmsh 3?
I am thinking about OpenCASCADE boolean operations.
Of course those would be properly guarded to avoid errors with gmsh 2.

IGES

Not a real issue, I just wanted to answer a query about IGES that had been appended irrelevantly to #108.

From @tady57, a few hours ago

My question is how to connect .igs file with python. Inside of that file I have geometry which I have to use to create mesh. I do not know what all text in .igs file mean so I have not tried to read line by line from .igs file yet.

I expect .igs is an MS-DOS-style three-letter suffix for IGES, a CAD-format largely superseded by STEP. The latter was discussed and dismissed in #121 ; much of that applies to IGES too: it can be opened directly in Gmsh with

Merge "thing.igs";

but not much done with it thereafter, so it doesn't have much to do with pygmsh, I wouldn't think.

New versions of gmsh don't support -optimize_lloyd

New versions of gmsh have removed the "-optimize_lloyd" option, which breaks the default behavior of pygmsh (e.g. the two examples in README.md don't run without modification).

Recombine Surface

I tried to create a quadrilateral two-dimensional mesh by inserting

Recombine Surface {s0};

into the Gmsh .geo code to be emitted, but it didn't seem to work.

Specifically, I started a branch recsurf and in test copied and modified test_rectangle.py to test_quads.py by:

• catching the Polygon returned by add_rectangle
• putting its surface.id into Recombine Surface using add_raw_code (for the moment, with the idea of maybe adding a specialized method later)
• saving the output of get_code

    rectangle = geom.add_rectangle(...)
    geom.add_raw_code('Recombine Surface {%s};' % rectangle.surface.id)
    with open('quads.geo', 'w') as fgeo:
        fgeo.write(geom.get_code())

Then on comparing the meshes generated by

• pygmsh.generate_mesh and meshio.write, as in test_rectangle.py and
• running Gmsh manually on the quads.geo file

I find that the latter has quadrilaterals but the former only triangles.

Is the directive to Recombine Surface ignored or lost in going via meshio.read and meshio.write?

Is there a more direct way to get unstructured quadrilateral surface meshes?

(One application of this is to combine it with the extrude(layers=...) option added in #109 to get hexahedral meshes on extrusions of fairly arbitrary surfaces.)

Need ability to reverse line orientation in a line loop

I am getting around this by manually entering a - sign in front of

for k in range(nslices-1):
    ll = geom.add_line_loop([
        gmsh_hlines[n][k],
        gmsh_lines[n][k],
        '-' + gmsh_hlines[n][k+1],
        '-' + gmsh_lines[n-1][k]
        ])
    s = geom.add_plane_surface(ll)

Simplify process of creating planar shapes with holes

Hi @nschloe,

Great work! I have just started to work with Gmsh and thought a Python interface would be quite nice. Initially I wanted to develop it myself but then I found your package. I have some suggestions for a few small additions which would be much appreciated.

Currently, I want to mesh a rectangular domain with a polygonal-shaped hole (e.g. flow around a cylinder, airfoil, etc.). For this I need two Line Loops: one for the rectangle and one for the polygon. Then I define a Planar Surface using these two loops, et voila!

This can be done using your package quite straightforwardly by first creating points, then lines, then line loops and finally the planar surface. No problem. However, because you already provide the very sweet method add_polygon, I suggest you to also add something like add_polygon_loop. This way, the user does not have to create points and lines manually just for having a hole in a surface.

So you already have this:

def add_polygon(self, X, lcar):
        # Create points.
        p = [self.add_point(x, lcar) for x in X]
        # Create lines
        e = [self.add_line(p[k], p[k+1]) for k in range(len(p)-1)]
        e.append(self.add_line(p[-1], p[0]))
        ll = self.add_line_loop(e)
        s = self.add_plane_surface(ll)
        return s

from that you can create add_polygon_loop with no effort at all:

def add_polygon_loop(self, X, lcar):
        # Create points.
        p = [self.add_point(x, lcar) for x in X]
        # Create lines
        e = [self.add_line(p[k], p[k+1]) for k in range(len(p)-1)]
        e.append(self.add_line(p[-1], p[0]))
        ll = self.add_line_loop(e)
        return ll

Similar things can also be done for all other methods that provide planar shapes (e.g. for the circle, rectangle etc.).

I am going to fork your project and add this for the polygon case. I will send you a pull request when I'm done. Let me know what you think of this.

Disk Example: Changing Order of Elements in Mesh?

Hi again,

This is somewhat of an extension of my previous question, but lies at the crux of my problems. If there is a better medium through which I should ask questions, please let me know.

Question 3
For the following simple example of a disk (code that follows works)
"Generate Disk Geometry"

geom = pygmsh.built_in.Geometry()

#add points for circle
x0 = (0.,0., 0.)
radius = 1
lcar = 10
#approximate circle by polygon with #edges = 2*num_sections
num_sections = 40
poly = geom.add_circle(x0, radius, lcar, R=None, compound=False, num_sections=num_sections, holes=None, make_surface=True)

#generate mesh from geom
points, cells, point_data, cell_data, field_data = pygmsh.generate_mesh(geom)

import meshio
meshio.write('disk.vtu', points, cells, cell_data=cell_data)

How do I change the 'order' of the elements (i.e. the number of nodes on a given triangle, in this example).

Thanks, and Happy Thanksgiving
Charlie