Mpirun performs mesh generation

[ZhangQi]

Hello, Dear Professor. This is a function that I wrote about generating moving meshes base on Can’t create 1D mesh in parallel, where Ulh is the numerical result of calculations on the original mesh, mesh modification with consistent topology, and this processing is performed serially when running, as well as when running on the original mesh, performance is good. But when I use mpirun -np n python3 my_file.py, If n is greater than 1, it will result in a new grid derangement. What I understand is that the transfer takes place only in the main process, and the data from the other processes is not added to the main process, what should I do? Looking forward to your guidance, thank you.

def get_displacement_at_geometry_nodes(msh, displacement_func):
   

    geom_nodes = msh.geometry.x.copy()  # (N_geom, 3)
    print(f"geom_nodes-Number of nodes: {len(geom_nodes)}")
    
    V_linear = fem.functionspace(msh, ("Lagrange", 1, (3,)))

    disp_linear = fem.Function(V_linear)
    disp_linear.interpolate(displacement_func)
    displacement_at_nodes = disp_linear.x.array.reshape(-1, 3)
    print(f"displacement_at_nodes-Number of nodes: {len(displacement_at_nodes)}")

    if MPI.COMM_WORLD.rank == 0:
        deformed_coords = msh.geometry.x.copy()
        deformed_coords[:] += displacement_at_nodes  # 

        #
        num_nodes = len(deformed_coords)
        print(f"deformed_coords-Number of nodes: {num_nodes}")

        # 
        mesh_obj = displacement_func.function_space.mesh

    
        num_nodes_global = mesh_obj.topology.index_map(mesh_obj.topology.dim).size_global

        cells = msh.topology.connectivity(3, 0)  #
        cell_vertices = cells.array.reshape((-1, 8))  # 
        print(f"deformed_coords---{num_nodes_global}")
    else:
        # 
        cell_vertices = np.empty((0, 8), dtype=np.int64)
        deformed_coords = np.empty((0, 3), dtype=np.float64)

    #
    from basix import CellType as BasixCellType

    #
    # 
    tet_element = ufl.Mesh(element("Lagrange", BasixCellType.hexahedron, 1, shape=(3,)) )

    #
    deformed_mesh = mesh.create_mesh(
        MPI.COMM_WORLD,
        cell_vertices,
        deformed_coords, #dian
        tet_element
    )
    cell_index_map = deformed_mesh.topology.index_map(deformed_mesh.topology.dim)
    print(
        f"Num cells local: {cell_index_map.size_local}\n Num cells global: {cell_index_map.size_global}"
    )

    #
    num_nodes = len(deformed_mesh.geometry.x)
    print(f"number of the nodes: {num_nodes}")
    #
    return disp_linear, deformed_mesh

disp_linear, deformed_mesh = get_displacement_at_geometry_nodes(msh, Ulh)

[ZhangQi]

Can anyone help me? Thanks a million!

[dokken]

If you want to keep the connectivity, and only move the nodes, please consider:

I need to create stuctured uniform, structuted non-uniform, and unstructured mesh in 2D

Dear @Monisaddam, Please note that you have only provided images, leaving it to me to interpret what the spacing between each cell is. Here I have extracted the relevant code from the post I previously referred to to illustrate how to do this: """ Move mesh based on geometrical function Author: Jørgen S. Dokken and Henrik N.T. Finsberg SPDX-License-Identifier: MIT """ from mpi4py import MPI import dolfinx import numpy as np def vertex_to_dof_map_vectorized(V): """Create a map from the…

which explains how to move the current mesh nodes, rather than recreating the whole mesh.

[ZhangQi]

Thank you for your reply. This part is also what I want to learn, but I need to describe my own problem more clearly. My problem is a linear elastic problem; when there is external pressure, the whole region will be deformed, rather than simply internal node position changes. I think this is part of the content moving grid, not just the internal interface moving. Looking forward to your guidance!

[dokken]

The example above moves the internal nodes based on the data in a dolfinx.fem.Function, which I assume is what you will have as input to your own problem. It doesn’t have any assumptions that would keep the exterior boundary fixed, that is just due to the specifics of the prescribed functions expr and expr2

@dokken

def expr(x):
        move_1 = x[1] * (0.5 - x[1]) * (x[1] < 0.5)
        move_2 = (x[1] - 1) * (x[1] - 0.5) * (x[1] > 0.5)
        return (np.zeros(x.shape[1]), move_1 + move_2)

    u = dolfinx.fem.Function(V)
    u.interpolate(expr)
    perturb_mesh(u)

    # Create the dof -> vertex map and vertex->node map
    with dolfinx.io.XDMFFile(mesh.comm, "moved_mesh.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)

    # Reset mesh
    u.x.array[:] *= -1
    perturb_mesh(u)

    def expr2(x, eps=1e-8):
        xi = np.floor(x[0] * N)
        yi = np.floor(x[1] * N)
        move_x = (
            0.2
            / N
            * np.cos(10 * np.pi * x[0])
            * (x[1] > eps)
            * (x[1] < 1 - eps)
            * (-1) ** yi
        )
        move_y = (
            0.2
            / N
            * np.cos(10 * np.pi * x[1])
            * (x[0] > eps)
            * (x[0] < 1 - eps)
            * (-1) ** xi
        )
        return (move_y, move_x)

If u is a solution from a elasticity probloem, you can simply pass that to perturb_mesh(u)

[ZhangQi]

I’m sorry. Yes, you are right. In addition, I also have a question: I would like to ask if FEniCSx 0.9.0 has a similar FEniCS support ALE function.

[dokken]

We do not have a specific ALE module in FEniCSx. As you can observe, moving the mesh given a function can be done as above.

If you want to make a smoothen function, for instance by a Laplace smoothing, this is covered in:

Changing the shape in the surface normal directions

The simplest mesh deformation can be described as follows: import dolfinx from mpi4py import MPI from typing import List import numpy as np def deform_mesh(V, bcs: List[dolfinx.fem.DirichletBCMetaClass]): mesh = V.mesh u = dolfinx.fem.Function(V) dolfinx.fem.petsc.set_bc(u.vector, bcs) deformation_array = u.x.array.reshape((-1, mesh.geometry.dim)) mesh.geometry.x[:, :mesh.geometry.dim] += deformation_array mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10) …

NOTE: In the post above it doesn’t use the perturb_mesh, but assumes a 1-1 relationship between the mesh nodes and dofs. I would recommend adapting the code such that it looks like

def deform_mesh(V, bcs: List[dolfinx.fem.DirichletBCMetaClass]):
    mesh = V.mesh
    uh = dolfinx.fem.Function(V)
    dolfinx.fem.petsc.set_bc(uh.vector, bcs)
    u = ufl.TrialFunction(V)
    v = ufl.TestFunction(V)
    a = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
    L = ufl.inner(dolfinx.fem.Constant(mesh, (0., 0.)), v)*ufl.dx
    problem = dolfinx.fem.petsc.LinearProblem(a, L, bcs, uh)
    problem.solve()
    perturb_mesh(uh) 

[ZhangQi]

Thank you for your guidance, which is very helpful to me; then I will do as you say. Thank you again for your patient guidance.