I am aware of that, but I think i should get a bunch of 1 eigenvalues
(number depending on the problem size and thus number of Dirichlet BCs) and
then the correct eigenvalues for my problem with the posted code.
On Tue, Sep 1, 2015 at 6:38 PM, David Knezevic <david.kneze...@akselos.com>
wrote:
> On Tue, Sep 1, 2015 at 12:35 PM, Julian Andrej <j...@tf.uni-kiel.de>
> wrote:
>
>> 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.
>>
>
>
> Regarding BCs: Constraint rows introduce spurious eigenvalues. For
> Dirichlet boundary conditions, I would recommend following the approach in
> eigenproblems_ex3 which uses a CondensedEigenSystem, which eliminates
> Dirichlet dofs from the system before doing the solve.
>
> David
>
>
------------------------------------------------------------------------------
_______________________________________________
Libmesh-users mailing list
Libmesh-users@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/libmesh-users
