Laplace operator

There is no D command, maybe you mean Dx? You can also do this:

def laplace(f):
    return dot(grad(f))

Martin

On 23 May 2012 13:45, Sebastian Rau
<email address hidden> wrote:
> New question #198186 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/198186
>
> Hello,
>
> i would like to use the Laplace operator in a variational formulation, but it seems to that it's not predefined yet, right?
>
> And then i tried to define it by myself with the D( y, i) command, but my fenics version (how can I see which one i have) does not recognize this command either. Does anyone of you has a clue what's the problem?
>
> Thanks a lot!
>
> Regards
> Sebastian
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

Yes, i meant the Dx command, sorry!

But dot() needs two inputs!? So if nabla denotes the gradient operator, i should put something like

return dot(nabla, grad(f)) ?

Martin probably meant

  div(grad(f))

(div instead of dot)

As for your question on how to see which version you are using, just write

dolfin-version

in a terminal. Or in python:

import dolfin
print dolfin.__version__

Anders

Anders, thanks!

Marie and Martin, thanks, too, that seems to work somehow, but now i get a shape mismatch, i can't figure out yet where it comes from. If f is a function in let's say FunctionSpace(mesh, "CG", 1) (with the mesh being two dimensional), then div(grad(f)) is also in FunctionSpace(mesh, "CG", 0 or 1) ?

Please send a minimal running code example and the error message you get, then we might provide more help.

Ok, here is the shortest version i could extract:

import numpy
import numpy as np
from dolfin import *
from scitools.BoxField import *
import os
import shutil

# Create mesh and define function space
meshsize = 32
timesteps = 20
mesh = UnitSquare(meshsize,meshsize)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
Q = VectorFunctionSpace(mesh, "CG", 1)
W = MixedFunctionSpace([V,Q,V])

# Define Constants
dt = 10**(-4) # Time step
eps = 10**(-1)
lamda = 10**(-1)
Estart = 0

#######################################################

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS \
        or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

######################################################

# Define the Laplace operator
def lapl(f):
 return div(grad(f))

class State(object):
 def __init__(self):
  self.uk = Function(W)
  self.u0 = Function(W)
  self.E = Function(V)
  self.Vext = Function(V)
  self.problem = None

 def solve(self, u0, E, massn0):
  # self.E.interpolate(E)
  self.u0 = u0
  self.massn0 = massn0
  self.uk = Function(W, self.u0)
  self.du = TrialFunction(W)

  self.J = (4.0/dt)*(self.uk[0]*self.du[0] - self.du[0]*self.u0[0])*v[0]*dx \
   - 2.0*inner( div( dot(self.uk[1],self.uk[1])*self.du[0]/self.uk[0]**3), v[1] )*dx \
   - 2.0*( self.uk[0]*self.du[0]*(self.uk[2] + self.E*Vext[0]) + self.du[0]*lapl(self.uk[0]) + self.uk[0]*lapl(self.du[0]) )*div( v[1] )*dx \
   - 2.0*inner( grad(self.uk[0]*self.du[0])*(self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0])), v[1] )*dx \
   + eps**2*inner( grad(self.uk[0]**2)*(lapl(self.du[0])*self.uk[0] - lapl(self.uk[0])*self.du[0])/(self.uk[0]**2), v[1] )*dx \
   - 2.0*self.uk[0]*self.du[0]*dx \
   - inner(self.du[1], grad(v[0]))*dx + 1.0/dt*inner(self.uk[1],v[1])*dx + ( div( 2*dot(self.uk[1], self.du[1])/(self.uk[0]**2) ), v[1]) \
   - self.uk[0]**2*self.du[2]*div(v[1])*dx - inner( grad(self.uk[0]**2)*self.du[2],v[1] )*dx + lamda**2*inner( grad(self.du[2]),grad(v[2]) )*dx (***)

self.F = (2.0/dt)*(self.uk[0]**2 - self.uk[0]*self.u0[0])*v[0]*dx - inner(self.uk[1], grad(v[0]))*dx \
   + 1.0/dt*inner(self.uk[1] - self.u0[1], v[1])*dx + inner(div(1/(self.uk[0]**2)*dot(self.uk[1], self.uk[1])), v[1] )*dx \
   - self.uk[0]**2*( self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0]))*div(v[1])*dx \
   - inner( grad(self.uk[1])*(self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0])), v[1])*dx \
   + lamda**2*inner(grad(self.uk[1]), grad(self.uk[2]))*dx - self.uk[0]**2*v[2]*dx

  self.du = Function(W)
  self.err = 1.0
  self.i = 0
  while (self.err > 1e-5) and (self.i < maxiter):
   self.A = assemble(self.J)
   self.b = assemble(self.F)
   [bc.apply(self.A, self.b) for bc in bcfwd]
   solve(self.A, self.du.vector(), self.b)
   self.err = numpy.linalg.norm(self.du.vector().array(), ord=numpy.Inf)
   self.uk.vector()[:] = self.uk.vector() + 1.0*self.du.vector()
   self.i += 1
   print 'Error: ', self.err

  return self.uk

#######################################################

class StateT(object):
 def __init__(self):
  self.state = State()
  self.u = [Function(W)]*(timesteps + 1)

 def solve(self, uinit, E):
  self.u[0].vector()[:] = uinit.vector()
  self.massn0int = self.u[0][0]*self.u[0][0]*dx
  self.massn0 = assemble(self.massn0int, mesh = mesh)
  self.E = E
  for i in range(0,timesteps):
   self.u[i+1] = Function(self.state.solve(self.u[i], self.E[i+1], self.massn0))

  return self.u

#######################################################

E = [Function(R)]*(timesteps + 1)

# Initialize Control:
for i in range(0,timesteps+1):
 E[i] = Constant(Estart)

# Define Test Function for Forward Problem:
v = TestFunction(W)
# Define boundary condition for Forward Problem:
bvalue = Constant(0.0)
bvalue2 = (Constant(0.0), Constant(0.0))
bcfwd = [ DirichletBC(W.sub(0), bvalue, boundary), DirichletBC(W.sub(1), bvalue2, boundary), DirichletBC(W.sub(2), bvalue, boundary)]

Vext = Function(W)
V_ext = ExternalPotential()
Vext.interpolate(V_ext)

uinit = Function(W)
u_init = InitialConditions()
uinit.interpolate(u_init)
statet = StateT()
uopt = statet.solve(uinit, E)

And the following error message points to the line denoted with (***):
ufl.log.UFLException: Shape mismatch.

That is not short...

Johan

On 05/23/2012 05:45 PM, Sebastian Rau wrote:
> Question #198186 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/198186
>
> Status: Answered => Open
>
> Sebastian Rau is still having a problem:
> Ok, here is the shortest version i could extract:
>
> import numpy
> import numpy as np
> from dolfin import *
> from scitools.BoxField import *
> import os
> import shutil
>
> # Create mesh and define function space
> meshsize = 32
> timesteps = 20
> mesh = UnitSquare(meshsize,meshsize)
> V = FunctionSpace(mesh, "CG", 1)
> R = FunctionSpace(mesh, "R", 0)
> Q = VectorFunctionSpace(mesh, "CG", 1)
> W = MixedFunctionSpace([V,Q,V])
>
> # Define Constants
> dt = 10**(-4) # Time step
> eps = 10**(-1)
> lamda = 10**(-1)
> Estart = 0
>
> #######################################################
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> def boundary(x):
> return x[0]< DOLFIN_EPS or x[0]> 1.0 - DOLFIN_EPS \
> or x[1]< DOLFIN_EPS or x[1]> 1.0 - DOLFIN_EPS
>
> ######################################################
>
> # Define the Laplace operator
> def lapl(f):
> return div(grad(f))
>
> class State(object):
> def __init__(self):
> self.uk = Function(W)
> self.u0 = Function(W)
> self.E = Function(V)
> self.Vext = Function(V)
> self.problem = None
>
> def solve(self, u0, E, massn0):
> # self.E.interpolate(E)
> self.u0 = u0
> self.massn0 = massn0
> self.uk = Function(W, self.u0)
> self.du = TrialFunction(W)
>
> self.J = (4.0/dt)*(self.uk[0]*self.du[0] - self.du[0]*self.u0[0])*v[0]*dx \
> - 2.0*inner( div( dot(self.uk[1],self.uk[1])*self.du[0]/self.uk[0]**3), v[1] )*dx \
> - 2.0*( self.uk[0]*self.du[0]*(self.uk[2] + self.E*Vext[0]) + self.du[0]*lapl(self.uk[0]) + self.uk[0]*lapl(self.du[0]) )*div( v[1] )*dx \
> - 2.0*inner( grad(self.uk[0]*self.du[0])*(self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0])), v[1] )*dx \
> + eps**2*inner( grad(self.uk[0]**2)*(lapl(self.du[0])*self.uk[0] - lapl(self.uk[0])*self.du[0])/(self.uk[0]**2), v[1] )*dx \
> - 2.0*self.uk[0]*self.du[0]*dx \
> - inner(self.du[1], grad(v[0]))*dx + 1.0/dt*inner(self.uk[1],v[1])*dx + ( div( 2*dot(self.uk[1], self.du[1])/(self.uk[0]**2) ), v[1]) \
> - self.uk[0]**2*self.du[2]*div(v[1])*dx - inner( grad(self.uk[0]**2)*self.du[2],v[1] )*dx + lamda**2*inner( grad(self.du[2]),grad(v[2]) )*dx (***)
>
> self.F = (2.0/dt)*(self.uk[0]**2 - self.uk[0]*self.u0[0])*v[0]*dx - inner(self.uk[1], grad(v[0]))*dx \
> + 1.0/dt*inner(self.uk[1] - self.u0[1], v[1])*dx + inner(div(1/(self.uk[0]**2)*dot(self.uk[1], self.uk[1])), v[1] )*dx \
> - self.uk[0]**2*( self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0]))*div(v[1])*dx \
> - inner( grad(self.uk[1])*(self.uk[2] + self.E*Vext[0] - eps**2*lapl(self.uk[0])/(self.uk[0])), v[1])*dx \
> + lamda**2*inner(grad(self.uk[1]), grad(self.uk[2]))*dx - self.uk[0]**2*v[2]*dx
>
>
> self.du = Function(W)
> self.err = 1.0
> self.i = 0
> while (self.err> 1e-5) and (self.i< maxiter):
> self.A = assemble(self.J)
> self.b = assemble(self.F)
> [bc.apply(self.A, self.b) for bc in bcfwd]
> solve(self.A, self.du.vector(), self.b)
> self.err = numpy.linalg.norm(self.du.vector().array(), ord=numpy.Inf)
> self.uk.vector()[:] = self.uk.vector() + 1.0*self.du.vector()
> self.i += 1
> print 'Error: ', self.err
>
> return self.uk
>
> #######################################################
>
> class StateT(object):
> def __init__(self):
> self.state = State()
> self.u = [Function(W)]*(timesteps + 1)
>
> def solve(self, uinit, E):
> self.u[0].vector()[:] = uinit.vector()
> self.massn0int = self.u[0][0]*self.u[0][0]*dx
> self.massn0 = assemble(self.massn0int, mesh = mesh)
> self.E = E
> for i in range(0,timesteps):
> self.u[i+1] = Function(self.state.solve(self.u[i], self.E[i+1], self.massn0))
>
> return self.u
>
> #######################################################
>
> E = [Function(R)]*(timesteps + 1)
>
> # Initialize Control:
> for i in range(0,timesteps+1):
> E[i] = Constant(Estart)
>
> # Define Test Function for Forward Problem:
> v = TestFunction(W)
> # Define boundary condition for Forward Problem:
> bvalue = Constant(0.0)
> bvalue2 = (Constant(0.0), Constant(0.0))
> bcfwd = [ DirichletBC(W.sub(0), bvalue, boundary), DirichletBC(W.sub(1), bvalue2, boundary), DirichletBC(W.sub(2), bvalue, boundary)]
>
> Vext = Function(W)
> V_ext = ExternalPotential()
> Vext.interpolate(V_ext)
>
> uinit = Function(W)
> u_init = InitialConditions()
> uinit.interpolate(u_init)
> statet = StateT()
> uopt = statet.solve(uinit, E)
>
>
> And the following error message points to the line denoted with (***):
> ufl.log.UFLException: Shape mismatch.
>

Hi,

i think i figured out what the problem is, it's not related to the Laplace operator. So i opened a new Thread for it.

Thanks for your help!

Regards
Sebastian