On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon < [email protected]> 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 < > [email protected]> 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
