Assemble an optimality system

Hi Ole,

I have done something similar using cbc.block (
https://launchpad.net/cbc.block).
I attach the code.

Kent

On 19 January 2012 11:35, Ole Elvetun
<email address hidden>wrote:

> New question #185175 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/185175
>
> Hi,
>
> I'm working on a problem in optimization with pde constraints, but I have
> some difficulties in trying to implement it.
> Computing the stiffness and mass matrix in FEniCS works fine, but I need
> to collect them in an optimality system.
>
> Some specifically, what I need is a 3x3 block matrix consisting of some
> blocks of stiffness and mass matrices
> (as well as their transpose), a block of the identity matrix as well as
> some blocks of a zero matrix.
> Any idea how to do this in FEniCS?
>
> I have made a code in Octave which solves the problem above in general, so
> I may dump the matrices to file,
> and solve it there, but my guess is that I could solve this faster staying
> in Python.
>
> However, when I try to dump the file, I get the following error :
> "Could not convert _00c2470200000000_p_boost__shared_ptrT_dolfin__Matrix_t
> (type <type 'SwigPyObject'>) to array"
>
> The code leading up to the error is:
> a = (inner(grad(v), grad(u)) - v*u)*dx
> L = v*g*dx
>
> A = assemble(a)
> b = assemble(L)
> scipy.io.savemat('Ab.mat', {'A': A, 'b': b})
>
>
> Also, I wonder if there is an easy way to drop specific columns from, say
> the mass matrix ?
> (in Octave I just did "B = B(:,[x != 0])" to drop the columns where a
> corresponding vector x was equal to zero.)
>
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

Here is a routine for saving to octave (when your not
having to large matrices)

def dump_matrix(filename, name, AA):
    f = open(filename, 'w')
    AA = AA.array()
    for i in range(AA.shape[0]):
        for j in range(AA.shape[1]):
            f.write("%s (%d, %d) = %e;\n " % (name,i+1,j+1,AA[i,j]))

Kent

On 19 January 2012 11:46, Kent-Andre Mardal <email address hidden> wrote:

> Hi Ole,
>
> I have done something similar using cbc.block (
> https://launchpad.net/cbc.block).
> I attach the code.
>
> Kent
>
>
> On 19 January 2012 11:35, Ole Elvetun <
> <email address hidden>> wrote:
>
>> New question #185175 on DOLFIN:
>> https://answers.launchpad.net/dolfin/+question/185175
>>
>> Hi,
>>
>> I'm working on a problem in optimization with pde constraints, but I have
>> some difficulties in trying to implement it.
>> Computing the stiffness and mass matrix in FEniCS works fine, but I need
>> to collect them in an optimality system.
>>
>> Some specifically, what I need is a 3x3 block matrix consisting of some
>> blocks of stiffness and mass matrices
>> (as well as their transpose), a block of the identity matrix as well as
>> some blocks of a zero matrix.
>> Any idea how to do this in FEniCS?
>>
>> I have made a code in Octave which solves the problem above in general,
>> so I may dump the matrices to file,
>> and solve it there, but my guess is that I could solve this faster
>> staying in Python.
>>
>> However, when I try to dump the file, I get the following error :
>> "Could not convert
>> _00c2470200000000_p_boost__shared_ptrT_dolfin__Matrix_t (type <type
>> 'SwigPyObject'>) to array"
>>
>> The code leading up to the error is:
>> a = (inner(grad(v), grad(u)) - v*u)*dx
>> L = v*g*dx
>>
>> A = assemble(a)
>> b = assemble(L)
>> scipy.io.savemat('Ab.mat', {'A': A, 'b': b})
>>
>>
>> Also, I wonder if there is an easy way to drop specific columns from, say
>> the mass matrix ?
>> (in Octave I just did "B = B(:,[x != 0])" to drop the columns where a
>> corresponding vector x was equal to zero.)
>>
>>
>>
>> --
>> You received this question notification because you are a member of
>> DOLFIN Team, which is an answer contact for DOLFIN.
>>
>
>

It seems that you forgot to attach the code for using cbc. /stig

It is in the file ole.py in the first mail. But you have to download
cbc.block to make it work.

Kent

On 19 January 2012 12:10, Stig Larsson <<email address hidden>
> wrote:

> Question #185175 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/185175
>
> Stig Larsson posted a new comment:
> It seems that you forgot to attach the code for using cbc. /stig
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.
>

It is not possible to attach files to the Launchpad answer tracker.

Thanks Kent,

that Octave dump file worked fine! I can't find the example file though. However, the cbc-block seems promising for my purpose, so I might try to have a look at that first.

Ole

Instead of assembling the parts and putting them together as blocks,
you can also define the whole problem using a mixed function space.
Then you assemble the whole thing at once.

--
Anders

On Thu, Jan 19, 2012 at 10:35:40AM -0000, Ole Elvetun wrote:
> New question #185175 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/185175
>
> Hi,
>
> I'm working on a problem in optimization with pde constraints, but I have some difficulties in trying to implement it.
> Computing the stiffness and mass matrix in FEniCS works fine, but I need to collect them in an optimality system.
>
> Some specifically, what I need is a 3x3 block matrix consisting of some blocks of stiffness and mass matrices
> (as well as their transpose), a block of the identity matrix as well as some blocks of a zero matrix.
> Any idea how to do this in FEniCS?
>
> I have made a code in Octave which solves the problem above in general, so I may dump the matrices to file,
> and solve it there, but my guess is that I could solve this faster staying in Python.
>
> However, when I try to dump the file, I get the following error :
> "Could not convert _00c2470200000000_p_boost__shared_ptrT_dolfin__Matrix_t (type <type 'SwigPyObject'>) to array"
>
> The code leading up to the error is:
> a = (inner(grad(v), grad(u)) - v*u)*dx
> L = v*g*dx
>
> A = assemble(a)
> b = assemble(L)
> scipy.io.savemat('Ab.mat', {'A': A, 'b': b})
>
>
> Also, I wonder if there is an easy way to drop specific columns from, say the mass matrix ?
> (in Octave I just did "B = B(:,[x != 0])" to drop the columns where a corresponding vector x was equal to zero.)
>
>
>