Volume form not what I expected.

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 = x
values = x
values = sqrt(abs(1-x*x-x*x))
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:
Last query:
 Revision history for this message Jan Blechta (blechta) said on 2013-02-24: #1

That's probably discretization error:

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

 Revision history for this message Christopher Laing (9e9o1k-chris) said on 2013-02-24: #2

Excellent, cheers.