Manipulate Mixed Functions

Asked by Marcos Samudio

I am trying to run this simple code:

mesh = UnitSquare(20,20)

V = VectorFunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
W = MixedFunctionSpace([V,Q])

u, p = TrialFunctions(W)
v, q = TestFunctions(W)

up_ = Function(W)
u_, p_ = up_.split()

Suppose I perform one iteration of some variational form, and get my results stored in up_.
What I want to do next is to pick and use some values of u_, so I convert it to a vector instance:

U_ = u_.vector()

But I get the following message:

Size of vector: 1323
Size of function space: 882
---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)

I am actually doing u_.vector() before solving, but this should work anyway, right?

Thanks for your time,

Marcos

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Marcos Samudio
Solved:
Last query:
Last reply:
Revision history for this message
Marcos Samudio (marcossamudio) said :
#1

Just for the record, doing

u_ = Function(Q)
U_ = u_.vector()

works fine, but I would like to use the mixed formulation.

Thanks!!

Revision history for this message
Johan Hake (johan-hake) said :
#2

If you are only going to use the values for some calculations, and not
change the corresponding values in the mixed function you can cal split
with deepcopy=True:

   u_, p_ = up_.split(deepcopy=True)

Now the vector of u_ is a copy of the vector corresponding to the first
  subfunction of up.

Johan

On 04/20/2012 10:30 PM, Marcos Samudio wrote:
> New question #194259 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/194259
>
> I am trying to run this simple code:
>
> mesh = UnitSquare(20,20)
>
> V = VectorFunctionSpace(mesh, 'CG', 1)
> Q = FunctionSpace(mesh, 'CG', 1)
> W = MixedFunctionSpace([V,Q])
>
> u, p = TrialFunctions(W)
> v, q = TestFunctions(W)
>
> up_ = Function(W)
> u_, p_ = up_.split()
>
> Suppose I perform one iteration of some variational form, and get my results stored in up_.
> What I want to do next is to pick and use some values of u_, so I convert it to a vector instance:
>
> U_ = u_.vector()
>
> But I get the following message:
>
> Size of vector: 1323
> Size of function space: 882
> ---------------------------------------------------------------------------
> RuntimeError Traceback (most recent call last)
>
> I am actually doing u_.vector() before solving, but this should work anyway, right?
>
> Thanks for your time,
>
> Marcos
>
>
>

Revision history for this message
Marcos Samudio (marcossamudio) said :
#3

Hello Johan,

Well, I am actually thinking on modifying some values as well. How could that be achieved?

When using Dirichlet boundary conditions this is done internally by FEniCS, right? How could I do it manually? I am trying to set some wall values in a flow simulation.

Thanks

Revision history for this message
Anders Logg (logg) said :
#4

On Sat, Apr 21, 2012 at 04:25:48PM -0000, Marcos Samudio wrote:
> Question #194259 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/194259
>
> Marcos Samudio posted a new comment:
> Hello Johan,
>
> Well, I am actually thinking on modifying some values as well. How could
> that be achieved?
>
> When using Dirichlet boundary conditions this is done internally by
> FEniCS, right? How could I do it manually? I am trying to set some wall
> values in a flow simulation.

The DirichletBC class does something like this:

  collect values to set in array <values>
  collect positions to set in array <dofs>
  call b.set(values, size, dofs)
  call b.apply("insert")

For details, check the file DirichletBC.cpp. It's quite complex but
well documented. You will need to work throught he DofMap object of a
(sub)FunctionSpace to get the positions in the vector to modify.

The vector to modify will be the vector for the entire mixed space.

--
Anders

Revision history for this message
Marcos Samudio (marcossamudio) said :
#5

Hello Anders,

I will take a look at the DirichletBC class. Anyway I think I have figured a way to do it using only Python..
 I would just do

Revision history for this message
Marcos Samudio (marcossamudio) said :
#6

Sorry for that, I accidentally pressed the comment button.

I changed my mind and will do it differently, I will modify the up_ vector instead, which should not be too hard, since the velocities would always be the first mesh.geometry().dim()*N dofs of each node, where N is the number of nodes (at least when using CG1 elements, I guess...)

So my code is something like:

for i in vertices_on_boundary:
      j = bndary_to_inner_mapping[i]
      for k in range(mesh.geometry().dim()):
             up_.vector()[i + k*N] = some_function(up_.vector()[j +k*N])

Anyway, if something seems wrong please let me know!