<html>
<head>
<meta content="text/html; charset=ISO-8859-1"
http-equiv="Content-Type">
</head>
<body bgcolor="#FFFFFF" text="#000000">
<div class="moz-cite-prefix">On 08/28/2013 06:08 PM, Matthew Knepley
wrote:<br>
</div>
<blockquote
cite="mid:CAMYG4GniHr+omf7ru4-BUnrNPzkNbu8iX0UzxeOgZRrn45Za_w@mail.gmail.com"
type="cite">
<div dir="ltr">On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon
<span dir="ltr"><<a moz-do-not-send="true"
href="mailto:olivier.bonnefon@avignon.inra.fr"
target="_blank">olivier.bonnefon@avignon.inra.fr</a>></span>
wrote:<br>
<div class="gmail_extra">
<div class="gmail_quote">
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
<div>
<div class="h5">
<div>On 08/28/2013 10:28 AM, Olivier Bonnefon wrote:<br>
</div>
<blockquote type="cite">
<div>On 08/23/2013 04:42 PM, Matthew Knepley
wrote:<br>
</div>
<blockquote type="cite">
<div dir="ltr">On Fri, Aug 23, 2013 at 9:35 AM,
Olivier Bonnefon <span dir="ltr"><<a
moz-do-not-send="true"
href="mailto:olivier.bonnefon@avignon.inra.fr"
target="_blank">olivier.bonnefon@avignon.inra.fr</a>></span>
wrote:<br>
<div class="gmail_extra">
<div class="gmail_quote">
<blockquote class="gmail_quote"
style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Hello,<br>
<br>
Thanks for your answers, I'm now able to
import and distribute a mesh:<br>
</blockquote>
<div><br>
</div>
<div>You might simplify this to</div>
<div><br>
</div>
<div>if (rank) {obNbCells = 0; obNbVertex
= 0;}</div>
<div>ierr =
DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);<br>
</div>
<div> </div>
<blockquote class="gmail_quote"
style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">if
(!rank){<br>
ierr =
DMPlexCreateFromCellList(comm,dim,obNbCells,obNbVertex,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);<br>
for (i=0;i<obNbBound;i++){<br>
ierr
=DMPlexSetLabelValue(*dm, "marker",
obBoundary[i]+obNbCells,
1);CHKERRQ(ierr);<br>
}<br>
}else {<br>
ierr =
DMPlexCreateFromCellList(comm,dim,0,0,3,0,obCells,2,obVertex,dm);CHKERRQ(ierr);<br>
}<br>
<br>
ierr = DMPlexDistribute(*dm,
partitioner, 0,
&distributedMesh);CHKERRQ(ierr);<br>
if (distributedMesh) {<br>
ierr =
DMDestroy(dm);CHKERRQ(ierr);<br>
*dm = distributedMesh;<br>
}<br>
<br>
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)?<br>
I need this to write an efficient
implementation of the FEM struct
functions f0 and g0, space depending.<br>
</blockquote>
<div><br>
</div>
<div>Yes, but I really do not think you
want to do things that way. I am
assuming you want different material
models or something</div>
<div>in different places. The way I
envision that is using a DMLabel to mark
up parts of the domain. All labels are
automatically</div>
<div>distributed with the mesh. Is that
what you want?</div>
</div>
</div>
</div>
</blockquote>
Hello,<br>
<br>
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:<br>
<br>
for (i=0;i<obNbCells;i++){<br>
if (labelCells[i]==1){<br>
ierr =DMPlexSetLabelValue(*dm, "marker1",
i, 1);CHKERRQ(ierr);<br>
}else{<br>
ierr =DMPlexSetLabelValue(*dm, "marker2",
i, 1);CHKERRQ(ierr);<br>
}<br>
}<br>
<br>
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 ?</blockquote>
</div>
</div>
</div>
</blockquote>
<div><br>
</div>
<div>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</div>
<div>equation, so I do a loop over all cells. This loop
takes place in DMPlexComputeResidualFEM(). You would
instead like</div>
<div>to do a few loops over sets of cells with different
material models, using different f0/f1. Is this correct?</div>
<div> </div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000">
<div>
<div class="h5">
<blockquote type="cite"> Thanks a lot for your help.<br>
<br>
Olivier B<br>
</blockquote>
</div>
</div>
Hello,<br>
<br>
This is the solution I implemented to get the label
level in the plugins "fem.f0Funcs" and "fem.g0Funcs":<br>
<br>
I need the DM and the index element, so i do:<br>
1) I add some static variables:<br>
static DM * spDM[128];<br>
static int scurElem[128];<br>
</div>
</blockquote>
<div><br>
</div>
<div>Notice that the DM is available in
DMPlexComputeResidualFEM(). Here is what the function
does:</div>
<div><br>
</div>
<div> a) Batches up elements into groups</div>
<div><br>
</div>
<div> b) Integrates each group using a call to
FEMIntegrateResidualBatch(). Notice that in 'next' this
has</div>
<div> changed to PetscFEIntegrateResidual() since we
have added a few FEM classes to make things</div>
<div> simpler and more flexible.</div>
<div><br>
</div>
<div>What you can do, I think, to get what you want is:</div>
<div><br>
</div>
<div> a) Write a new MY_DMPlexComputeResidualFEM() to do a
few loops. This is supplied to your app using</div>
<div><br>
</div>
<div>
<div> ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode
(*)(DM,Vec,Vec,void*)) MY_DMPlexComputeResidualFEM,
&user);CHKERRQ(ierr);</div>
<div> ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode
(*)(DM,Vec,Mat,Mat,MatStructure*,void*))
MY_DMPlexComputeJacobianFEM, &user);CHKERRQ(ierr);</div>
</div>
<div><br>
</div>
<div>just as in the examples. You could use different f0/f1
for each loop somehow.</div>
<div><br>
</div>
<div> b) Write a new PetscFEIntegrateResidual() that does
what you want. The easiest way to do this is create</div>
<div> a new PetscFE subclass, since they only really
do one thing which is these integrals. I can help you.</div>
<div> </div>
<div>HOWEVER, if what you really want to do is get
coefficient information into f0/f1 instead of a different
physical model,</div>
<div>then you can do something easier that we just put in.
You can layout a coefficient, like nu in</div>
<div><br>
</div>
<div> \div \nu \grad u = \rho</div>
<div><br>
</div>
<div>and provide a DM for \nu. This will be passed all the
way down inside until f0 gets</div>
<div><br>
</div>
<div> f0Func(u, gradU, nu, gradNu, x, f0)</div>
<div><br>
</div>
<div>so that the pointwise values of the your coefficient
and its gradient are available to your physics.</div>
<div><br>
</div>
<div>I am sure there will be questions about this, but the
first thing to do is get entirely clear what you want to
do.</div>
</div>
</div>
</div>
</blockquote>
Hello,<br>
<br>
This is the 2D systems I want to simulate:<br>
<br>
system 1:<br>
0=\div \rho(x) \grad u + r(x)*u(1-u) <br>
<br>
corresponding to the stationary state of the time dependent problem:<br>
system 2:<br>
\partial_t u = \div \rho(x) \grad u + r(x)*u(1-u) <br>
<br>
It represents the diffusion of a specie throw a 2D space, u is the
density of this specie. <br>
<br>
I want also to study the effect of a predator p:<br>
system 3:<br>
\partial_t u = \div \rho_u(x) \grad u + r(x)*u(1-u) - \beta (x) p*u<br>
\partial_t p = \div \rho_p(x) \grad p - \delta (x) p + \gamma (x)
p*u<br>
<br>
I'm focused on the (system 1).<br>
<br>
About the geometry:<br>
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. <br>
<br>
I'm focused on the (system 1). I'm working from the ex12.c.
Therefore, the plungins functions are:<br>
<br>
f0_u(u,grad_u,x,f0){<br>
...<br>
f0= r(x)*u*(1-u)<br>
...<br>
}<br>
<br>
f1_u(u,grad_u,x,f1){<br>
...<br>
f1[comp*dim+d] = rho(x)*gradU[comp*dim+d];<br>
...<br>
}<br>
<br>
g0_u(u,grad_u,x,g0){<br>
...<br>
g0=-r(x)*(1-2*u)<br>
...<br>
}<br>
<br>
g3_uu(u,grad_u,x,g3){<br>
...<br>
g3[((compI*Ncomp+compI)*dim+d)*dim+d] = \rho(x);<br>
...<br>
}<br>
<br>
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.<br>
<br>
Regards,<br>
<br>
Olivier B<br>
<blockquote
cite="mid:CAMYG4GniHr+omf7ru4-BUnrNPzkNbu8iX0UzxeOgZRrn45Za_w@mail.gmail.com"
type="cite">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div><br>
</div>
<blockquote class="gmail_quote" style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">
<div bgcolor="#FFFFFF" text="#000000"> 2) I overload the
DMPlexComputeJacobianFEM with :<br>
PetscErrorCode MY_DMPlexComputeJacobianFEM(DM dm, Vec
X, Mat Jac, Mat JacP, MatStructure *str,void *user)<br>
{<br>
<br>
PetscMPIInt rank;<br>
PetscErrorCode ierr;<br>
<br>
PetscFunctionBeginUser;<br>
ierr = MPI_Comm_rank(PETSC_COMM_WORLD,
&rank);CHKERRQ(ierr);<br>
spDM[rank]=&dm;<br>
PetscFunctionReturn(DMPlexComputeJacobianFEM(dm, X,
Jac,JacP, str,user));<br>
<br>
}<br>
3) overload FEMIntegrateResidualBatch adding code:<br>
.<br>
.<br>
for (e = 0; e < Ne; ++e) {<br>
scurElem[rank]=e;//added ligne<br>
.<br>
.<br>
<br>
So that, I can get the label level using DMPlexHasLabel
and DMLabelGetValue<br>
<br>
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 ??<br>
<br>
Thanks,
<div class="im"><br>
Olivier B<br>
<blockquote type="cite">
<blockquote type="cite">
<div dir="ltr">
<div class="gmail_extra">
<div class="gmail_quote">
<div><br>
</div>
<div> Thanks,</div>
<div><br>
</div>
<div> Matt</div>
<div> </div>
<blockquote class="gmail_quote"
style="margin:0px 0px 0px
0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex">Regards,<br>
<br>
Olivier B<span><font color="#888888"><br>
<br>
-- <br>
Olivier Bonnefon<br>
INRA PACA-Avignon, Unité BioSP<br>
Tel: <a moz-do-not-send="true"
href="tel:%2B33%20%280%294%2032%2072%2021%2058"
value="+33432722158" target="_blank">+33
(0)4 32 72 21 58</a><br>
<br>
</font></span></blockquote>
</div>
<br>
<br clear="all">
<div><br>
</div>
-- <br>
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>
</blockquote>
<br>
<br>
<pre cols="72">--
Olivier Bonnefon
INRA PACA-Avignon, Unité BioSP
Tel: <a moz-do-not-send="true" href="tel:%2B33%20%280%294%2032%2072%2021%2058" value="+33432722158" target="_blank">+33 (0)4 32 72 21 58</a></pre>
</blockquote>
<br>
<br>
<pre cols="72">--
Olivier Bonnefon
INRA PACA-Avignon, Unité BioSP
Tel: <a moz-do-not-send="true" href="tel:%2B33%20%280%294%2032%2072%2021%2058" value="+33432722158" target="_blank">+33 (0)4 32 72 21 58</a></pre>
</div>
</div>
</blockquote>
</div>
<br>
<br clear="all">
<div><br>
</div>
-- <br>
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>
</blockquote>
<br>
<br>
<pre class="moz-signature" cols="72">--
Olivier Bonnefon
INRA PACA-Avignon, Unité BioSP
Tel: +33 (0)4 32 72 21 58</pre>
</body>
</html>