### How to use, as a finite element function, a vector obtained from a system with matrices formed with cbc.block

233

views

1

I am solving a Poisson equation with a mixed formulation where the second variable is the Hessian matrix. I have managed to use abc.block up to the point where I have a vector x which should contain the degrees of freedom of the approximate solution. I need to compare it with the exact solution.

The tricks used in the examples, to convert a vector x to a finite element function, don't seem to work (u0.vector() = x[:], u0 = map(Function, V,x ) or Function(V,u)).

I am probably not taking into account some property of the objects being manipulated in the code. Could you suggest a fix? Perhaps I can even get suggestions for a better alternative.

I am using dolfin version 2017.1.0 with fences installed via condo-forge.

A minimal code is copied below.

""" Poisson with jump terms P1-P0 Discontinuous discrete Hessian implemented

with unknown the scalar variable"""

from __future__ import division

from dolfin import *

from block import *

from block.dolfin_util import *

from block.iterative import *

from block.algebraic.petsc import *

import numpy as np

import os

dolfin.set_log_level(15) #

n=4;

mesh = UnitSquareMesh(n, n)

V = FunctionSpace(mesh, "CG", 1)

Sigma = TensorFunctionSpace(mesh, "DG", 0)

exact = Expression('exp((x[0]*x[0]+x[1]*x[1])/2.0)',degree=5)

f = Expression('(2+x[0]*x[0]+x[1]*x[1])*exp((x[0]*x[0]+x[1]*x[1])/2)',degree=5)

class Boundary(SubDomain): # define the Dirichlet boundary

def inside(self, x, on_boundary):

return on_boundary

exact_boundary = Boundary()

bc = DirichletBC(V, exact, exact_boundary)

bcs = block_bc([None, bc], True)

v, u = TestFunction(V), TrialFunction(V)

tau, sigma = TestFunction(Sigma), TrialFunction(Sigma)

n = FacetNormal(mesh)

gu=grad(u)

a11 = inner(sigma, tau)*dx

a12 = avg(tau[0,0])*(gu[0]('+')*n[0]('+')+gu[0]('-')*n[0]('-'))*dS\

+ avg(tau[0,1])*(gu[1]('+')*n[0]('+')+gu[1]('-')*n[0]('-'))*dS \

+avg(tau[1,0])*(gu[0]('+')*n[1]('+')+gu[0]('-')*n[1]('-'))*dS \

+avg(tau[1,1])*(gu[1]('+')*n[1]('+')+gu[1]('-')*n[1]('-'))*dS

a21 = v*(sigma[0,0]+sigma[1,1])*dx

L2 = f * v * dx

AA = block_assemble([[a11, a12],

[a21, 0 ]])

bb = block_assemble([0, L2])

bcs.apply(AA).apply(bb)

# Extract the individual submatrices

[[A, B],

[C, _]] = AA

[G, F] = bb

Ainv = ConjGrad(A, tolerance=1e-10, show=2)

R = C*Ainv*B

L=C*A

u0 = Function(V)

sigma0 = Function(Sigma)

Rinv = ConjGrad(R, tolerance=1e-10, show=2)

x=Rinv*L

#print x.shape

#u0.vector() = x[:]

#u0 = map(Function, V,x )

#Function(V,u)

vertex_values_exact = exact.compute_vertex_values(mesh)

The tricks used in the examples, to convert a vector x to a finite element function, don't seem to work (u0.vector() = x[:], u0 = map(Function, V,x ) or Function(V,u)).

I am probably not taking into account some property of the objects being manipulated in the code. Could you suggest a fix? Perhaps I can even get suggestions for a better alternative.

I am using dolfin version 2017.1.0 with fences installed via condo-forge.

A minimal code is copied below.

""" Poisson with jump terms P1-P0 Discontinuous discrete Hessian implemented

with unknown the scalar variable"""

from __future__ import division

from dolfin import *

from block import *

from block.dolfin_util import *

from block.iterative import *

from block.algebraic.petsc import *

import numpy as np

import os

dolfin.set_log_level(15) #

n=4;

mesh = UnitSquareMesh(n, n)

V = FunctionSpace(mesh, "CG", 1)

Sigma = TensorFunctionSpace(mesh, "DG", 0)

exact = Expression('exp((x[0]*x[0]+x[1]*x[1])/2.0)',degree=5)

f = Expression('(2+x[0]*x[0]+x[1]*x[1])*exp((x[0]*x[0]+x[1]*x[1])/2)',degree=5)

class Boundary(SubDomain): # define the Dirichlet boundary

def inside(self, x, on_boundary):

return on_boundary

exact_boundary = Boundary()

bc = DirichletBC(V, exact, exact_boundary)

bcs = block_bc([None, bc], True)

v, u = TestFunction(V), TrialFunction(V)

tau, sigma = TestFunction(Sigma), TrialFunction(Sigma)

n = FacetNormal(mesh)

gu=grad(u)

a11 = inner(sigma, tau)*dx

a12 = avg(tau[0,0])*(gu[0]('+')*n[0]('+')+gu[0]('-')*n[0]('-'))*dS\

+ avg(tau[0,1])*(gu[1]('+')*n[0]('+')+gu[1]('-')*n[0]('-'))*dS \

+avg(tau[1,0])*(gu[0]('+')*n[1]('+')+gu[0]('-')*n[1]('-'))*dS \

+avg(tau[1,1])*(gu[1]('+')*n[1]('+')+gu[1]('-')*n[1]('-'))*dS

a21 = v*(sigma[0,0]+sigma[1,1])*dx

L2 = f * v * dx

AA = block_assemble([[a11, a12],

[a21, 0 ]])

bb = block_assemble([0, L2])

bcs.apply(AA).apply(bb)

# Extract the individual submatrices

[[A, B],

[C, _]] = AA

[G, F] = bb

Ainv = ConjGrad(A, tolerance=1e-10, show=2)

R = C*Ainv*B

L=C*A

u0 = Function(V)

sigma0 = Function(Sigma)

Rinv = ConjGrad(R, tolerance=1e-10, show=2)

x=Rinv*L

#print x.shape

#u0.vector() = x[:]

#u0 = map(Function, V,x )

#Function(V,u)

vertex_values_exact = exact.compute_vertex_values(mesh)

Community: FEniCS Project

Please use formatting for your code.

written
10 months ago by
Garth Wells

### 5 Answers

0

I did not go through your code (see https://www.allanswered.com/post/emre/please-read-how-to-get-your-question-answered/ ), but if you need to set the values of a vector, you need to add [:] behind it:

If you will do that in parallel, you have to use

`u0.vector()[:] = x`

I think, you don't even need the [:] behind the x, if it is an numpy array or a dolfin vector.If you will do that in parallel, you have to use

`u0.vector().set_local(your_array)`

where your array has to be in vertex ordering, not nodal ordering. (Or is it the other way around? I'm always confused with this..)
0

```
""" Poisson with jump terms P1-P0 Discontinuous discrete Hessian implemented
with unknown the scalar variable"""
from __future__ import division
from dolfin import *
from block import *
from block.dolfin_util import *
from block.iterative import *
from block.algebraic.petsc import *
import numpy as np
import os
dolfin.set_log_level(15) #
n=4;
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG", 1)
Sigma = TensorFunctionSpace(mesh, "DG", 0)
exact = Expression('exp((x[0]*x[0]+x[1]*x[1])/2.0)',degree=5)
f = Expression('(2+x[0]*x[0]+x[1]*x[1])*exp((x[0]*x[0]+x[1]*x[1])/2)',degree=5)
class Boundary(SubDomain): # define the Dirichlet boundary
def inside(self, x, on_boundary):
return on_boundary
exact_boundary = Boundary()
bc = DirichletBC(V, exact, exact_boundary)
bcs = block_bc([None, bc], True)
v, u = TestFunction(V), TrialFunction(V)
tau, sigma = TestFunction(Sigma), TrialFunction(Sigma)
n = FacetNormal(mesh)
gu=grad(u)
a11 = inner(sigma, tau)*dx
a12 = avg(tau[0,0])*(gu[0]('+')*n[0]('+')+gu[0]('-')*n[0]('-'))*dS\
+ avg(tau[0,1])*(gu[1]('+')*n[0]('+')+gu[1]('-')*n[0]('-'))*dS \
+avg(tau[1,0])*(gu[0]('+')*n[1]('+')+gu[0]('-')*n[1]('-'))*dS \
+avg(tau[1,1])*(gu[1]('+')*n[1]('+')+gu[1]('-')*n[1]('-'))*dS
a21 = v*(sigma[0,0]+sigma[1,1])*dx
L2 = f * v * dx
AA = block_assemble([[a11, a12],
[a21, 0 ]])
bb = block_assemble([0, L2])
bcs.apply(AA).apply(bb)
# Extract the individual submatrices
[[A, B],
[C, _]] = AA
[G, F] = bb
Ainv = ConjGrad(A, tolerance=1e-10, show=2)
R = C*Ainv*B
L=C*A
u0 = Function(V)
sigma0 = Function(Sigma)
Rinv = ConjGrad(R, tolerance=1e-10, show=2)
x=Rinv*L
#print x.shape
#u0.vector() = x[:]
#u0 = map(Function, V,x )
#Function(V,u)
vertex_values_exact = exact.compute_vertex_values(mesh)
```

0

Thanks for your replies.

In the previous post, I have formatted the code.

Using u0.vector()[:] = x, I get the error

raise IndexError("can only set full slices v[:]")
Hey,

This happens for me, if I try to set a part of the dolfin vector:

Are you sure, x and u0 have the same shape? What type does x have?

In the previous post, I have formatted the code.

Using u0.vector()[:] = x, I get the error

raise IndexError("can only set full slices v[:]")

This happens for me, if I try to set a part of the dolfin vector:

`u0.vector()[0:20] =some_np_array[0:20] # Do not do that, it won't work`

Setting the values works fine in the following example:```
import fenics as fnx
import numpy as np
mesh = fnx.UnitSquareMesh(4,4)
V = fnx.FunctionSpace(mesh, "CG", 1)
u0 = fnx.Function(V)
#create numpy array with same shape
u0_shape = u0.vector().array().shape[0]
u0_numpy = np.arange(0, u0_shape)
u0.vector()[:] = u0_numpy
```

Are you sure, x and u0 have the same shape? What type does x have?

written
10 months ago by
Baltasar

0

Thank you very much.

x has type block.block_compose.block_mul

and u0.vector has type dolfin.cpp.la.GenericVector

No wonder. But does this suggest a way forward?
Hey,

I think I cannot help you here. I have not worked with the block package. Further things, I would look into:

1) Have you checked if it is possible to set the dolfin vector class with the block_mul class?

2) Additionally, have you checked the shapes? Perhaps you have to retrieve a numpy array out of the block_mul class to be able to set the dolfin vector.

Good luck

Baltasar

x has type block.block_compose.block_mul

and u0.vector has type dolfin.cpp.la.GenericVector

No wonder. But does this suggest a way forward?

I think I cannot help you here. I have not worked with the block package. Further things, I would look into:

1) Have you checked if it is possible to set the dolfin vector class with the block_mul class?

2) Additionally, have you checked the shapes? Perhaps you have to retrieve a numpy array out of the block_mul class to be able to set the dolfin vector.

Good luck

Baltasar

written
10 months ago by
Baltasar

0

Thank you very much. I will look into that later. For now, I will just code it directly. Again, thank you.

Please login to add an answer/comment or follow this question.