Deforming a mesh

mesh
1

Hi everyone!

I am trying to deform the mesh from some sensitivity results. I have used this code:

mesh = domain
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
boundary_facets = dolfinx.mesh.exterior_facet_indices(mesh.topology)
boundary_vertices = dolfinx.mesh.compute_incident_entities(mesh, boundary_facets, mesh.topology.dim-1, 0)
vertex_to_geometry = dolfinx.cpp.mesh.entities_to_geometry(mesh._mesh, 0, boundary_vertices, False)


c_to_v = mesh.topology.connectivity(mesh.topology.dim, 0)
mesh.topology.create_connectivity(0, mesh.topology.dim)
v_to_c = mesh.topology.connectivity(0, mesh.topology.dim)

num_dofs = V.dofmap.index_map.size_local+ V.dofmap.index_map.num_ghosts
dof_to_geometry_map = np.full(num_dofs, -1, dtype=np.int32)
dofmap = V.dofmap
layout = dofmap.dof_layout

for i, (vertex, node) in enumerate(zip(boundary_vertices, vertex_to_geometry)):
    cell = v_to_c.links(vertex)[0]
    cell_dofs = dofmap.cell_dofs(cell)
    cvs = c_to_v.links(cell)
    local_index = np.flatnonzero(cvs == vertex)[0]
    dof = cell_dofs[layout.entity_dofs(0, local_index)]
    dof_to_geometry_map[dof] = node

# Create boundary perturbation vector
bs = V.dofmap.bs
perturbation_data = np.zeros((len(boundary_vertices), 3), dtype=np.float64)
geom = np.zeros(len(boundary_vertices), dtype=np.int32)
c = 0
for i, node in enumerate(dof_to_geometry_map):
    if node != -1:
        perturbation_data[c, :bs] = sensi_func1.x.array[bs*i:bs*(i+1)]
        geom[c] = node
        c+=1


mesh.geometry.x[geom]+= perturbation_data
with dolfinx.io.XDMFFile(mesh.comm, "Plot_mesh/mesh_def_70.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

and I get

mesh_def
mesh_def1158×784 302 KB

Is there any way to solve it?
This is the code of how I generate the original mesh

gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0

gmsh.initialize()
     
    #Points
n_point = 70
str_point = str(n_point)


##################################################################################

    # Coordinates
    
##################################################################################


cx = [0,0.225,0.252,0.292,0.3737,0.9450,0.9450,0.3737,0.292,0.252,0.225,0]
cy = [0.025, 0.025,0.02,0.02,0.025,0.025,-0.025, -0.025,-0.02,-0.02,-0.025,-0.025]


bump = 0.05



a = 2.0*(bump-1)/(n_point*1.3)+1 

n_bump = int(math.log(bump,a))
c = 1e-1;
gmsh.model.geo.add_point(0, 0.025, 0, c,1)
gmsh.model.geo.add_point(0.225, 0.025, 0, c,2)
gmsh.model.geo.add_point(0.252, 0.02, 0, c,3)
gmsh.model.geo.add_point(0.292, 0.02, 0, c, 4)
gmsh.model.geo.add_point(0.3737, 0.025, 0, c, 5)
gmsh.model.geo.add_point(0.9450, 0.025, 0, c, 6)
gmsh.model.geo.add_point(0.9450, -0.025, 0, c, 7)
gmsh.model.geo.add_point(0.3737, -0.025, 0, c,8)
gmsh.model.geo.add_point(0.292, -0.02, 0, c, 9)
gmsh.model.geo.add_point(0.252, -0.02, 0, c, 10)
gmsh.model.geo.add_point(0.225, -0.025, 0, c, 11)
gmsh.model.geo.add_point(0, -0.025, 0, c, 12)

    #Lines
    
gmsh.model.geo.addLine(1, 2 ,1)
gmsh.model.geo.addLine(2, 3, 2)
gmsh.model.geo.addLine(3, 4,3)
gmsh.model.geo.addLine(4, 5,4)
gmsh.model.geo.addLine(5, 6,5)
gmsh.model.geo.addLine(6, 7,6)
gmsh.model.geo.addLine(7, 8,7)
gmsh.model.geo.addLine(8, 9,8)
gmsh.model.geo.addLine(9, 10,9)
gmsh.model.geo.addLine(10, 11,10)
gmsh.model.geo.addLine(11, 12,11)
gmsh.model.geo.addLine(12, 1,12)


gmsh.model.geo.add_curve_loop([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12],1)
gmsh.model.geo.add_plane_surface([1],1)
gmsh.model.geo.synchronize()
volumes = gmsh.model.getEntities(dim=gdim)


gmsh.model.addPhysicalGroup(1, [12], 1)
gmsh.model.addPhysicalGroup(1, [6], name = 'inlet')
gmsh.model.addPhysicalGroup(1, [1,2,3,4,5,7,8,9,10,11], name= 'walls')
gmsh.model.addPhysicalGroup(2, [1], 1)
gmsh.model.occ.synchronize()


gmsh.model.geo.mesh.setTransfiniteCurve(6, n_bump, "Bump", bump)
gmsh.model.geo.mesh.setTransfiniteCurve(12, n_bump, "Bump", bump)
gmsh.model.geo.mesh.setTransfiniteCurve(1, int(5.2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(2, int(0.8*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(3, int(1*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(4, int(2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(5, int(10*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(7, int(10*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(8, int(2*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(9, int(1*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(10, int(0.8*n_point))
gmsh.model.geo.mesh.setTransfiniteCurve(11, int(5.2*n_point))

gmsh.model.geo.mesh.setTransfiniteSurface(1, "Left", [1,6,7,12])


gmsh.model.geo.mesh.setRecombine(2, 1)
gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(2)

domain, _, ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=2)
ft.name = "Facet markers"
gmsh.finalize()

Thank you so much

2

If you are only perturbing boundary nodes, you will need to apply smoothing to the interior nodes to make sure that you Get a valid deformation, see

1 Like
3

Thank you so much @dokken! I have another question. Is it possible to save the new mesh in a format that can read gmsh?

4

(User trying to help.) For what I know, there is an io module which uses meshio in the background to export (?), but you may already know this so I would check what meshio provides

5

meshio is what you are looking for. You can read an xdmf into gmsh in the following manner.

import gmsh
import meshio

mesh = mesio.read('PATH/FILENAME.xdmf')
mesh.write('PATH/FILENAME_TO_WRITE.gmsh', file_format='gmsh')

gmsh. initialize()
gmsh.open('PATH/FILENAME.gmsh')
#.....
# gmsh code
#....
gmsh.finalize()

You may have also seen references to gmshio somewhere, and this is a wrapper for some of this functionality.

If you are looking to avoid the file I/O, have a look at this tutorial as well.

6

I get this error:

File h5py/_objects.pyx:54, in h5py._objects.with_phil.wrapper()

File h5py/_objects.pyx:55, in h5py._objects.with_phil.wrapper()

File h5py/h5p.pyx:1296, in h5py.h5p.PropFAID.set_libver_bounds()

ValueError: High bound is not valid (high bound is not valid)
7

Thank you @jkrokowski
It’s worked for me. I had problems with HDF5, but I solved too.

1 Like