<div dir="ltr">On Wed, Aug 28, 2013 at 10:01 AM, Olivier Bonnefon <span dir="ltr"><<a 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 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><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 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 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 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>