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


122
views
0
4 months ago by
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)





Community: FEniCS Project
Please use formatting for your code.
written 4 months ago by Garth Wells  

5 Answers


0
4 months ago by
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:
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
4 months ago by
""" 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
4 months ago by
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:
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 4 months ago by Baltasar  
0
4 months ago by
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
written 4 months ago by Baltasar  
0
4 months ago by
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.

Similar posts:
Search »