Mesheditor in MPI


239
views
1
7 months ago by
Hello,

In serial, to create a mesh with the mesheditor i do (it works in Fenics 2016.1):

Node_ = spio.loadmat(Node_Electrolyte_path)['node_region'] # node_region is the variable name saved in Matlab
Tetrahedron_ = spio.loadmat(Tetrahedron_Electrolyte_path)['elem_'] # elem_ is the variable name saved in Matlab

# Number of vertices
number_vertices=len(Node_)
# Number of cells
number_cells=len(Tetrahedron_)

# Create an empty mesh object
mesh_ = Mesh();
# Call the editor
editor = MeshEditor();
# Open the mesh (mesh to edit, topological dimension, geometrical dimension)
editor.open(mesh_, 3, 3)
# Set number of vertices
editor.init_vertices(number_vertices)
# Set number of cells
editor.init_cells(number_cells)

# Set vortex location (vortex id,x,y)
for vertex_id in range(number_vertices):
	x_ = Node_[vertex_id][0]        # Vertex coordinate
	y_ = Node_[vertex_id][1]        # Vertex coordinate
	z_ = Node_[vertex_id][2]        # Vertex coordinate
	# The vertex id must be an integer
	vertex_id=int(vertex_id)
	# Add the vertex
	editor.add_vertex(vertex_id, x_, y_, z_)

# Set cell (cell id, vertices id)
for cell_id in range(number_cells):
	v0 = Tetrahedron_[cell_id][0] # Vertex indice
	v1 = Tetrahedron_[cell_id][1] # Vertex indice
	v2 = Tetrahedron_[cell_id][2] # Vertex indice
	v3 = Tetrahedron_[cell_id][3] # Vertex indice
	# The id must be an integer
	cell_id=int(cell_id)
	v0=int(v0); v1=int(v1);	v2=int(v2); v3=int(v3)
	# Add the cell
	editor.add_cell(cell_id, v0, v1, v2, v3)

# Close the editor
editor.close()

plot(mesh_, title='mesh', interactive=True)

​


Now, how am i suppose to do the same thing with MPI instructions ?
It seems possible to build a mesh in MPI, with the mesheditor (with MeshEditor.init_cells_global, MeshEditor.init_vertices_global MeshEditor.add_vertex_global plus MeshEditor.add_cell), but i didnt find examples.

I have tried:

# Number of vertices
number_vertices=len(Node_)
# Number of cells
number_cells=len(Tetrahedron_)

# Verteces distribution among process
verteces_process = np.linspace(0, number_vertices-1, num=number_process+1)
if rank_process==0:
	range_verteces=[int(verteces_process[0]), int(verteces_process[1])]
else:
	range_verteces=[int(verteces_process[rank_process])+1, int(verteces_process[rank_process+1])]
number_vertices_process=range_verteces[1]-range_verteces[0]+1

# Distribute cell among processor
cell_process = np.linspace(0, number_cells-1, num=number_process+1)
if rank_process==0:
	range_cell=[int(cell_process[0]), int(cell_process[1])]
else:
	range_cell=[int(cell_process[rank_process])+1, int(cell_process[rank_process+1])]	
number_cell_process=range_cell[1]-range_cell[0]+1

# Create an empty mesh object
mesh_ = Mesh();
# Call the editor
editor = MeshEditor();
# Open the mesh (mesh to edit, topological dimension, geometrical dimension)
editor.open(mesh_, 3, 3)
# Set number of vertices
editor.init_vertices_global(number_vertices_process,number_vertices)
# Set number of cells
editor.init_cells_global(number_cell_process,number_cells)
	
mpi_comm.Barrier() # Blocks until all processes in the communicator have reached this routine 
	
# Set vortex location (vortex local index, vortex global index, point)
vertex_id_local=-1
global_to_local=np.zeros((number_vertices,1)) # Global index to local index
global_to_local_allreduce=np.zeros((number_vertices,1))
for vertex_global_id in range(range_verteces[0], range_verteces[1]):
	vertex_id_local+=1
	global_to_local[vertex_global_id][0]=vertex_id_local # Global to local index
	x_ = Node_Electrolyte[vertex_global_id][0]        # Vertex coordinate
	y_ = Node_Electrolyte[vertex_global_id][1]        # Vertex coordinate
	z_ = Node_Electrolyte[vertex_global_id][2]        # Vertex coordinate
	# Create point
	p=Point(x_, y_, z_)
	# The vertex id must be an integer
	vertex_global_id=int(vertex_global_id)
	# Add the vertex
	editor.add_vertex_global(vertex_id_local,vertex_global_id,p)
	
mpi_comm.Barrier() # Blocks until all processes in the communicator have reached this routine 
mpi_comm.Allreduce(global_to_local,global_to_local_allreduce, op=MPI.SUM)	
mpi_comm.Barrier() # Blocks until all processes in the communicator have reached this routine 

# Set cell (cell local index, cell global index, vertices array)
cell_id_local=0
for cell_id_global in range(range_cell[0], range_cell[1]):
	cell_id_local+=1
	v0 = Tetrahedron_Electrolyte[cell_id_global][0] # Vertex global indice
	v1 = Tetrahedron_Electrolyte[cell_id_global][1] # Vertex global indice
	v2 = Tetrahedron_Electrolyte[cell_id_global][2] # Vertex global indice
	v3 = Tetrahedron_Electrolyte[cell_id_global][3] # Vertex global indice
	vo_local=global_to_local_allreduce[v0][0] # Vertex local indice
	v1_local=global_to_local_allreduce[v1][0] # Vertex local indice
	v2_local=global_to_local_allreduce[v2][0] # Vertex local indice
	v3_local=global_to_local_allreduce[v3][0] # Vertex local indice
	
	# The id must be an integer
	cell_id_global=int(cell_id_global)
	v0=int(vo_local); v1=int(v1_local);	v2=int(v2_local); v3=int(v3_local)
	v_=[v0, v1, v2, v3]
	
	# Add the cell
	editor.add_cell(cell_id_local, cell_id_global, v_)

mpi_comm.Barrier() # Blocks until all processes in the communicator have reached this routine 

# Close the editor
editor.close()

mpi_comm.Barrier() # Blocks until all processes in the communicator have reached this routine 

plot(mesh_, title='mesh_, process {}'.format(rank_process), interactive=True)	


My code crashed at the line: editor.add_cell(cell_id_local, cell_id_global, v_) with the error 'TypeError: (size_t) expected positive 'int' for argument 4'

Two questions:
- I think the way i distributed the verteces and cell for each process is not correct. Do you have any guidance here?
- What is wrong with the add_cell instruction in the code?

Thanks.




Community: FEniCS Project
Check a related thread at issue tracker. Might be helpful for you https://bitbucket.org/fenics-project/dolfin/issues/403/mesheditor-and-meshinit-not-initializing.

Regarding the typecheck problem, try numpy.array([...], dtype=numpy.uintp).
written 7 months ago by Jan Blechta  
Thanks for the answer.
So, it seems the way it works is one process builds the whole mesh, then it is partitioned, as shown in this example i have taken from the link you have sent me:

#create the mesh by mesheditor 
m=Mesh()

if rank==0 :
        editor = MeshEditor()
        editor.open(m,2,2)
        editor.init_vertices_global(4,4)
        editor.init_cells_global(2,2)
        editor.add_vertex_global(0,0, Point(0.0,0.0))                
        editor.add_vertex_global(1,1, Point(1.0,0.0))                
        editor.add_vertex_global(2,2, Point(1.0,1.0))                
        editor.add_vertex_global(3,3, Point(0.0,1.0))
        editor.add_cell(0,0,np.array([0,1,2],dtype=np.uintp))
        editor.add_cell(1,1,np.array([0,2,3],dtype=np.uintp))                
        editor.close()
        m.init()
        m.order()
        info(m)

build_distributed_mesh(m)
info(m)
plot(m,interactive=True)​

written 7 months ago by Francois Usseglio Viretta  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »