divergenz matrix computed as vector

On 07/20/11 13:46, Melanie Jahny wrote:
> New question #165449 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/165449
>
> I need to compute the following matrix C containing the divergence of trial functions:
>
> C_{ij} = \int_{K_j} div(\psi_i) dx
>
> K_j are the triangles of the trianglulation and \psi_i are the basis functions
> of the raviart-thomas function space.
> In this definition, C_ij should be a matrix, but when I'm computing this,
> I only get a Vector-object and not a matrix.
> (I need to compute the QR-decomposition of this matrix)
>
> I implemented:
>
> mesh = Rectangle(0,0,1,1,10,10)
> mesh.order()
>
> U = FunctionSpace(mesh, "Raviart-Thomas", 1)
> psi = TrialFunction(U)
>
> int_psi = div(psi)*dx
> mat_psi = assemble(int_psi)
>
> The code works, but mat_psi is a Vector - object and so it is not possible to compute
> the QR-decomposition.

It is a vector because that is what you specify: the above reads as

     \int_{triangulation} div (\psi_i) dx

for i = 1, ..., n.

The following might be more along the lines of what you want:

-----
from dolfin import *

mesh = Rectangle(0,0,1,1,10,10)

U = FunctionSpace(mesh, "Raviart-Thomas", 1)
psi = TestFunction(U) # TestFunction has lowest index -> rows

# Define space of piecewise constants over mesh. This space is spanned
# by a set of basis functions, each basis function is associated with
# a cell T and _happens_ to be 1 on T and 0 else where:
Q = FunctionSpace(mesh, "DG", 0)
v = TrialFunction(Q) # TrialFunction has second lowest index -> columns

# Because of comments above this will evaluate as
# \int_{T_h} v_j * div(psi_i) dx = \int_{K_j} div(psi_i)
int_psi = v*div(psi)*dx
mat_psi = assemble(int_psi)

print "Matrix dimensions: ", (mat_psi.size(0), mat_psi.size(1))
---------

--
Marie

> My problem is that I'm note sure if I compute the matrix correctly or if I do have to
> specify the TrialFunctions in another way? Perhaps this is a stupid question but
> I can't find my mistake!
> Any advice would be helpful! Thanks!
>
>
>
>
>