<html><head><meta http-equiv="Content-Type" content="text/html charset=utf-8"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space;" class="">Hi Lawrence, <div class=""><br class=""></div><div class="">Thank you. Now it works!</div><div class=""><br class=""></div><div class="">I’ve looked at the source code for DMPlexMatSetClosure() and it uses MatSetValues() so I am wondering in which situation should we use MatSetValuesLocal() ? and how to scatter the local matrices assembled with this ?</div><div class=""><br class=""></div><div class=""><br class=""></div><div class="">Nestor</div><div class=""><br class=""></div><div class=""><br class=""><div><blockquote type="cite" class=""><div class="">On Jun 6, 2016, at 2:38 AM, Lawrence Mitchell <<a href="mailto:lawrence.mitchell@imperial.ac.uk" class="">lawrence.mitchell@imperial.ac.uk</a>> wrote:</div><br class="Apple-interchange-newline"><div class=""><br class=""><blockquote type="cite" class="">On 6 Jun 2016, at 01:43, Nestor Cerpa <<a href="mailto:ncerpagi@umn.edu" class="">ncerpagi@umn.edu</a>> wrote:<br class=""><br class="">Hello,<br class=""><br class="">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.<br class=""><br class="">I made a simple example to explain the issue :<br class="">- I read a simple mesh (generated by gmsh) with 4 quadrilaterls and generate the corresponding DMPlex<br class="">- I define a scalar field on vertices<br class="">- Vertices on Bottom and Top boundaries are contrained (so I have 3 unknowns)<br class="">- I create a DM Matrix<br class="">- For the purpose of my example during the assembly I define an “elementary matrix” containing only ones<br class=""><br class="">When I run the program with one process I obtain what the matrix is supposed to be :<br class=""><br class="">Mat Object: 1 MPI processes<br class=""> type: seqaij<br class="">row 0: (0, 2.) (2, 2.)<br class="">row 1: (1, 2.) (2, 2.)<br class="">row 2: (0, 2.) (1, 2.) (2, 4.)<br class=""><br class="">But running with two processes gives me this :<br class=""><br class="">Mat Object: 2 MPI processes<br class=""> type: mpiaij<br class="">row 0: (0, 2.) (2, 0.)<br class="">row 1: (1, 2.) (2, 2.)<br class="">row 2: (0, 0.) (1, 2.) (2, 2.)<br class=""><br class="">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.<br class=""><br class=""><br class="">For the assembly I am doing something like this (the complete example is attached) :<br class=""><br class=""><br class=""> !*************************************<br class=""> !--- Matrix Assembly<br class=""> !*************************************<br class=""><br class=""> call DMCreateMatrix(dm,A,ierr)<br class=""> call DMGetLocalToGlobalMapping(dm,ltog,ierr)<br class=""> call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr);<br class=""> call ISLocalToGlobalMappingDestroy(ltog,ierr);<br class=""><br class=""><br class="">do icell = !--- Loop over cells<br class=""><br class=""> do ic = ! Loop over closure (somehow only vertices)<br class=""> ipoin = closure((ic-1)*2+1)<br class=""> call DMPlexGetPointLocal(dm,ipoin,pStart,pEnd,ierr)<br class=""> !<br class=""> ! Something to get location of point data in a local Vec<br class=""> !<br class=""> enddo<br class=""><br class=""> elstf(:,:) = 1.0<br class=""><br class=""> do i = !Loop over rows corresponding to dofs which are not constrained<br class=""> call MatSetValuesLocal(A,1,local_ldof(i)-1,ndof*nodes,local_ldof(:)-1, elstf(i,:),ADD_VALUES,ierr)<br class=""> enddo<br class=""><br class="">enddo<br class=""><br class=""> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)<br class=""> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)<br class=""><br class=""> !*************************************<br class=""> !--- End Assembly<br class=""> !*************************************<br class=""><br class=""><br class=""><br class="">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.<br class=""></blockquote><br class=""><br class="">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.<br class=""><br class="">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.<br class=""><br class="">So in pseudo-code you'd probably do:<br class=""><br class="">for c in cells:<br class=""> DMPlexVecGetClosure(dm, NULL, coeff, c, &csize, &values)<br class=""> ! build element matrix (using values)<br class=""> ...<br class=""> DMPlexVecRestoreClosure(dm, NULL, coeff, c, &size, &values)<br class=""> DMPlexMatSetClosure(dm, NULL, NULL, mat, c, element_matrix, ADD_VALUES)<br class=""><br class="">Cheers,<br class=""><br class="">Lawrence<br class=""><br class=""></div></blockquote></div><br class=""><div apple-content-edited="true" class="">
<div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;"><br class="Apple-interchange-newline"><br class=""></div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">---------------------------------------------------</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">---------------------------------------------------</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;"><br class=""></div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">Nestor Cerpa Gilvonio</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">Postdoctoral research associate</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;"><br class=""></div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">Department of Earth Sciences</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">University of Minnesota</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">310 Pillsbury Drive SE</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">Minnepolis, MN 55455</div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;"><br class=""></div><div class="" style="color: rgb(0, 0, 0); font-family: Helvetica; font-size: 12px; font-style: normal; font-variant: normal; font-weight: normal; letter-spacing: normal; line-height: normal; orphans: auto; text-align: start; text-indent: 0px; text-transform: none; white-space: normal; widows: auto; word-spacing: 0px; -webkit-text-stroke-width: 0px;">E-mail : <a href="mailto:cerpa003@umn.edu" class="">ncerpagi@umn.edu</a></div>
</div>
<br class=""></div></body></html>