FEniCS / dolfinx

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

[BUG]: Can't create mesh #2680

Closed ruf10 closed 1 year ago

ruf10 commented 1 year ago

How to reproduce the bug

import numpy as np import warnings from dolfinx.io import XDMFFile, gmshio from mpi4py import MPI import gmsh import os comm = MPI.COMM_WORLD rank = comm.Get_rank()

Need to figure out how to get hmin and hmax

print("Gmsh version:", gmsh.version)

Create "taylor_couette" mesh

gmsh.option.setNumber("General.Terminal", 0) warnings.filterwarnings("ignore") gmsh.initialize() gmsh.model.add("taylor_couette")

height = 6 # height of the cylinder R = 1 r = 0.5 cyl_in = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, height, r) cyl_out = gmsh.model.occ.addCylinder(0, 0, 0, 0,0, height, R) fluid = gmsh.model.occ.cut([(3, cyl_out)], [(3, cyl_in)],removeTool=True) gmsh.model.occ.synchronize() gmsh.model.mesh.generate(3) gmsh.model.mesh.refine() gmsh.model.mesh.refine()

gmsh.model.mesh.generate(3)

volumes = gmsh.model.getEntities(dim=3) assert(volumes == fluid[0]) fluid_marker = 11 gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker) gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid volume")

impose periodic boundary condtion on z direction

surfaces = gmsh.model.getEntities(dim=2) # Get 2D surfaces print("surfaces") print(surfaces)

Find the surfaces corresponding to the top and bottom regions

top_surface_tag = -1 bottom_surface_tag = -1

Loop through the surfaces to find the top and bottom surfaces

for dim, tag in surfaces: vertices = gmsh.model.mesh.getNodes(dim, tag) for node_index in vertices[0]: node_coords = gmsh.model.mesh.getNode(node_index) z_values = node_coords[0][2] #z values if np.allclose(z_values, height): # Check if all z vals are close to 2.0 (top surface) top_surface_tag = tag break elif np.allclose(z_values, 0.0): # Check if all z vals are close to 0.0 (bottom surface) bottom_surface_tag = tag break

Check if the top and bottom surfaces have been found

if top_surface_tag == -1 or bottom_surface_tag == -1: raise ValueError("Unable to find the top and bottom surfaces")

Define the periodic boundary condition

Define the translation vector

translation_vector = [0., 0., height]

Create the affine transformation matrix

affine_matrix = [1., 0.0, 0.0, translation_vector[0], 0., 1., 0., translation_vector[1], 0., 0., 1., translation_vector[2], 0, 0., 0., 1.] gmsh.model.mesh.setPeriodic(2, [top_surface_tag], [bottom_surface_tag], affine_matrix)

Set the physical names for the top and bottom surfaces

top_surface_marker = 13 bottom_surface_marker = 14

here surafces[0][0]=2, the dim

gmsh.model.addPhysicalGroup(2, [top_surface_tag], top_surface_marker) gmsh.model.addPhysicalGroup(2, [bottom_surface_tag], bottom_surface_marker) gmsh.model.setPhysicalName(2, top_surface_marker, "Top Surface") gmsh.model.setPhysicalName(2, bottom_surface_marker, "Bottom Surface") gmsh.write("mesh3D.msh")

plot mesh

import pyvista as pv p = pv.Plotter() meshio_mesh = pv.read("mesh3D.msh") p.add_mesh(meshio_mesh, color=True)

Save the plot to a file

p.show(screenshot="tc_new.png")

msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_SELF, 0) msh.name = "Sphere" cell_markers.name = f"{msh.name}_cells" facet_markers.name = f"{msh.name}_facets"

Minimal Example (Python)

import numpy as np
import warnings
from dolfinx.io import XDMFFile, gmshio
from mpi4py import MPI
import gmsh
import os
comm = MPI.COMM_WORLD
rank = comm.Get_rank()

#Need to figure out how to get hmin and hmax
print("Gmsh version:", gmsh.__version__)

#Create "taylor_couette" mesh
gmsh.option.setNumber("General.Terminal", 0)
warnings.filterwarnings("ignore")
gmsh.initialize()
gmsh.model.add("taylor_couette")

height = 6 # height of the cylinder
R = 1
r = 0.5
cyl_in = gmsh.model.occ.addCylinder(0, 0, 0, 0, 0, height, r)
cyl_out = gmsh.model.occ.addCylinder(0, 0, 0, 0,0, height, R)
fluid = gmsh.model.occ.cut([(3, cyl_out)], [(3, cyl_in)],removeTool=True)
gmsh.model.occ.synchronize()
gmsh.model.mesh.generate(3)
gmsh.model.mesh.refine()
gmsh.model.mesh.refine()

#gmsh.model.mesh.generate(3)
volumes = gmsh.model.getEntities(dim=3)
assert(volumes == fluid[0])
fluid_marker = 11
gmsh.model.addPhysicalGroup(volumes[0][0], [volumes[0][1]], fluid_marker)
gmsh.model.setPhysicalName(volumes[0][0], fluid_marker, "Fluid volume")
#impose periodic boundary condtion on z direction 
surfaces = gmsh.model.getEntities(dim=2) # Get 2D surfaces
print("surfaces")
print(surfaces)
# Find the surfaces corresponding to the top and bottom regions
top_surface_tag = -1
bottom_surface_tag = -1
# Loop through the surfaces to find the top and bottom surfaces
for dim, tag in surfaces:
    vertices = gmsh.model.mesh.getNodes(dim, tag)
    for node_index in vertices[0]:
        node_coords = gmsh.model.mesh.getNode(node_index)
        z_values = node_coords[0][2] #z values
        if np.allclose(z_values, height):  # Check if all z vals are close to 2.0 (top surface)
            top_surface_tag = tag
            break
        elif np.allclose(z_values, 0.0):  # Check if all z vals are close to 0.0 (bottom surface)
            bottom_surface_tag = tag
            break

# Check if the top and bottom surfaces have been found
if top_surface_tag == -1 or bottom_surface_tag == -1:
    raise ValueError("Unable to find the top and bottom surfaces")

# Define the periodic boundary condition
# Define the translation vector
translation_vector = [0., 0., height]
# Create the affine transformation matrix
affine_matrix = [1., 0.0, 0.0, translation_vector[0],
                 0., 1., 0., translation_vector[1],
                 0., 0., 1., translation_vector[2],
                 0, 0., 0., 1.]
gmsh.model.mesh.setPeriodic(2, [top_surface_tag], [bottom_surface_tag], affine_matrix)
# Set the physical names for the top and bottom surfaces
top_surface_marker = 13
bottom_surface_marker = 14 
#here surafces[0][0]=2, the dim
gmsh.model.addPhysicalGroup(2, [top_surface_tag], top_surface_marker)
gmsh.model.addPhysicalGroup(2, [bottom_surface_tag], bottom_surface_marker)
gmsh.model.setPhysicalName(2, top_surface_marker, "Top Surface")
gmsh.model.setPhysicalName(2, bottom_surface_marker, "Bottom Surface")
gmsh.write("mesh3D.msh")
#plot mesh
import pyvista as pv
p = pv.Plotter()
meshio_mesh = pv.read("mesh3D.msh")
p.add_mesh(meshio_mesh, color=True)
# # Save the plot to a file
p.show(screenshot="tc_new.png")

msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_SELF, 0)
msh.name = "Sphere"
cell_markers.name = f"{msh.name}_cells"
facet_markers.name = f"{msh.name}_facets"

Output (Python)

msh, cell_markers, facet_markers = gmshio.model_to_mesh(gmsh.model, MPI.COMM_SELF, 0)
  File "/ix/wlayton/rui/miniconda/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/io/gmshio.py", line 250, in model_to_mesh
    mesh = create_mesh(comm, cells, x[:, :gdim], ufl_domain, partitioner)
  File "/ix/wlayton/rui/miniconda/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/mesh.py", line 177, in create_mesh
    cmap = _cpp.fem.CoordinateElement(_uflcell_to_dolfinxcell[cell_shape], cell_degree, variant)
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfinx.cpp.fem.CoordinateElement(celltype: dolfinx.cpp.mesh.CellType, degree: int)
    2. dolfinx.cpp.fem.CoordinateElement(celltype: dolfinx.cpp.mesh.CellType, degree: int, variant: basix::element::lagrange_variant)

Version

main branch

DOLFINx git commit

No response

Installation

No response

Additional information

No response

ruf10 commented 1 year ago

also, I want to impose periodic boundary condition on the z direction, customize the size of my mesh, I can only do refine, but not specify which part I want to refine. Can anyone help me with that? Thank you so much!

I am trying to create a concentric 3 D cylinders

jorgensd commented 1 year ago

also, I want to impose periodic boundary condition on the z direction, customize the size of my mesh, I can only do refine, but not specify which part I want to refine. Can anyone help me with that? Thank you so much!

I am trying to create a concentric 3 D cylinders

Please note that parts of your post above is not properly formatted. Please use 3x` encapsulation on that code.

For the error message you are printing above, it looks similar to: https://fenicsproject.discourse.group/t/error-converting-mesh-from-gmsh-to-dolfinx/11554/9?u=dokken

See the temporary fix there.

For the last part of your question, please post this on discourse, as this is not an issue with your code, but a question/request for help.