DOLFIN

Assemble diagonal matrix

Asked by Claas Abert

Hi,

I'd like to solve a PDE with constraints on the nodes of the mesh.
The resulting matrix would look like this

|A B|
|B^t 0|

where A is the result of a fenics assembly A=assemble(a) and B is diagonal.
Is there a way/what is the best way to assemble a diagonal matrix in FEniCS?

So far my idea would be to
a) use a mixed formulation, but I couldn't figure out how to formulate the constraint on the nodes (without integration)
b) use CBC.Block and the diag command, but the resulting B seems to be stored as dense numpy matrix

Any hints are appreciated.

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Joachim Haga
Solved:
Last query:
Last reply:
Revision history for this message
Best Joachim Haga (jobh) said : 		#1

I can comment on the cbc.block issue only.

> b) use CBC.Block and the diag command, but the resulting B seems to be
stored as dense numpy matrix

There's no diag command in cbc.block, possibly you are hitting numpy.diag.

cbc.block can create a diagonal matrix, but only when using the trilinos
backend. Two ways, both require a Dolfin/trilinos vector (with the correct
dimension and parallel layout) as input:

block.algebraic.trilinos.create_identity(v) -- creates a real identity
matrix (the matrix elements may be manipulated afterwards)

block.algebraic.trilinos.diag_op(down_cast(v).vec()) -- wraps the vector so
that it acts like a diagonal matrix in matrix-vector products etc.

Revision history for this message
Claas Abert (cabert) said : 		#2

Thanks Joachim Haga, that solved my question.

Revision history for this message
Claas Abert (cabert) said : 		#3

Thank you Joachim Haga,

this looks like the perfect solution for me. The only problem was that
down_cast(v).vec() returns an Epetra.Epetra_FEVector which is not
accepted by diag_op.

I'm now using:
diag_op(Epetra.MultiVector(down_cast(v).vec()))
which seems to do the job.

Revision history for this message
Joachim Haga (jobh) said : 		#4

Good! I'll fix diag_op so that it accepts Epetra.Epetra_FEVector and
dolfin.Vector.

One thing to be aware of is that if you need boundary conditions, you will
have to find and zero the corresponding rows of B yourself. If it's a real
matrix instead, you can use DirichletBC.zero. If this is a concern, the
easiest way to convert to a real matrix is:

d = diag_op(...)
I = create_identity(...)
D = collapse(d*I)

-j.

On 7 June 2012 16:45, Claas Abert <email address hidden>wrote:

> Question #199682 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/199682
>
> Claas Abert posted a new comment:
> Thank you Joachim Haga,
>
> this looks like the perfect solution for me. The only problem was that
> down_cast(v).vec() returns an Epetra.Epetra_FEVector which is not
> accepted by diag_op.
>
> I'm now using:
> diag_op(Epetra.MultiVector(down_cast(v).vec()))
> which seems to do the job.
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Revision history for this message
Claas Abert (cabert) said : 		#5

Right now I don't have boundary conditions, but thanks for the hint!
Claas