# Regarding the vector reordering problem caused by parallel grid partitioning?

**URL:** <https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682>\
**Category:** General\
**Created:** [January 7, 2025, 1:17pm UTC](https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682 "2025-01-07T13:17:04Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![wanghao](https://avatars.discourse-cdn.com/v4/letter/w/96bed5/32.png) [@wanghao](https://fenicsproject.discourse.group/u/wanghao)\
**Post date:** [January 7, 2025, 1:17pm UTC](https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682/1 "2025-01-07T13:17:04Z")

</div>

Dear all,  
Grid partitioning can lead to the problem of reordering the same vector. How can I obtain this index?（global to local map?）  
**My question is :**  
**How can I extract rows belonging to u from vector b while running in parallel?**  
The smallest example is as follows：

```auto
from mpi4py import MPI

try:
    from petsc4py import PETSc

    import dolfinx

    if not dolfinx.has_petsc:
        print("This demo requires DOLFINx to be compiled with PETSc enabled.")
        exit(0)
except ModuleNotFoundError:
    print("This demo requires petsc4py.")
    exit(0)

import numpy as np

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type, fem, la
from dolfinx.fem import (
    Constant,
    Function,
    dirichletbc,
    extract_function_spaces,
    form,
    functionspace,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import assemble_matrix_block, assemble_vector_block
from dolfinx.io import XDMFFile
from dolfinx.mesh import CellType, create_rectangle, locate_entities_boundary
from ufl import div, dx, grad, inner
from dolfinx.fem.petsc import (apply_lifting, assemble_matrix_block, assemble_vector_block, 
                                 create_vector_block, set_bc,create_matrix_block,create_matrix,create_vector)

comm = MPI.COMM_WORLD 
rank=MPI.COMM_WORLD.Get_rank() 
size = comm.Get_size() 

msh = create_rectangle(
    MPI.COMM_WORLD, [np.array([0, 0]), np.array([1, 1])], [32, 32], CellType.triangle
)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=default_real_type)
V, Q = functionspace(msh, P2), functionspace(msh, P1)

noslip = np.zeros(msh.geometry.dim, dtype=PETSc.ScalarType) 

# Define variational problem
(u, p) = ufl.TrialFunction(V), ufl.TrialFunction(Q)
(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)

x = ufl.SpatialCoordinate(msh)
f1=ufl.as_vector([x[0],x[1]])
f2=ufl.as_vector([x[1]**2,x[0]**2])
a = form([[inner(grad(u), grad(v)) * dx, inner(p, div(v)) * dx], [inner(div(u), q) * dx, None]])
L = form([inner(f1, v) * dx, inner(f1, v) * dx])  

bb= create_vector_block(L)      
assemble_vector_block(bb, L, a, bcs=[])  

ccb=bb.getArray()[:] 

print(f'rank{rank} size is {len(ccb)} ,values is',ccb)

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 9, 2025, 8:48am UTC](https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682/2 "2025-01-09T08:48:55Z")

</div>

As explained in [DOLFINx: The next generation FEniCS problem solving environment](https://zenodo.org/records/10447666) (Chapter 4.1).  
DOLFINx uses the concept of an IndexMap to describe the parallel layout of data.

In the case of a single function space, you can look at the index map of the dofmap, and in particular its local\_range (multiplied by the number of vector components, i.e. the multiplicative product of all entries in value shape (`np.prod(V.element.value_shape)`).

Example:

```python
from mpi4py import MPI

import ufl
import dolfinx.fem.petsc

from basix.ufl import element
import numpy as np
comm = MPI.COMM_WORLD 
rank=MPI.COMM_WORLD.Get_rank() 
size = comm.Get_size() 

msh = dolfinx.mesh.create_unit_square(comm, 10, 10)

P2 = element("Lagrange", msh.basix_cell(), 2, shape=(msh.geometry.dim,), dtype=dolfinx.default_real_type)
P1 = element("Lagrange", msh.basix_cell(), 1, dtype=dolfinx.default_real_type)
V = dolfinx.fem.functionspace(msh, P2)
Q = dolfinx.fem.functionspace(msh, P1)
local_V_range = np.asarray(V.dofmap.index_map.local_range)
local_Q_range = np.asarray(Q.dofmap.index_map.local_range)
print(MPI.COMM_WORLD.rank, local_V_range*V.dofmap.index_map_bs, local_Q_range*Q.dofmap.index_map_bs)

```

yields

```bash
0 [0 298] [0 40]
1 [298 584] [40 77]
2 [584 882] [77 121]

```

A blocked space is ordered such that all the dofs of the first block comes first, then all the dofs of the second block, i.e.  
(dofs\_V, dofs\_Q).  
However, when you can `getArray()` on a process, you only get the owned entries:

```python
num_owned_V_dofs = (local_V_range[1] - local_V_range[0])*V.dofmap.index_map_bs
num_owned_Q_dofs = (local_Q_range[1] - local_Q_range[0])*Q.dofmap.index_map_bs

(v, q) = ufl.TestFunction(V), ufl.TestFunction(Q)

one = dolfinx.fem.Constant(msh, dolfinx.default_scalar_type(1.0))
one_v = ufl.as_vector([one, one])
L = dolfinx.fem.form([ufl.inner(one_v, v) * ufl.dx, ufl.inner(one, q) * ufl.dx])  

bb = dolfinx.fem.petsc.create_vector_block(L)

ccb=bb.getArray()[:] 

assert len(ccb) == num_owned_V_dofs + num_owned_Q_dofs

```

i.e. `ccb=(V_dofs_local, Q_dofs_local)`.  
Thus, you can get the global index by:

```python
V_global_start = 0
Q_global_start = V.dofmap.index_map_bs * V.dofmap.index_map.size_global

global_indices = np.zeros( num_owned_V_dofs + num_owned_Q_dofs, dtype=np.int64)
global_indices[:num_owned_V_dofs] = np.arange(V_global_start + local_V_range[0]*V.dofmap.index_map_bs, V_global_start + local_V_range[1]*V.dofmap.index_map_bs, dtype=np.int64)
global_indices[num_owned_V_dofs:] = np.arange(Q_global_start + local_Q_range[0]*Q.dofmap.index_map_bs, Q_global_start + local_Q_range[1]*Q.dofmap.index_map_bs, dtype=np.int64)
print(global_indices)

```

---

<div class="post-metadata">

**Author:** ![Stein](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/stein/32/8751_2.png) [@Stein](https://fenicsproject.discourse.group/u/Stein)\
**Post date:** [January 9, 2025, 10:51am UTC](https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682/3 "2025-01-09T10:51:55Z")

</div>

Coincidentally I was just reading the cited section as this reply came up, so maybe I can ask a related follow-up question here.

How is it possible that `index_maps` are always contiguous? This seemed like a ‘special case’ to me, but from the paper I understand this is a basic property. E.g., the nodes in a mesh are partitioned and ownership is distributed. Naturally, an optimal partitioning does not yield contiguously sets of (a priorily defined) node numbers. Is the numbering of the nodes in the mesh then updated a-posteriori? It seems similar issues should arise all the time.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 9, 2025, 1:15pm UTC](https://fenicsproject.discourse.group/t/regarding-the-vector-reordering-problem-caused-by-parallel-grid-partitioning/16682/4 "2025-01-09T13:15:44Z")

</div>

> [@Stein](#):
>
> How is it possible that `index_maps` are always contiguous? This seemed like a ‘special case’ to me, but from the paper I understand this is a basic property. E.g., the nodes in a mesh are partitioned and ownership is distributed. Naturally, an optimal partitioning does not yield contiguously sets of (a priorily defined) node numbers. Is the numbering of the nodes in the mesh then updated a-posteriori? It seems similar issues should arise all the time.

For creating structures in DOLFINx, there are several things at play.

1. Mesh data is read in (cell-node connectivity and node coordinates), and these are distributed over M processes with a graph partitioner. This will try to find an optimal partitioning of the cells, minimizing communication.
2. Once the cells have been partitioned (and ghost information has been distributed), one has to assign ownership (and renumber) the mesh nodes, such that those owned by process 0 is numbered [0, M-1] those owned by process 2 is numbered [M, N-1], etc.  
Nodes that are on a process and not owned is assigned to the index map as ghosts.
3. Similar procedure for the vertices (a subset of mesh nodes if the geometry is higher order).
4. Any other structure that needs an index map is based of these initial partitions, with suitable amounts of ghosts included.

Thus, not every index-map is optimally partitioned, they are results of the initial mesh partitioning. In particular, you may create any index map you want in DOLFINx by using the relevant constructors:

> **[dolfinx.common — DOLFINx 0.9.0 documentation](https://docs.fenicsproject.org/dolfinx/v0.9.0/python/generated/dolfinx.common.html#dolfinx.common.IndexMap)**

This is for instance done for creating point cloud meshes in:

> <https://github.com/FEniCS/dolfinx/pull/3455#issue-2567705770>
>
> Allows for submeshes made of vertices.
> 
> Coupled with: https://github.com/FEniC…S/basix/pull/851
> 
> Would also allow for creating of point-cloud meshes with a simple Python or C++ code of 20 lines, as shown in https://github.com/FEniCS/dolfinx/issues/2593
> \`\`\`python
> \# Create a mesh consisting of points only
> \# Author: Jørgen S. Dokken
> \# SPDX-License-Identifier: MIT
> 
> from mpi4py import MPI
> import basix.ufl
> 
> import dolfinx
> import numpy as np
> import numpy.typing as npt
> import ufl
> 
> 
> 
> def create\_point\_mesh(comm: MPI.Intracomm, points:npt.NDArray\[np.float32\]|npt.NDArray\[np.float64\])-\>dolfinx.mesh.Mesh:
> """
> Create a mesh consisting of points only.
> 
> Note:
> No nodes are shared between processes.
>     
> Args:
> comm: MPI communicator to create the mesh on.
> points: Points local to the process in the mesh.
> """
> # Create mesh topology
> cells = np.arange(points.shape\[0\], dtype=np.int32).reshape(-1, 1)
> topology = dolfinx.cpp.mesh.Topology(MPI.COMM\_WORLD, dolfinx.mesh.CellType.point)
> num\_nodes\_local = cells.shape\[0\]
> imap = dolfinx.common.IndexMap(MPI.COMM\_WORLD, num\_nodes\_local)
> local\_range = imap.local\_range\[0\]
> igi = np.arange(num\_nodes\_local, dtype=np.int64)+local\_range
> topology.set\_index\_map(0, imap)
> topology.set\_connectivity(dolfinx.graph.adjacencylist(cells), 0, 0)
> 
> # Create mesh geometry
> e = basix.ufl.element("Lagrange", "point", 0, shape=(points.shape\[1\],))
> c\_el = dolfinx.fem.coordinate\_element(e.basix\_element)
> geometry = dolfinx.mesh.create\_geometry(imap, cells, c\_el.\_cpp\_object, points, igi)
> 
> # Create DOLFINx mesh
> if points.dtype == np.float64:
> cpp\_mesh = dolfinx.cpp.mesh.Mesh\_float64(comm, topology, geometry.\_cpp\_object)
> elif points.dtype == np.float32:
> cpp\_mesh = dolfinx.cpp.mesh.Mesh\_float32(comm, topology, geometry.\_cpp\_object)
> else:
> raise RuntimeError(f"Unsupported dtype for mesh {points.dtype}")
> # Wrap as Python object
> return dolfinx.mesh.Mesh(cpp\_mesh, domain = ufl.Mesh(e))
> 
> \`\`\`

where there is no ghosting in construction.  
This is also used (with ghosting) in the internals of DOLFINx MPC, where one adds additional ghost information between processes depending on the multi-point constraint and partitioning.

To get the original input numbering of cells and indices, one uses `original_cell_index` and `input_global_indices`, as for instance sketched out in: [How to apply a boundary condition to a range of DOFs? - #7 by dokken](https://fenicsproject.discourse.group/t/how-to-apply-a-boundary-condition-to-a-range-of-dofs/15758/7)

Let me know if this helps your understanding.