DOLFIN

Questions in defining function and variational problem

Asked by Houdong Hu

# r, c, rs, cx are the solutions of an eigenvalue solver, what is the datastructure of rx, cx? a vector?
r, c, rx, cx = eigensolver.get_eigenpair(0)

# rho is a new function I defined, I want to set rho = 2rx^2, is "rho.vector()[:] = 2*rx*rx" right?
rho = Function(V)
rho.vector()[:] = 2*rx*rx

# I want to solve "grad^2 V_rho = rho", do I need to interpolate rho first?
V_rho = Function(V)
b = -0.25*pi*inner(nabla_grad(u), nabla_grad(v))*dx
L = rho*v*dx
solve(b == L, V_rho, bc)

# I want to set "V_xc = (3/pi*V_rho)**(1/3)", is the following right?
V_xc=Function(V)
V_xc = (3/pi*V_rho)**(1/3)

# Is the setting of this variational problem right? Do we need to intepolate V_rho or V_xc?
c = 0.5*inner(nabla_grad(u), nabla_grad(v))*dx-u*v_ext*v*dx+u*V_rho*v*dx+u*V_xc*v*dx

All the questions are listed in the comments.

Thanks

Question information

Language:
English Edit question
Status:
Answered
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Houdong Hu (vincehouhou) said : 		#1

I test what I posted above sentence by sentence

everthing works fine till

# I want to set "V_xc = (3/pi*V_rho)**(1/3)", is the following right?
V_xc=Function(V)
V_xc = (3/pi*V_rho)**(1/3)

I can not do plot(V_xc)?

I want to try give the value to V_xc point by point. Thus I try to do "V_xc.vector()[:] = (3/pi*V_rho.vector())**(1/3)", but it seems vector doesnot not support ** or pow. How could I give the value to V_xc point by point?

Revision history for this message
Jan Blechta (blechta) said : 		#2

1. project or
2. interpolate or
3. set dof by dof

ad 3.: if function spaces of V_xc and V_rho are same then
    V_xc.vector()[:] = (3./pi*V_rho.vector()[:])**(1./3.)
should work as raising powers of numpy arrays is ok.

Revision history for this message
Houdong Hu (vincehouhou) said : 		#3

I tried
"V_xc.vector()[:] = (3./pi*V_rho.vector()[:])**(1./3.)"
and got the error
"TypeError: unsupported operand type(s) for ** or pow(): 'GenericVector' and 'float'"

Do you think we could do it directly like this
"V_xc = (3/pi*V_rho)**(1/3)"?

This could be run without errors, but I can not do plot(V_xc).

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

On 04/03/2013 07:51 AM, Houdong Hu wrote:
> Question #225321 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/225321
>
> Status: Answered => Open
>
> Houdong Hu is still having a problem:
> I tried
> "V_xc.vector()[:] = (3./pi*V_rho.vector()[:])**(1./3.)"

Try:

  V_xc.vector()[:] = 3./pi*V_rho.array()**(1./3.)

The error comes from

  V_rho.vector()[:]

just returning a copy of the GenericVector, which as the error message
says, does not support element wise power.

  V_rho.array()

on the other hand returns a NumPy array, which does support element wise
power and

  V_xc.vector()[:]

supports assignment with a NumPy array.

Johan

> and got the error
> "TypeError: unsupported operand type(s) for ** or pow(): 'GenericVector' and 'float'"
>
> Do you think we could do it directly like this
> "V_xc = (3/pi*V_rho)**(1/3)"?
>
> This could be run without errors, but I can not do plot(V_xc).
>

Revision history for this message
Jan Blechta (blechta) said : 		#5

On Wed, 03 Apr 2013 06:41:23 -0000
Johan Hake <email address hidden> wrote:
> Question #225321 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/225321
>
> Status: Open => Answered
>
> Johan Hake proposed the following answer:
> On 04/03/2013 07:51 AM, Houdong Hu wrote:
> > Question #225321 on DOLFIN changed:
> > https://answers.launchpad.net/dolfin/+question/225321
> >
> > Status: Answered => Open
> >
> > Houdong Hu is still having a problem:
> > I tried
> > "V_xc.vector()[:] = (3./pi*V_rho.vector()[:])**(1./3.)"
>
> Try:
>
> V_xc.vector()[:] = 3./pi*V_rho.array()**(1./3.)

Rather
   V_xc.vector()[:] = 3./pi*V_rho.vector().array()**(1./3.)

>
> The error comes from
>
> V_rho.vector()[:]
>
> just returning a copy of the GenericVector, which as the error message
> says, does not support element wise power.
>
> V_rho.array()

Rather
   V_rho.vector().array()

>
> on the other hand returns a NumPy array, which does support element
> wise power and
>
> V_xc.vector()[:]
>
> supports assignment with a NumPy array.
>
> Johan
>
> > and got the error
> > "TypeError: unsupported operand type(s) for ** or pow():
> > 'GenericVector' and 'float'"
> >
> > Do you think we could do it directly like this
> > "V_xc = (3/pi*V_rho)**(1/3)"?

This creates UFL expression which can be projected or used in any
variational form but DOLFIN has a priori no clue how to evaluate it
directly.

> >
> > This could be run without errors, but I can not do plot(V_xc).
> >
>

Can you help with this problem?

Provide an answer of your own, or ask Houdong Hu for more information if necessary.

To post a message you must log in.