Volume form not what I expected.

Asked by Christopher Laing on 2013-02-22

Hello everyone,

Let us suppose I want to create a finite element representation of half-sphere of radius one, where the surface is described by a 3-component vector phi = (X, Y, Z). I then wish to compute the volume of this half-sphere.

Here is a minimal example:

# Begin
from dolfin import *

mesh = UnitCircle(7) # I'm aware that this is deprecated, but I can't recall the new syntax

space= VectorFunctionSpace(mesh, "CG", 2, dim=3)

class phiexpression(Expression):
 def eval(self,values,x):
  values[0] = x[0]
  values[1] = x[1]
  values[2] = sqrt(abs(1-x[0]*x[0]-x[1]*x[1]))
 def value_shape(self):
  return (3,)

phi = Function(space)
phiexp = phiexpression()
phi.interpolate(phiexp)

V = (1./3)*dot(phi,cross(phi.dx(0),phi.dx(1))) # volume form

initial_guess_volume = assemble(V*dx)

print initial_guess_volume
# End

Now, this ought to give me 2 pi / 3, but instead I get 1.92919724808

What am I doing wrong? Is this just round-off error?

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Christopher Laing
Solved:
2013-02-24
Last query:
2013-02-24
Last reply:
2013-02-24
Jan Blechta (blechta) said : #1

That's probably discretization error:

mesh = UnitCircle(70) ... 2.07779064251
mesh = UnitCircle(700) ... 2.09275075086

Excellent, cheers.