[petsc-users] distribute and cells mapping.
Matthew Knepley
knepley at gmail.com
Wed Aug 28 11:08:59 CDT 2013
On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon <
olivier.bonnefon at avignon.inra.fr> wrote:
> On 08/28/2013 10:28 AM, Olivier Bonnefon wrote:
>
> On 08/23/2013 04:42 PM, Matthew Knepley wrote:
>
> On Fri, Aug 23, 2013 at 9:35 AM, Olivier Bonnefon <
> olivier.bonnefon at avignon.inra.fr> wrote:
>
>> Hello,
>>
>> Thanks for your answers, I'm now able to import and distribute a mesh:
>>
>
> You might simplify this to
>
> if (rank) {obNbCells = 0; obNbVertex = 0;}
> ierr =
> DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>
>
>> if (!rank){
>> ierr =
>> DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>> for (i=0;i<obNbBound;i++){
>> ierr =DMPlexSetLabelValue(*dm, "marker",
>> obBoundary[i]+obNbCells, 1);CHKERRQ(ierr);
>> }
>> }else {
>> ierr =
>> DMPlexCreateFromCellList(comm,dim,0,0,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);
>> }
>>
>> ierr = DMPlexDistribute(*dm, partitioner, 0,
>> &distributedMesh);CHKERRQ(ierr);
>> if (distributedMesh) {
>> ierr = DMDestroy(dm);CHKERRQ(ierr);
>> *dm = distributedMesh;
>> }
>>
>> Is it possible to known the resulting partition ? ie, What is the mapping
>> between the initial cell number and the local cell (used in
>> DMPlexComputeResidualFEM)?
>> I need this to write an efficient implementation of the FEM struct
>> functions f0 and g0, space depending.
>>
>
> Yes, but I really do not think you want to do things that way. I am
> assuming you want different material models or something
> in different places. The way I envision that is using a DMLabel to mark up
> parts of the domain. All labels are automatically
> distributed with the mesh. Is that what you want?
>
> Hello,
>
> It is exactly what I need: I'm mobilized a landscape, and the parameters
> of the model depend of the type of crop. Therefore, I have created a label
> for each type of crop and I have labeled each triangle with the
> corresponding label:
>
> for (i=0;i<obNbCells;i++){
> if (labelCells[i]==1){
> ierr =DMPlexSetLabelValue(*dm, "marker1", i, 1);CHKERRQ(ierr);
> }else{
> ierr =DMPlexSetLabelValue(*dm, "marker2", i, 1);CHKERRQ(ierr);
> }
> }
>
> So, I'm able to mark the triangles, but I'm not able to get this label in
> the plugin "fem.f0Funcs" and "fem.g0Funcs": These plugins are called by
> looping on the triangles in the function "FEMIntegrateResidualBatch", but
> the dm is not available, so I can't use the functions DMPlexGetLabel,
> DMLabelGetStratumSize and DMLabelGetStratumIS. What is the good way to get
> the labels in the user plugins of the fem struct ?
>
>
So lets start with the abstract problem so that I can see exactly what you
want to do. In ex12 (or ex62, etc.) I have a single
equation, so I do a loop over all cells. This loop takes place in
DMPlexComputeResidualFEM(). You would instead like
to do a few loops over sets of cells with different material models, using
different f0/f1. Is this correct?
> Thanks a lot for your help.
>
> Olivier B
>
> Hello,
>
> This is the solution I implemented to get the label level in the plugins
> "fem.f0Funcs" and "fem.g0Funcs":
>
> I need the DM and the index element, so i do:
> 1) I add some static variables:
> static DM * spDM[128];
> static int scurElem[128];
>
Notice that the DM is available in DMPlexComputeResidualFEM(). Here is what
the function does:
a) Batches up elements into groups
b) Integrates each group using a call to FEMIntegrateResidualBatch().
Notice that in 'next' this has
changed to PetscFEIntegrateResidual() since we have added a few FEM
classes to make things
simpler and more flexible.
What you can do, I think, to get what you want is:
a) Write a new MY_DMPlexComputeResidualFEM() to do a few loops. This is
supplied to your app using
ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM,Vec,Vec,void*))
MY_DMPlexComputeResidualFEM, &user);CHKERRQ(ierr);
ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode
(*)(DM,Vec,Mat,Mat,MatStructure*,void*)) MY_DMPlexComputeJacobianFEM,
&user);CHKERRQ(ierr);
just as in the examples. You could use different f0/f1 for each loop
somehow.
b) Write a new PetscFEIntegrateResidual() that does what you want. The
easiest way to do this is create
a new PetscFE subclass, since they only really do one thing which is
these integrals. I can help you.
HOWEVER, if what you really want to do is get coefficient information into
f0/f1 instead of a different physical model,
then you can do something easier that we just put in. You can layout a
coefficient, like nu in
\div \nu \grad u = \rho
and provide a DM for \nu. This will be passed all the way down inside until
f0 gets
f0Func(u, gradU, nu, gradNu, x, f0)
so that the pointwise values of the your coefficient and its gradient are
available to your physics.
I am sure there will be questions about this, but the first thing to do is
get entirely clear what you want to do.
Thanks,
Matt
2) I overload the DMPlexComputeJacobianFEM with :
> PetscErrorCode MY_DMPlexComputeJacobianFEM(DM dm, Vec X, Mat Jac, Mat
> JacP, MatStructure *str,void *user)
> {
>
> PetscMPIInt rank;
> PetscErrorCode ierr;
>
> PetscFunctionBeginUser;
> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank);CHKERRQ(ierr);
> spDM[rank]=&dm;
> PetscFunctionReturn(DMPlexComputeJacobianFEM(dm, X, Jac,JacP, str,user));
>
> }
> 3) overload FEMIntegrateResidualBatch adding code:
> .
> .
> for (e = 0; e < Ne; ++e) {
> scurElem[rank]=e;//added ligne
> .
> .
>
> So that, I can get the label level using DMPlexHasLabel and
> DMLabelGetValue
>
> I'm sure this solution is awful, and works only in this version, but i
> didn't find a better way to get the label in the plugins fem struc. Do you
> know the correct way to do that ??
>
> Thanks,
>
> Olivier B
>
>
> Thanks,
>
> Matt
>
>
>> Regards,
>>
>> Olivier B
>>
>> --
>> Olivier Bonnefon
>> INRA PACA-Avignon, Unité BioSP
>> Tel: +33 (0)4 32 72 21 58 <%2B33%20%280%294%2032%2072%2021%2058>
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
>
>
> --
> Olivier Bonnefon
> INRA PACA-Avignon, Unité BioSP
> Tel: +33 (0)4 32 72 21 58
>
>
>
> --
> Olivier Bonnefon
> INRA PACA-Avignon, Unité BioSP
> Tel: +33 (0)4 32 72 21 58
>
>
--
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20130828/f24cc358/attachment-0001.html>
More information about the petsc-users
mailing list