← Back to team overview

dolfin team mailing list archive

[Bug 531207] [NEW] Problem with DOLFIN_EPS

 

Public bug reported:

I am having problems defining boundary conditions using DOLFIN_EPS.
See example below, which works when stretch=1.0, but fails if stretch> 4.5
I suppose this could be something to do with machine precision, so would
be interesting to see if it is reproducable on different systems...


from dolfin import *

stretch = 4.5

mesh = UnitSquare(40,20)
for x in mesh.coordinates():
    x[0] *= stretch

Vs = VectorFunctionSpace(mesh, "CG", 1)

E  = 10.0
nu = 0.35

mu    = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)


# Define Dirichlet boundaries (x = 0 or x = stretch)
def boundary0(x):
    return x[0] < DOLFIN_EPS
  
def boundary1(x):
    return x[0] > stretch - DOLFIN_EPS


# Define boundary condition
u0 = Expression(("0.0","0.0"))
bc0 = DirichletBC(Vs, u0, boundary0)
u1 = Expression(("0.1","0.0"))
bc1 = DirichletBC(Vs, u1, boundary1)
bcs = [bc0,bc1]

# Define variational problem
v = TestFunction(Vs)
u = TrialFunction(Vs)

f = Constant((0.0,0.0))
a = inner(grad(v), sigma(u))*dx
L = inner(v, f)*dx

# Compute solution
problem = VariationalProblem(a, L, bcs)
u = problem.solve()

# Plot solution
plot(u, mode="displacement", interactive=True)

mesh.move(u)
plot(mesh,interactive=True)

** Affects: dolfin
     Importance: Undecided
         Status: New

-- 
Problem with DOLFIN_EPS
https://bugs.launchpad.net/bugs/531207
You received this bug notification because you are a member of DOLFIN
Team, which is subscribed to DOLFIN.

Status in DOLFIN: New

Bug description:
I am having problems defining boundary conditions using DOLFIN_EPS.
See example below, which works when stretch=1.0, but fails if stretch> 4.5
I suppose this could be something to do with machine precision, so would
be interesting to see if it is reproducable on different systems...


from dolfin import *

stretch = 4.5

mesh = UnitSquare(40,20)
for x in mesh.coordinates():
    x[0] *= stretch

Vs = VectorFunctionSpace(mesh, "CG", 1)

E  = 10.0
nu = 0.35

mu    = E / (2.0*(1.0 + nu))
lmbda = E*nu / ((1.0 + nu)*(1.0 - 2.0*nu))

def sigma(v):
    return 2.0*mu*sym(grad(v)) + lmbda*tr(sym(grad(v)))*Identity(v.cell().d)


# Define Dirichlet boundaries (x = 0 or x = stretch)
def boundary0(x):
    return x[0] < DOLFIN_EPS
  
def boundary1(x):
    return x[0] > stretch - DOLFIN_EPS


# Define boundary condition
u0 = Expression(("0.0","0.0"))
bc0 = DirichletBC(Vs, u0, boundary0)
u1 = Expression(("0.1","0.0"))
bc1 = DirichletBC(Vs, u1, boundary1)
bcs = [bc0,bc1]

# Define variational problem
v = TestFunction(Vs)
u = TrialFunction(Vs)

f = Constant((0.0,0.0))
a = inner(grad(v), sigma(u))*dx
L = inner(v, f)*dx

# Compute solution
problem = VariationalProblem(a, L, bcs)
u = problem.solve()

# Plot solution
plot(u, mode="displacement", interactive=True)

mesh.move(u)
plot(mesh,interactive=True)





Follow ups

References