hii I have this error:"expression needs an integer argument 'degree' if no 'element' is provided

Question

the following code is presented

from dolfin import *
from math import *

External load

class Traction(Expression):
def init(self, end):

  Expression.__init__(self)
  self.t = 0.0
  self.end = end

def eval(self, values, x) :
values[0] = 0.0
values[1] = 0.0
if x[0] > 1.0 - DOLFIN_EPS:
values[0] = self.t/self.end if self.t < self.end else 1.0

def value_shape(self):
return (2)

def update(u, u0, v0, a0, beta, gamma, dt):
# Get vectors (references)
u_vec, u0_vec = u.vector(), u0.vector()
v0_vec, a0_vec = v0.vector(), a0.vector()

 # Update acceleration and velocity
 a_vec = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/

(0.5dtdt) -(1.0-2.0beta)a0_vec )

 # v = dt * ((1-gamma)*a0 + gamma*a) + v0
 v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec

 # Update (t(n) <-- t(n+1))
 v0.vector()[:], a0.vector()[:] = v_vec, a_vec
 u0.vector()[:] = u.vector()

Load mesh and define function space

mesh = UnitSquareMesh(32, 32)

Define function space

V = VectorFunctionSpace(mesh, "Lagrange", 1)

Test and trial functions

u1, w = TrialFunction(V), TestFunction(V)

E, nu = 10.0, 0.3
mu, lmbda = E/(2.0(1.0 + nu)), Enu/((1.0 + nu)(1.0 - 2.0nu))

Mass density and viscous damping coefficient

rho, eta = 1.0, 0.2

Time stepping parameters

beta, gamma = 0.25, 0.5
dt = 0.1
t, T = 0.0, 20*dt

Fields from previous time step (displacement, velocity, acceleration)

u0, v0, a0 = Function(V), Function(V), Function(V)
h = Traction(T/4.0)

Velocity and acceleration at t_(n+1)

v1 = (gamma/(betadt))(u1 - u0) - (gamma/beta - 1.0)v0 - dt(gamma/

(2.0beta)- 1.0)a0
a1 = (1.0/(beta*dt**2))*(u1 - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0

Stress tensor

def sigma(u, v):
return 2.0musym(grad(u)) + (lmbdatr(grad(u)) +etatr(grad(v)))

*Identity(u.cell().d)

Governing equation

F = (rhodot(a1, w) + inner(sigma(u1, v1), sym(grad(w))))dx - dot(h,

w)*ds

Extract bilinear and linear forms

a, L = lhs(F), rhs(F)

Set up boundary condition at left end

zero = Constant((0.0, 0.0))
def left(x):
return x[0] < DOLFIN_EPS
bc = DirichletBC(V, zero, left)

Set up PDE, advance in time and solve

u = Function(V)
problem = LinearVariationalProblem(a, L, u, bcs=bc)
solver = LinearVariationalSolver(problem)

Save solution in VTK format

while t <= T:
t += dt
h.t = t
solver.solve()
vtkfile = File('displacement_test/solution.pvd')
vtkfile << u
update(u, u0, v0, a0, beta, gamma, dt)

Plz I need your help to solve this problem

Thanks

Answer 1

Because when you are creation a Traction instance, it inherits the Expression constructor, so it still has the same requirements. Expressions need a degree to figure how many points of quadrature to use, which must be given when you create your instance. You could do

Expression("x[0]", degree = 1)

or you could give an element, which gives the expression information about the geometry you are using

el = V.ufl_element()
Expression("x[0]", element = el)

This should be added to your Traction custom expression. Best regards.

Comments

So where should I add this? Thanks

---

actually, I added this but i got an error of an invalid syntax...

---

It is not quite clear where to add it, please format your code correctly (adding spaces befores code lines for auto formatting). Probably, as you are overriding the init method, the element (or degree) is not being set correctly. For this, override init correctly as shown in the docs (https://fenicsproject.org/olddocs/dolfin/1.5.0/python/programmers-reference/functions/expression/Expression.html). Also, traction is considered as a boundary term, so instead of your 'if' statement, you could use

g = Expression ("min ( t/end, 1) ", degree = 1)

and then integrate over surface with a boundary differential:

g*v*ds

---

I tried to enhance the code quality by adding spaces, actually I didn't wrote this code I found it in a fenics developpement thesis and I think that this code can be very useful to me so please I will be gratefull if some one helps

from dolfin import *

External load

class Traction(Expression):
def init(self, end):
Expression.init(self)
self.t = 0.0
self.end = end

def eval(self, values, x):
values[0] = 0.0
values[1] = 0.0
if x[0] > 1.0 - DOLFIN_EPS:
values[0] = self.t/self.end if self.t < self.end else 1.0

def value_shape(self):
return (2,)

def update(u, u0, v0, a0, beta, gamma, dt):

Get vectors (references)

u_vec, u0_vec = u.vector(), u0.vector()

v0_vec, a0_vec = v0.vector(), a0.vector()

Update acceleration and velocity

a_vec = (1.0/(2.0beta))( (u_vec - u0_vec - v0_vecdt)/(0.5dtdt) -(1.0-2.0beta)*a0_vec )

v = dt * ((1-gamma)a0 + gammaa) + v0

v_vec = dt((1.0-gamma)a0_vec + gamma*a_vec) + v0_vec

Update (t(n) <-- t(n+1))

v0.vector()[:], a0.vector()[:] = v_vec, a_vec

u0.vector()[:] = u.vector()

Load mesh and define function space

mesh = UnitSquareMesh(32, 32)

Define function space

V = VectorFunctionSpace(mesh, "Lagrange", 1)

Test and trial functions

u1, w = TrialFunction(V), TestFunction(V)
E, nu = 10.0, 0.3
mu, lmbda = E/(2.0(1.0 + nu)), Enu/((1.0 + nu)(1.0 - 2.0nu))

Mass density and viscous damping coefficient

rho, eta = 1.0, 0.2

Time stepping parameters

beta, gamma = 0.25, 0.5

dt = 0.1

t, T = 0.0, 20*dt

Fields from previous time step (displacement, velocity, acceleration)

u0, v0, a0 = Function(V), Function(V), Function(V)

h = Traction(T/4.0)

Velocity and acceleration at t_(n+1)

v1 = (gamma/(betadt))(u1 - u0) - (gamma/beta - 1.0)v0 - dt(gamma/(2.0beta)- 1.0)a0

a1 = (1.0/(beta*dt**2))*(u1 - u0 - dt*v0) - (1.0/(2.0*beta) - 1.0)*a0

Stress tensor

def sigma(u, v):
return 2.0musym(grad(u)) + (lmbdatr(grad(u)) +etatr(grad(v)))*Identity(u.cell().d)

Governing equation

F = (rhodot(a1, w) + inner(sigma(u1, v1), sym(grad(w))))dx - dot(h, w)*ds

Extract bilinear and linear forms

a, L = lhs(F), rhs(F)

Set up boundary condition at left end

zero = Constant((0.0, 0.0))
def left(x):
return x[0] < DOLFIN_EPS

bc = DirichletBC(V, zero, left)

Set up PDE, advance in time and solve

u = Function(V)

problem = LinearVariationalProblem(a, L, u, bcs=bc)

solver = LinearVariationalSolver(problem)

Save solution in VTK format

file = File("displacement.pvd")

while t <= T:
t += dt
h.t = t
solver.solve()
update(u, u0, v0, a0, beta, gamma, dt)
file << u