On Tue, Sep 1, 2015 at 6:15 PM, John Peterson <jwpeter...@gmail.com> wrote:
>
>
> On Tue, Sep 1, 2015 at 8:29 AM, Julian Andrej <j...@tf.uni-kiel.de> wrote:
>
>> I'm trying to get into libMesh, so i tried to convert an existing project.
>>
>> So far i could reproduce the lid driven cavity case for the Stokes
>> equations like in
>> example examples/systems_of_equations/systems_of_equations_ex1.
>>
>> Now i'm trying to build a mass matrix so i can solve the generalized
>> eigenproblem in order to retrieve the physical eigenvalues in the
>> following
>> domain:
>>
>> build_square (mesh,
>> 16, 16,
>>
>> -1., 1.,
>>
>> -1., 1.,
>>
>> TRI6);
>>
>> So the eigenvalues should be (rounded to int) [13, 23, 23, 32, 38, 41...].
>>
>> My approach is that i add a matrix
>>
>> system.add_matrix("mass");
>>
>> in the assemble loop i add:
>>
>> DenseMatrix<Number> Me;
>> Me.resize (n_dofs, n_dofs);
>>
>> and to the integration parts:
>> Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
>>
>> then i close the matrices
>> system.get_matrix("mass").close();
>>
>> and convert them to PETSc so i can use my own SLEPc object
>> Mat A = (static_cast<PetscMatrix<Number> *>(system.matrix))->mat();
>>
>> Mat M = (static_cast<PetscMatrix<Number>
>> &>(system.get_matrix("mass"))).mat();
>>
>> The problem is, the calculated eigenvalues are not the ones i expect (so
>> they're wrong).
>> This method worked for me for the poisson equation.
>> Do i miss something?
>>
>
> So are you trying to solve the generalized eigenvalue problem
>
> A*x = lambda*B*x,
>
> where A corresponds to the Stokes operator and B to the "mass" matrix,
> which should have a zero block for the pressure equation?
>
Exactly
> If yes, then
> .) Are you building both matrices, since you only talk about the mass
> matrix?
>
Yes i build the Stiffness matrix based on the example code
> .) Are you building the correct mass matrix for Stokes?
>
That is the main question i guess. I'm adding the code
Me(i,j) += JxW[qp]*(phi[i][qp]*phi[j][qp]);
to the Kuu and Kvv loops. I've also tried adding it to the Kup, Kvp, Kpv,
Kvp, as well but the quadrature seems to be zero there (using mat view from
petsc and viewing the matrix structure).
>
> Another question is BCs: what are they, and how are you imposing them
> (e.g. penalty method, constraints, etc.)? I can't remember if BCs should
> be imposed in both the A and B matrices or not...
>
I'm imposing pure zero dirichlet using the penalty method as in the example
code. I tried applying them just to A and also tried A and M.
You can take a look at what i've got in that paste
http://hastebin.com/raw/ekasibazan. I don't know if it is desired that i'd
attach the code (i could do so for further reference). It is basically the
example code with a few changes.
>
> --
> John
>
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
