search

FEniCS / dolfinx

Next generation FEniCS problem solving environment
https://fenicsproject.org
GNU Lesser General Public License v3.0
739 stars 178 forks source link

[BUG]: extract_geometry expects no point tags in the gmsh geometry #2460

Closed niravshah241 closed 1 year ago

niravshah241 commented 1 year ago

extract_geometry in gmshio assumes the indices to start with 1. However, this assumption is hidden from the user. Also, it is not always possible to create meshes without any point or edge tags.

Adding a good error message such as (assuming openCASCADE kernel):

assert min(indices) == 1, "Consider removing all point tags in gmsh using gmsh.model.occ.remove([0,point_tag_number])"

after 

indices, points, _ = model.mesh.getNodes()

in extract_geometry in gmshio will be helpful to users for understanding and for resolving the error.

Also similar changes for edge tags in extract_topology_and_markers should be considered.

How to reproduce the bug

Reproduce the error by commenting/removing below lines in the minimal example.

for i in range(1,48):
    gmsh.model.occ.remove([(0,i)])

Without removing point tags, the indices start at 48 (i.e. the indices start after the tags in gmsh geometry and do not start at 1).

Minimal Example (Python)

import gmsh 

gmsh.initialize('',False) 

gmsh.model.occ.addPoint(0.000000,1.7700E-2,0.,0.1,1)
gmsh.model.occ.addPoint(0.002300,0.030900,0.,0.1,2)
gmsh.model.occ.addPoint(0.005000,0.037200,0.,0.1,3)
gmsh.model.occ.addPoint(0.007600,0.041500,0.,0.1,4)
gmsh.model.occ.addPoint(0.014300,0.049900,0.,0.1,5)
gmsh.model.occ.addPoint(0.024900,0.058200,0.,0.1,6)
gmsh.model.occ.addPoint(0.049500,0.073000,0.,0.1,7)
gmsh.model.occ.addPoint(0.074000,0.081400,0.,0.1,8)
gmsh.model.occ.addPoint(0.099000,0.086600,0.,0.1,9)
gmsh.model.occ.addPoint(0.153000,0.090700,0.,0.1,10)
gmsh.model.occ.addPoint(0.196100,0.090500,0.,0.1,11)
gmsh.model.occ.addPoint(0.250400,0.088700,0.,0.1,12)
gmsh.model.occ.addPoint(0.309400,0.085800,0.,0.1,13)
gmsh.model.occ.addPoint(0.352000,0.083300,0.,0.1,14)
gmsh.model.occ.addPoint(0.391900,0.080400,0.,0.1,15)
gmsh.model.occ.addPoint(0.447700,0.075600,0.,0.1,16)
gmsh.model.occ.addPoint(0.503400,0.069600,0.,0.1,17)
gmsh.model.occ.addPoint(0.559300,6.2600E-2,0.,0.1,18)
gmsh.model.occ.addPoint(0.596500,0.057500,0.,0.1,19)
gmsh.model.occ.addPoint(0.648800,0.049800,0.,0.1,20)
gmsh.model.occ.addPoint(0.835100,0.022400,0.,0.1,21)
gmsh.model.occ.addPoint(0.910900,0.013200,0.,0.1,22)
gmsh.model.occ.addPoint(1.000000,0.000300,0.,0.1,23)
gmsh.model.occ.addPoint(0.000000,0.017700,0.,0.1,25)
gmsh.model.occ.addPoint(0.002200,0.003800,0.,0.1,26)
gmsh.model.occ.addPoint(0.004900,-0.001800,0.,0.1,27)
gmsh.model.occ.addPoint(0.007200,-0.005300,0.,0.1,28)
gmsh.model.occ.addPoint(0.011900,-0.010600,0.,0.1,29)
gmsh.model.occ.addPoint(0.024300,-0.020400,0.,0.1,30)
gmsh.model.occ.addPoint(0.048600,-0.034200,0.,0.1,31)
gmsh.model.occ.addPoint(0.071600,-0.045700,0.,0.1,32)
gmsh.model.occ.addPoint(0.097900,-0.051600,0.,0.1,33)
gmsh.model.occ.addPoint(0.148800,-0.060700,0.,0.1,34)
gmsh.model.occ.addPoint(0.195300,-0.063200,0.,0.1,35)
gmsh.model.occ.addPoint(0.250100,-0.063200,0.,0.1,36)
gmsh.model.occ.addPoint(0.294500,-0.062600,0.,0.1,37)
gmsh.model.occ.addPoint(0.357900,-0.061000,0.,0.1,38)
gmsh.model.occ.addPoint(0.396500,-0.059500,0.,0.1,39)
gmsh.model.occ.addPoint(0.454300,-0.056300,0.,0.1,40)
gmsh.model.occ.addPoint(0.505000,-0.052700,0.,0.1,41)
gmsh.model.occ.addPoint(0.555600,-0.048200,0.,0.1,42)
gmsh.model.occ.addPoint(0.606300,-0.042700,0.,0.1,43)
gmsh.model.occ.addPoint(0.648500,-0.037500,0.,0.1,44)
gmsh.model.occ.addPoint(0.831700,-0.014900,0.,0.1,45)
gmsh.model.occ.addPoint(0.941000,-0.005300,0.,0.1,46)
gmsh.model.occ.addPoint(1.000000,-0.000300,0.,0.1,47)

gdim = 2 

spline_curves = gmsh.model.occ.addSpline([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 1])
spline_curve_loop = gmsh.model.occ.addCurveLoop([spline_curves])
spline_surface = gmsh.model.occ.addPlaneSurface([spline_curve_loop])

for i in range(1,48):
    gmsh.model.occ.remove([(0,i)])

gmsh.model.occ.synchronize()

gmsh.option.setNumber("Mesh.Algorithm", 8) # 8=Frontal-Delaunay for Quads (See section 7.4,  https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2) # 2=simple full-quad (See section 7.4,  https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options)
gmsh.option.setNumber("Mesh.RecombineAll", 1) # Apply recombination algorithm to all surfaces, ignoring per-surface spec
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1) # Mesh subdivision algorithm (0: none, 1: all quadrangles, 2: all hexahedra, 3: barycentric)
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.1) # Minimum characteristic element size
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.3) # Maximum characteristic element size
gmsh.model.mesh.generate(gdim) # Mesh generation
gmsh.model.mesh.setOrder(1) # Mesh order
gmsh.model.mesh.optimize("Netgen") # Mesh optimisation or improving quality of mesh

# Extract edges and surfaces to add physical groups
surfaces = gmsh.model.getEntities(dim=gdim)
edges = gmsh.model.getBoundary(surfaces) # Gives 'list' of boundaries in the form [(gdim-1),marker] with length = number of boundaries
# edges = gmsh.model.getEntities(dim=gdim-1) # Alternate to getting boundaries
print(surfaces,edges)

for i in range(1,len(surfaces)+1):
    gmsh.model.addPhysicalGroup(gdim,[surfaces[i-1][1]],surfaces[i-1][1])
for i in range(1,len(edges)+1):
    gmsh.model.addPhysicalGroup(gdim-1,[edges[i-1][1]],edges[i-1][1])

gmsh.model.occ.remove(gmsh.model.getEntities(dim=gdim-1))
gmsh.model.occ.remove(gmsh.model.getEntities(dim=gdim))

gmsh.write("domain_geometry.msh")
gmsh.fltk.run()

import os
import numpy as np
import matplotlib.pyplot as plt

from mpi4py import MPI
from petsc4py import PETSc

from dolfinx.mesh import *
from dolfinx.fem import *
from dolfinx.io import *
from ufl import *

import ufl
import dolfinx

gdim = 2
#mesh, cell_tags, facet_tags = gmshio.read_from_msh("domain_geometry.msh", MPI.COMM_WORLD, 0, gdim=gdim)

#with XDMFFile(MPI.COMM_WORLD, "mesh_data/mesh.xdmf", "w") as mesh_file_xdmf:
    #mesh_file_xdmf.write_mesh(domain)
    #mesh_file_xdmf.write_meshtags(facet_tags)

# Import mesh in dolfinx
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD

domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)

assert domain.topology.dim == gdim, "Geometric dimensions of mesh and domain are not same"

with XDMFFile(MPI.COMM_WORLD, "mesh_data/mesh.xdmf", "w") as mesh_file_xdmf:
    mesh_file_xdmf.write_mesh(domain)
    mesh_file_xdmf.write_meshtags(facet_markers)

gmsh.finalize()

Output (Python)

[0]PETSC ERROR: ------------------------------------------------------------------------
[0]PETSC ERROR: Caught signal number 11 SEGV: Segmentation Violation, probably memory access out of range
[0]PETSC ERROR: Try option -start_in_debugger or -on_error_attach_debugger
[0]PETSC ERROR: or see https://petsc.org/release/faq/#valgrind and https://petsc.org/release/faq/
[0]PETSC ERROR: ---------------------  Stack Frames ------------------------------------
[0]PETSC ERROR: No error traceback is available, the problem could be in the main program. 
application called MPI_Abort(MPI_COMM_WORLD, 59) - process 0
[unset]: write_line error; fd=-1 buf=:cmd=abort exitcode=59
:
system msg for write_line failure : Bad file descriptor

or

Traceback (most recent call last):
  File "test.py", line 121, in <module>
    domain, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
  File "***********/python3.8/site-packages/dolfinx/io/gmshio.py", line 255, in model_to_mesh
    local_entities, local_values = _cpp.io.distribute_entity_data(mesh, mesh.topology.dim, cells, cell_values)
MemoryError: std::bad_alloc

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

jorgensd commented 1 year ago

@niravshah241 if you had used the gmsh python interface to save the file it would automatically prune these, as Mesh.SaveAll by default is 0. See https://gmsh.info/doc/texinfo/gmsh.html#Geometry-options 

Mesh.SaveAll

    Save all elements, even if they don’t belong to physical groups (for some mesh formats, this removes physical groups altogether)
    Default value: 0
    Saved in: -

The following modified version of your code, using either read_from_msh or model_to_mesh works:

from dolfinx.io import gmshio, XDMFFile

from mpi4py import MPI
import gmsh

gmsh.initialize('', False)

gmsh.model.occ.addPoint(0.000000, 1.7700E-2, 0., 0.1, 1)
gmsh.model.occ.addPoint(0.002300, 0.030900, 0., 0.1, 2)
gmsh.model.occ.addPoint(0.005000, 0.037200, 0., 0.1, 3)
gmsh.model.occ.addPoint(0.007600, 0.041500, 0., 0.1, 4)
gmsh.model.occ.addPoint(0.014300, 0.049900, 0., 0.1, 5)
gmsh.model.occ.addPoint(0.024900, 0.058200, 0., 0.1, 6)
gmsh.model.occ.addPoint(0.049500, 0.073000, 0., 0.1, 7)
gmsh.model.occ.addPoint(0.074000, 0.081400, 0., 0.1, 8)
gmsh.model.occ.addPoint(0.099000, 0.086600, 0., 0.1, 9)
gmsh.model.occ.addPoint(0.153000, 0.090700, 0., 0.1, 10)
gmsh.model.occ.addPoint(0.196100, 0.090500, 0., 0.1, 11)
gmsh.model.occ.addPoint(0.250400, 0.088700, 0., 0.1, 12)
gmsh.model.occ.addPoint(0.309400, 0.085800, 0., 0.1, 13)
gmsh.model.occ.addPoint(0.352000, 0.083300, 0., 0.1, 14)
gmsh.model.occ.addPoint(0.391900, 0.080400, 0., 0.1, 15)
gmsh.model.occ.addPoint(0.447700, 0.075600, 0., 0.1, 16)
gmsh.model.occ.addPoint(0.503400, 0.069600, 0., 0.1, 17)
gmsh.model.occ.addPoint(0.559300, 6.2600E-2, 0., 0.1, 18)
gmsh.model.occ.addPoint(0.596500, 0.057500, 0., 0.1, 19)
gmsh.model.occ.addPoint(0.648800, 0.049800, 0., 0.1, 20)
gmsh.model.occ.addPoint(0.835100, 0.022400, 0., 0.1, 21)
gmsh.model.occ.addPoint(0.910900, 0.013200, 0., 0.1, 22)
gmsh.model.occ.addPoint(1.000000, 0.000300, 0., 0.1, 23)
gmsh.model.occ.addPoint(0.000000, 0.017700, 0., 0.1, 25)
gmsh.model.occ.addPoint(0.002200, 0.003800, 0., 0.1, 26)
gmsh.model.occ.addPoint(0.004900, -0.001800, 0., 0.1, 27)
gmsh.model.occ.addPoint(0.007200, -0.005300, 0., 0.1, 28)
gmsh.model.occ.addPoint(0.011900, -0.010600, 0., 0.1, 29)
gmsh.model.occ.addPoint(0.024300, -0.020400, 0., 0.1, 30)
gmsh.model.occ.addPoint(0.048600, -0.034200, 0., 0.1, 31)
gmsh.model.occ.addPoint(0.071600, -0.045700, 0., 0.1, 32)
gmsh.model.occ.addPoint(0.097900, -0.051600, 0., 0.1, 33)
gmsh.model.occ.addPoint(0.148800, -0.060700, 0., 0.1, 34)
gmsh.model.occ.addPoint(0.195300, -0.063200, 0., 0.1, 35)
gmsh.model.occ.addPoint(0.250100, -0.063200, 0., 0.1, 36)
gmsh.model.occ.addPoint(0.294500, -0.062600, 0., 0.1, 37)
gmsh.model.occ.addPoint(0.357900, -0.061000, 0., 0.1, 38)
gmsh.model.occ.addPoint(0.396500, -0.059500, 0., 0.1, 39)
gmsh.model.occ.addPoint(0.454300, -0.056300, 0., 0.1, 40)
gmsh.model.occ.addPoint(0.505000, -0.052700, 0., 0.1, 41)
gmsh.model.occ.addPoint(0.555600, -0.048200, 0., 0.1, 42)
gmsh.model.occ.addPoint(0.606300, -0.042700, 0., 0.1, 43)
gmsh.model.occ.addPoint(0.648500, -0.037500, 0., 0.1, 44)
gmsh.model.occ.addPoint(0.831700, -0.014900, 0., 0.1, 45)
gmsh.model.occ.addPoint(0.941000, -0.005300, 0., 0.1, 46)
gmsh.model.occ.addPoint(1.000000, -0.000300, 0., 0.1, 47)

gdim = 2

spline_curves = gmsh.model.occ.addSpline([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19,
                                         20, 21, 22, 23, 47, 46, 45, 44, 43, 42, 41, 40, 39, 38, 37, 36, 35, 34, 33, 32, 31, 30, 29, 28, 27, 26, 1])
spline_curve_loop = gmsh.model.occ.addCurveLoop([spline_curves])
spline_surface = gmsh.model.occ.addPlaneSurface([spline_curve_loop])

gmsh.model.occ.synchronize()

# 8=Frontal-Delaunay for Quads (See section 7.4,  https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options)
gmsh.option.setNumber("Mesh.Algorithm", 8)
# 2=simple full-quad (See section 7.4,  https://gmsh.info/doc/texinfo/gmsh.html#Mesh-options)
gmsh.option.setNumber("Mesh.RecombinationAlgorithm", 2)
# Apply recombination algorithm to all surfaces, ignoring per-surface spec
gmsh.option.setNumber("Mesh.RecombineAll", 1)
# Mesh subdivision algorithm (0: none, 1: all quadrangles, 2: all hexahedra, 3: barycentric)
gmsh.option.setNumber("Mesh.SubdivisionAlgorithm", 1)
# Minimum characteristic element size
gmsh.option.setNumber("Mesh.MeshSizeMin", 0.1)
# Maximum characteristic element size
gmsh.option.setNumber("Mesh.MeshSizeMax", 0.3)
gmsh.model.mesh.generate(gdim)  # Mesh generation
gmsh.model.mesh.setOrder(1)  # Mesh order
# Mesh optimisation or improving quality of mesh
gmsh.model.mesh.optimize("Netgen")

# Extract edges and surfaces to add physical groups
surfaces = gmsh.model.getEntities(dim=gdim)
# Gives 'list' of boundaries in the form [(gdim-1),marker] with length = number of boundaries
edges = gmsh.model.getBoundary(surfaces)
# edges = gmsh.model.getEntities(dim=gdim-1) # Alternate to getting boundaries

for i in range(1, len(surfaces)+1):
    gmsh.model.addPhysicalGroup(gdim, [surfaces[i-1][1]], surfaces[i-1][1])
for i in range(1, len(edges)+1):
    gmsh.model.addPhysicalGroup(gdim-1, [edges[i-1][1]], edges[i-1][1])

gmsh.model.mesh.generate(2)
gmsh.write("domain_geometry.msh")

# Import mesh in dolfinx
gmsh_model_rank = 0
mesh_comm = MPI.COMM_WORLD

domain, cell_markers, facet_markers = gmshio.model_to_mesh(
    gmsh.model, mesh_comm, gmsh_model_rank, gdim=gdim)
gmsh.finalize()

# domain, cell_markers, facet_markers = gmshio.read_from_msh(
#     "domain_geometry.msh", mesh_comm, gdim=gdim)

assert domain.topology.dim == gdim, "Geometric dimensions of mesh and domain are not same"

with XDMFFile(domain.comm, "mesh_data/mesh.xdmf", "w") as mesh_file_xdmf:
    mesh_file_xdmf.write_mesh(domain)
    mesh_file_xdmf.write_meshtags(facet_markers)

I do not see any reason of adding special handling to the interface (making it harder to maintain).

jorgensd commented 1 year ago

@niravshah241 I am going to close this issue, as special handling should not be required if one uses the gmsh interface with only saving physical quantities.

niravshah241 commented 1 year ago

@jorgensd thanks