<p dir="ltr">I am sorry about the delay. I have your example working but it exposed a bug in Plex so I need to push the fix first. I should have everything for you early next week.</p>
<p dir="ltr"> Thanks</p>
<p dir="ltr"> Matt</p>
<div class="gmail_quote">On Feb 25, 2016 2:04 AM, "Julian Andrej" <<a href="mailto:juan@tf.uni-kiel.de">juan@tf.uni-kiel.de</a>> wrote:<br type="attribution"><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">After a bit of rethinking the problem, the discrepancy between the size of matrix A and the mass matrix M arises because of the Dirichlet boundary conditions. So why aren't the BCs not imposed on the mass matrix? Do I need to handle Dirichlet BCs differently in this context (like zero rows and put one the diagonal?)<br>
<br>
On <a href="tel:24.02.2016%2020" value="+12402201620" target="_blank">24.02.2016 20</a>:54, juan wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
I attached another example which creates the correct mass matrix<br>
but also overwrites the DM for the SNES solve. Somehow i cannot manage<br>
to really copy the DM to dm_mass and use that. If i try to do that with<br>
DMClone(dm, &dm_mass) i get a smaller mass matrix (which is not of size A).<br>
<br>
Maybe this helps in the discussion.<br>
<br>
Relevant code starts at line 455.<br>
<br>
On 2016-02-24 15:03, Julian Andrej wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Thanks Matt,<br>
<br>
I attached the modified example.<br>
<br>
the corresponding code (and only changes to ex12) is starting at line<br>
832.<br>
<br>
It also seems that the mass matrix is of size 169x169 and the<br>
stiffness matrix is of dimension 225x225. I'd assume that if i<br>
multiply test and trial function i'd get a matrix of same size (if the<br>
space/quadrature is the same for the stiffness matrix)<br>
<br>
On <a href="tel:24.02.2016%2014" value="+12402201614" target="_blank">24.02.2016 14</a>:56, Matthew Knepley wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
On Wed, Feb 24, 2016 at 7:47 AM, Julian Andrej <<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a><br>
<mailto:<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a>>> wrote:<br>
<br>
I'm now using the petsc git master branch.<br>
<br>
I tried adding my code to the ex12<br>
<br>
DM dm_mass;<br>
PetscDS prob_mass;<br>
PetscFE fe;<br>
Mat M;<br>
PetscFECreateDefault(dm, user.dim, 1, PETSC_TRUE, NULL, -1, &fe);<br>
<br>
DMClone(dm, &dm_mass);<br>
DMGetDS(dm_mass, &prob_mass);<br>
PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);<br>
PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL, NULL,<br>
NULL);<br>
DMCreateMatrix(dm_mass, &M);<br>
<br>
MatSetOptionsPrefix(M, "M_";)<br>
<br>
and receive the error on running<br>
./exe -interpolate -refinement_limit 0.0125 -petscspace_order 2<br>
-M_mat_view binary<br>
<br>
WARNING! There are options you set that were not used!<br>
WARNING! could be spelling mistake, etc!<br>
Option left: name:-M_mat_view value: binary<br>
<br>
I don't know if the matrix is actually there and assembled or if the<br>
option is ommitted because something is wrong.<br>
<br>
<br>
Its difficult to know when I cannot see the whole code. You can always<br>
insert<br>
<br>
MatViewFromOptions(M, NULL, "-mat_view");<br>
<br>
Using<br>
MatView(M, PETSC_VIEWER_STDOUT_WORLD);<br>
<br>
gives me a reasonable output to stdout.<br>
<br>
<br>
Good.<br>
<br>
But saving the matrix and analysing it in matlab, results in an all<br>
zero matrix.<br>
<br>
PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Mout",FILE_MODE_WRITE,<br>
&viewer);<br>
MatView(M, viewer);<br>
<br>
<br>
I cannot explain this, but it has to be something like you are viewing<br>
the matrix before it is<br>
actually assembled. Feel free to send the code. It sounds like it is<br>
mostly working.<br>
<br>
Matt<br>
<br>
Any hints?<br>
<br>
<br>
On <a href="tel:24.02.2016%2013" value="+12402201613" target="_blank">24.02.2016 13</a> <tel:24.02.2016%2013>:58, Matthew Knepley wrote:<br>
<br>
On Wed, Feb 24, 2016 at 6:47 AM, Julian Andrej<br>
<<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a> <mailto:<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a>><br>
<mailto:<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a> <mailto:<a href="mailto:juan@tf.uni-kiel.de" target="_blank">juan@tf.uni-kiel.de</a>>>><br>
wrote:<br>
<br>
Hi,<br>
<br>
i'm trying to assemble a mass matrix with the<br>
PetscFE/DMPlex<br>
interface. I found something in the examples of TAO<br>
<br>
<a href="https://bitbucket.org/petsc/petsc/src/da8116b0e8d067e39fd79740a8a864b0fe207998/src/tao/examples/tutorials/ex3.c?at=master&fileviewer=file-view-default" rel="noreferrer" target="_blank">https://bitbucket.org/petsc/petsc/src/da8116b0e8d067e39fd79740a8a864b0fe207998/src/tao/examples/tutorials/ex3.c?at=master&fileviewer=file-view-default</a><br>
<br>
<br>
but using the lines<br>
<br>
DMClone(dm, &dm_mass);<br>
DMSetNumFields(dm_mass, 1);<br>
DMPlexCopyCoordinates(dm, dm_mass);<br>
DMGetDS(dm_mass, &prob_mass);<br>
PetscDSSetJacobian(prob_mass, 0, 0, mass_kernel, NULL,<br>
NULL, NULL);<br>
PetscDSSetDiscretization(prob_mass, 0, (PetscObject) fe);<br>
DMPlexSNESComputeJacobianFEM(dm_mass, u, M, M, NULL);<br>
DMCreateMatrix(dm_mass, &M);<br>
<br>
leads to errors in DMPlexSNESComputeJacobianFEM (u is a<br>
global vector).<br>
<br>
I don't can understand the necessary commands until<br>
DMPlexSNESComputeJacobianFEM. What does it do and why is it<br>
necessary? (especially why does the naming involve SNES?)<br>
<br>
Is there another/easier/better way to create a mass<br>
matrix (the<br>
inner product of the function space and the test space)?<br>
<br>
<br>
1) That example needs updating. First, look at SNES ex12 which<br>
is up to<br>
date.<br>
<br>
2) I assume you are using 3.6. If you use the development<br>
version, you<br>
can remove DMPlexCopyCoordinates().<br>
<br>
3) You need to create the matrix BEFORE calling the assembly<br>
<br>
4) Always always always send the entire error messge<br>
<br>
Matt<br>
<br>
Regards<br>
Julian Andrej<br>
<br>
<br>
<br>
<br>
--<br>
What most experimenters take for granted before they begin their<br>
experiments is infinitely more interesting than any results<br>
to which<br>
their experiments lead.<br>
-- Norbert Wiener<br>
<br>
<br>
<br>
<br>
--<br>
What most experimenters take for granted before they begin their<br>
experiments is infinitely more interesting than any results to which<br>
their experiments lead.<br>
-- Norbert Wiener<br>
</blockquote></blockquote></blockquote>
</blockquote></div>