# Access a vector of a subfunction

Hi,

I'm trying to access a vector of a subfunction to be used by normalize().

V = VectorFunctionS

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

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

W = MixedFunctionSp

w = Function(W)

normalize(

AttributeError: 'Indexed' object has no attribute 'vector'

normalize(

TypeError: expected an int, slice, list or numpy array of integers

#W.sub(

(u, p, theta) = w.split()

normalize(

*** Error: Unable to access vector of degrees of freedom fro function.

*** Reason: Cannot access a non-const vector from a subfunction.

*** Where: This error was encountered inside Function.cpp.

I've found few another ways. One could for example deepcopy p and then interpolate back to w. But it's an extra interpolate instead of simple algebraic operation and it could reduce accuraccy of normalization.

Or it has been suggested here that one could get indeces of w.vector() belonging to p-part by V.dim(), Q.dim() or by num_vertices(). But it will not work in parallel and also num_vertices work-around depends on chosen FunctionSpaces.

Is there any clean way to do it? Thanks

Jan Blechta

## Question information

- Language:
- English Edit question

- Status:
- Answered

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Last query:

- Last reply:

Hi Jan,

I used to have the same problem and my solution can be found in the extended_normalize class that is part of CBC.PDESys. I have copied it here for your convenience. It works in parallel as well.

from dolfin import *

class extended_normalize:

"""Normalize part or whole of vector.

V = Functionspace we normalize in

u = Function where part is normalized

part = The index of the part of the mixed function space

that we want to normalize.

For example. When solving for velocity and pressure coupled in the

Navier-Stokes equations we sometimes (when there is only Neuman BCs

on pressure) need to normalize the pressure.

Example of use:

mesh = UnitSquare(1, 1)

V = VectorFunctionS

Q = FunctionSpace(mesh, 'CG', 1)

VQ = V * Q

up = Function(VQ)

normalize_func = extended_

up.vector()[:] = 2.

print 'before ', up.vector(

normalize_

print 'after ', up.vector(

results in:

before [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

after [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0]

"""

def __init__(self, V, part='entire vector'):

self.part = part

if isinstance(part, int):

self.u = Function(V)

v = TestFunction(V)

self.c = assemble(

self.pp = ['0']*self.

self.u0 = interpolate(

self.x0 = self.u0.vector()

self.C1 = assemble(

def __call__(self, v):

if isinstance(

# assemble into c1 the part of the vector that we want to normalize

c1 = self.C1.inner(v)

if abs(c1) > 1.e-8:

# Perform normalization

else:

# normalize entire vector

dummy = normalize(v)

mesh = UnitSquare(2, 2)

V = VectorFunctionS

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

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

W = MixedFunctionSp

w = Function(W)

w.vector()[:] = 2.

print 'before ', w.vector(

normalize = extended_

normalize(

print 'after ', w.vector(

Hope it does the trick for you as well,

Best regards

Mikael

Den Aug 7, 2012 kl. 3:55 PM skrev Jan Blechta:

> New question #205221 on DOLFIN:

> https:/

>

> Hi,

> I'm trying to access a vector of a subfunction to be used by normalize().

>

>

> V = VectorFunctionS

> Q = FunctionSpace(mesh, "CG", 1)

> R = FunctionSpace(mesh, "CG", 1)

> W = MixedFunctionSp

>

> w = Function(W)

>

> normalize(

> AttributeError: 'Indexed' object has no attribute 'vector'

>

> normalize(

> TypeError: expected an int, slice, list or numpy array of integers

> #W.sub(

>

> (u, p, theta) = w.split()

> normalize(

> *** Error: Unable to access vector of degrees of freedom fro function.

> *** Reason: Cannot access a non-const vector from a subfunction.

> *** Where: This error was encountered inside Function.cpp.

>

>

> I've found few another ways. One could for example deepcopy p and then interpolate back to w. But it's an extra interpolate instead of simple algebraic operation and it could reduce accuraccy of normalization.

>

> Or it has been suggested here that one could get indeces of w.vector() belonging to p-part by V.dim(), Q.dim() or by num_vertices(). But it will not work in parallel and also num_vertices work-around depends on chosen FunctionSpaces.

>

> Is there any clean way to do it? Thanks

> Jan Blechta

>

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

Jan Blechta (blechta) said : | #2 |

Ok, thanks, this seems being a solution. I'll try it in few days..

But I think the problem is more general - there is no secure way to simply access and manipulate dofs of a subfunction, isn't it?

Den Aug 12, 2012 kl. 1:35 PM skrev Jan Blechta:

> Question #205221 on DOLFIN changed:

> https:/

>

> Status: Answered => Open

>

> Jan Blechta is still having a problem:

> Ok, thanks, this seems being a solution. I'll try it in few days..

>

> But I think the problem is more general - there is no secure way to

> simply access and manipulate dofs of a subfunction, isn't it?

Not sure it's the most efficient solution, but you can extract sub dofmaps by doing something like:

from dolfin import *

from numpy import array

parameters[

# Create mesh and finite element

mesh = UnitSquare(2, 2)

VV = VectorFunctionS

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

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

VQ = V * Q

# Create a function to normalize in

u = Function(VQ)

u.vector()[:] = 2.

# Get the subdofmap

sub_space = VQ.sub(1)

dofmap = sub_space.dofmap()

# Extract a list of dofs

dofs = dofmap.

# Remove off process dofs

off_dofs = dofmap.

for i in off_dofs.keys():

dofs.remove(i)

# Sum values on local process

total_sum = u.vector()[dofs]

total_sum = total_sum.sum()

# The above segfults for PETSc. The code below works in that case

#total_sum = 0

#for i in dofs:

# total_sum += u.vector()[i]

# Subtract from the original vector

v = Vector(u.vector())

v.zero()

v[dofs] = -MPI.sum(total_sum) / Q.dim()

print MPI.process_

u.vector(

print MPI.process_

For me this results in:

[mikael@ubuntu:~]$ mpirun -np 2 python test_dofs.py

Process 0: Number of global vertices: 9

Process 0: Number of global cells: 8

Process 1: Partitioned mesh, edge cut is 2.

Process 0: Partitioned mesh, edge cut is 2.

0 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2]

1 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

1 After: [2 0 0 0 2 0 2 2 2 2 2 2 2 2 2 2 2 2 0 0]

0 After: [2 2 2 2 2 0 2 0 2 2 0 2 2 2]

[mikael@ubuntu:~]$ mpirun -np 1 python test_dofs.py

0 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

0 After: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0]

Best regards

Mikael

>

> --

> You received this question notification because you are a member of

> DOLFIN Team, which is an answer contact for DOLFIN.

Johan Hake (johan-hake) said : | #4 |

On 08/12/2012 10:30 PM, Mikael Mortensen wrote:

> Question #205221 on DOLFIN changed:

> https:/

>

> Status: Open => Answered

>

> Mikael Mortensen proposed the following answer:

>

> Den Aug 12, 2012 kl. 1:35 PM skrev Jan Blechta:

>

>> Question #205221 on DOLFIN changed:

>> https:/

>>

>> Status: Answered => Open

>>

>> Jan Blechta is still having a problem:

>> Ok, thanks, this seems being a solution. I'll try it in few days..

>>

>> But I think the problem is more general - there is no secure way to

>> simply access and manipulate dofs of a subfunction, isn't it?

>

> Not sure it's the most efficient solution, but you can extract sub

> dofmaps by doing something like:

>

> from dolfin import *

> from numpy import array

> parameters[

>

> # Create mesh and finite element

> mesh = UnitSquare(2, 2)

> VV = VectorFunctionS

> V = FunctionSpace(mesh, "CG", 2)

> Q = FunctionSpace(mesh, "CG", 1)

> VQ = V * Q

>

> # Create a function to normalize in

> u = Function(VQ)

> u.vector()[:] = 2.

>

> # Get the subdofmap

> sub_space = VQ.sub(1)

> dofmap = sub_space.dofmap()

>

> # Extract a list of dofs

> dofs = dofmap.

>

> # Remove off process dofs

> off_dofs = dofmap.

> for i in off_dofs.keys():

> dofs.remove(i)

For larger data sets it might be faster to use sets:

dofs = set(dofmap.

dofs.

dofs = list(dofs)

Johan

> # Sum values on local process

> total_sum = u.vector()[dofs]

> total_sum = total_sum.sum()

> # The above segfults for PETSc. The code below works in that case

> #total_sum = 0

> #for i in dofs:

> # total_sum += u.vector()[i]

>

> # Subtract from the original vector

> v = Vector(u.vector())

> v.zero()

> v[dofs] = -MPI.sum(total_sum) / Q.dim()

> print MPI.process_

> u.vector(

> print MPI.process_

>

>

> For me this results in:

> [mikael@ubuntu:~]$ mpirun -np 2 python test_dofs.py

> Process 0: Number of global vertices: 9

> Process 0: Number of global cells: 8

> Process 1: Partitioned mesh, edge cut is 2.

> Process 0: Partitioned mesh, edge cut is 2.

> 0 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2]

> 1 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

> 1 After: [2 0 0 0 2 0 2 2 2 2 2 2 2 2 2 2 2 2 0 0]

> 0 After: [2 2 2 2 2 0 2 0 2 2 0 2 2 2]

> [mikael@ubuntu:~]$ mpirun -np 1 python test_dofs.py

> 0 Before: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2]

> 0 After: [2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 0 0 0 0 0 0 0 0 0]

>

>

> Best regards

>

> Mikael

>

>>

>> --

>> You received this question notification because you are a member of

>> DOLFIN Team, which is an answer contact for DOLFIN.

>

## Can you help with this problem?

Provide an answer of your own, or ask Jan Blechta for more information if necessary.