<div dir="ltr"><div class="gmail_extra"><div class="gmail_quote">On Tue, Jun 7, 2016 at 6:38 AM, Nestor Cerpa <span dir="ltr"><<a href="mailto:ncerpagi@umn.edu" target="_blank">ncerpagi@umn.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word">Hi Lawrence, <div><br></div><div>Thank you. Now it works!</div><div><br></div><div>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></blockquote><div><br></div><div>MatSetValuesLocal() just uses a mapping from local indices to global indices, and then calls MatSetValues(). It does not</div><div>assemble a different type of matrix.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div style="word-wrap:break-word"><div>Nestor</div><div><br></div><div><div><div class="h5"><br><div><blockquote type="cite"><div>On Jun 6, 2016, at 2:38 AM, Lawrence Mitchell <<a href="mailto:lawrence.mitchell@imperial.ac.uk" target="_blank">lawrence.mitchell@imperial.ac.uk</a>> wrote:</div><br><div><br><blockquote type="cite">On 6 Jun 2016, at 01:43, Nestor Cerpa <<a href="mailto:ncerpagi@umn.edu" target="_blank">ncerpagi@umn.edu</a>> wrote:<br><br>Hello,<br><br>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><br>I made a simple example to explain the issue :<br>- I read a simple mesh (generated by gmsh) with 4 quadrilaterls and generate the corresponding DMPlex<br>- I define a scalar field on vertices<br>- Vertices on Bottom and Top boundaries are contrained (so I have 3 unknowns)<br>- I create a DM Matrix<br>- For the purpose of my example during the assembly I define an “elementary matrix” containing only ones<br><br>When I run the program with one process I obtain what the matrix is supposed to be :<br><br>Mat Object: 1 MPI processes<br> type: seqaij<br>row 0: (0, 2.) (2, 2.)<br>row 1: (1, 2.) (2, 2.)<br>row 2: (0, 2.) (1, 2.) (2, 4.)<br><br>But running with two processes gives me this :<br><br>Mat Object: 2 MPI processes<br> type: mpiaij<br>row 0: (0, 2.) (2, 0.)<br>row 1: (1, 2.) (2, 2.)<br>row 2: (0, 0.) (1, 2.) (2, 2.)<br><br>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><br><br>For the assembly I am doing something like this (the complete example is attached) :<br><br><br> !*************************************<br> !--- Matrix Assembly<br> !*************************************<br><br> call DMCreateMatrix(dm,A,ierr)<br> call DMGetLocalToGlobalMapping(dm,ltog,ierr)<br> call MatSetLocalToGlobalMapping(A,ltog,ltog,ierr);<br> call ISLocalToGlobalMappingDestroy(ltog,ierr);<br><br><br>do icell = !--- Loop over cells<br><br> do ic = ! Loop over closure (somehow only vertices)<br> ipoin = closure((ic-1)*2+1)<br> call DMPlexGetPointLocal(dm,ipoin,pStart,pEnd,ierr)<br> !<br> ! Something to get location of point data in a local Vec<br> !<br> enddo<br><br> elstf(:,:) = 1.0<br><br> do i = !Loop over rows corresponding to dofs which are not constrained<br> call MatSetValuesLocal(A,1,local_ldof(i)-1,ndof*nodes,local_ldof(:)-1, elstf(i,:),ADD_VALUES,ierr)<br> enddo<br><br>enddo<br><br> call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr)<br> call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr)<br><br> !*************************************<br> !--- End Assembly<br> !*************************************<br><br><br><br>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></blockquote><br><br>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><br>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><br>So in pseudo-code you'd probably do:<br><br>for c in cells:<br> DMPlexVecGetClosure(dm, NULL, coeff, c, &csize, &values)<br> ! build element matrix (using values)<br> ...<br> DMPlexVecRestoreClosure(dm, NULL, coeff, c, &size, &values)<br> DMPlexMatSetClosure(dm, NULL, NULL, mat, c, element_matrix, ADD_VALUES)<br><br>Cheers,<br><br>Lawrence<br><br></div></blockquote></div><br></div></div><div>
<div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px"><br><br></div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">---------------------------------------------------</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">---------------------------------------------------</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px"><br></div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">Nestor Cerpa Gilvonio</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">Postdoctoral research associate</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px"><br></div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">Department of Earth Sciences</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">University of Minnesota</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">310 Pillsbury Drive SE</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">Minnepolis, MN 55455</div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px"><br></div><div 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;text-align:start;text-indent:0px;text-transform:none;white-space:normal;word-spacing:0px">E-mail : <a href="mailto:cerpa003@umn.edu" target="_blank">ncerpagi@umn.edu</a></div>
</div>
<br></div></div></blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature">What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div>
</div></div>