On 08/28/2013 06:08 PM, Matthew Knepley wrote:
On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon
<[email protected]
<mailto:[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]
<mailto:[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.
Hello,
This is the 2D systems I want to simulate:
system 1:
0=\div \rho(x) \grad u + r(x)*u(1-u)
corresponding to the stationary state of the time dependent problem:
system 2:
\partial_t u = \div \rho(x) \grad u + r(x)*u(1-u)
It represents the diffusion of a specie throw a 2D space, u is the
density of this specie.
I want also to study the effect of a predator p:
system 3:
\partial_t u = \div \rho_u(x) \grad u + r(x)*u(1-u) - \beta (x) p*u
\partial_t p = \div \rho_p(x) \grad p - \delta (x) p + \gamma (x) p*u
I'm focused on the (system 1).
About the geometry:
The geometry come from the landscape composed of crops. There are
different type of crops. The functions ( \rho(x), r(x), \delta (x),
\beta (x)) depend on this type of crops.
I'm focused on the (system 1). I'm working from the ex12.c. Therefore,
the plungins functions are:
f0_u(u,grad_u,x,f0){
...
f0= r(x)*u*(1-u)
...
}
f1_u(u,grad_u,x,f1){
...
f1[comp*dim+d] = rho(x)*gradU[comp*dim+d];
...
}
g0_u(u,grad_u,x,g0){
...
g0=-r(x)*(1-2*u)
...
}
g3_uu(u,grad_u,x,g3){
...
g3[((compI*Ncomp+compI)*dim+d)*dim+d] = \rho(x);
...
}
For an efficient implementation of theses plugins, I have to know the
type of crop. If I well understand your previous mail, you propose to me
to defined my own struct PetscFEM adding my useful parameters (crop type
for example). I have to overload the functions
DMPlexComputeResidualFEM, DMPlexComputeJacobianFEM,
FEMIntegrateResidualBatch, FEMIntegrateJacobianActionBatch. I agree,
thanks a lot.
Regards,
Olivier B
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
<tel:%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 <tel:%2B33%20%280%294%2032%2072%2021%2058>
--
Olivier Bonnefon
INRA PACA-Avignon, Unité BioSP
Tel:+33 (0)4 32 72 21 58 <tel:%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
