> On 6 Jun 2016, at 01:43, Nestor Cerpa <ncerp...@umn.edu> wrote: > > Hello, > > I am learning to use Petsc in order to modify my sequential FEM code to work > in parallel. First, I am trying to understand how a FEM assembly loop would > work in parallel using Petsc but I run into some problems. > > I made a simple example to explain the issue : > - I read a simple mesh (generated by gmsh) with 4 quadrilaterls and generate > the corresponding DMPlex > - I define a scalar field on vertices > - Vertices on Bottom and Top boundaries are contrained (so I have 3 unknowns) > - I create a DM Matrix > - For the purpose of my example during the assembly I define an “elementary > matrix” containing only ones > > When I run the program with one process I obtain what the matrix is supposed > to be : > > Mat Object: 1 MPI processes > type: seqaij > row 0: (0, 2.) (2, 2.) > row 1: (1, 2.) (2, 2.) > row 2: (0, 2.) (1, 2.) (2, 4.) > > But running with two processes gives me this : > > Mat Object: 2 MPI processes > type: mpiaij > row 0: (0, 2.) (2, 0.) > row 1: (1, 2.) (2, 2.) > row 2: (0, 0.) (1, 2.) (2, 2.) > > It seems that in the way I am doing the assembly, Process 0 is not giving the > contribution of its local assembly to process 1. > > > For the assembly I am doing something like this (the complete example is > attached) : > > > !************************************* > !--- Matrix Assembly > !************************************* > > call DMCreateMatrix(dm,A,ierr) > call DMGetLocalToGlobalMapping(dm,ltog,ierr) > call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr); > call ISLocalToGlobalMappingDestroy(ltog,ierr); > > > do icell = !--- Loop over cells > > do ic = ! Loop over closure (somehow only vertices) > ipoin = closure((ic-1)*2+1) > call DMPlexGetPointLocal(dm,ipoin,pStart,pEnd,ierr) > ! > ! Something to get location of point data in a local Vec > ! > enddo > > elstf(:,:) = 1.0 > > do i = !Loop over rows corresponding to dofs which are not > constrained > call > MatSetValuesLocal(A,1,local_ldof(i)-1,ndof*nodes,local_ldof(:)-1, > elstf(i,:),ADD_VALUES,ierr) > enddo > > enddo > > call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) > call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) > > !************************************* > !--- End Assembly > !************************************* > > > > There is something I am definitely missing because I don’t understand why it > does not work. I was thinking maybe I need a function like > DMLocallToGlobal…() but for matrices.
You can either do this "by hand" using DMPlexGetTransitiveClosure to get the points. Or else you use DMPlexMatSetClosure to scatter the element matrix to the appropriate points in the matrix. Note there is also a matching DMPlexVecSetClosure to set values in the residual vector, and DMPlexVecGetClosure to get the values in the closure of a point. So in pseudo-code you'd probably do: for c in cells: DMPlexVecGetClosure(dm, NULL, coeff, c, &csize, &values) ! build element matrix (using values) ... DMPlexVecRestoreClosure(dm, NULL, coeff, c, &size, &values) DMPlexMatSetClosure(dm, NULL, NULL, mat, c, element_matrix, ADD_VALUES) Cheers, Lawrence
signature.asc
Description: Message signed with OpenPGP using GPGMail
