How to restrict a Function to a part of a mesh given by a MeshFunction?

Asked by Heinz Zorn on 2013-06-24

Hello,

I would like to restrict an existing function to only a part of the mesh. The wanted part is given by a MeshFunction. Simply project the function to a submesh doesn't work in parallel. I did some experimentation with Restriction() but this led to nothing.

Here a minimal example to explain what I'm trying to do:

from dolfin import *

mesh = UnitCubeMesh( 10, 10, 10 )

V0 = VectorFunctionSpace( mesh, "DG", 0 )
v0 = project( Constant( ( 1, 1, 1 ) ), V0 )
plot( v0 )

domains = CellFunction( "size_t", mesh )
for cell in cells( mesh ):
 p = cell.midpoint()
 if p.x() < p.y():
  domains[ cell ] = 1
 else:
  domains[ cell ] = 2
plot( domains )

v0_vec = v0.vector().array()
for cell in cells( mesh ):
 if domains[ cell ] == 1:
  print( cell ) # dummy line to make the code executable
  # v0_vec[ k ] = 0 for the correct k

v0.vector().set_local( v0_vec )

interactive()

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Heinz Zorn
Solved:
2013-06-24
Last query:
2013-06-24
Last reply:
Heinz Zorn (zorn) said : #1

And of course thanks in advance for any help!

Heinz Zorn (zorn) said : #2

I found the right indices:
  v0_vec[ 3*cell.index() ] = 0
  v0_vec[ 3*cell.index()+1 ] = 0
  v0_vec[ 3*cell.index()+2 ] = 0