Compile error of a expression in the weak from


221
views
2
4 months ago by
Hi all,

 

I tried to implement the 2D case of the problem in Fenics but got compile errors about the term  $J=tr\left(\nabla u\right)$J=tr(u) in B(u,X). There is no compile error if I remove the term  $J=tr\left(\nabla u\right)$J=tr(u) in B(u,X) (although it does not converge if I remove it). The code is as follows.

from dolfin import *
import os, shutil, math, sys, sympy
import numpy as np
  
# Create mesh and define function space
mesh = UnitSquareMesh(10, 10)

d = mesh.topology().dim()
e1 = [Constant([1.,0.]),Constant((1.,0.,0.))][d-2]
e2 = [Constant([0.,1.]),Constant((0.,1.,0.))][d-2]

# Mark boundary subdomians
class top(SubDomain):
	def inside(self,x,on_boundary):
		return near(x[1], 1) # change the number 

top = top()
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
top.mark(boundaries, 1)
ds = Measure("ds")[boundaries]

# Define functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
du = TrialFunction(V)            # Incremental displacement
v  = TestFunction(V)             # Test function
u  = Function(V)                 # Displacement from previous iteration

# Define Dirichlet boundary 
utop = Expression(("0", "-0.01*t"), t = 0.0, degree = 1)
bcs = DirichletBC(V, utop, top)

# Elasticity parameters
#E, nu = 10.0, 0.3
#mu, lmbda = E/(2*(1 + nu)), E*nu/((1 + nu)*(1 - 2*nu))
mu, lmbda = 4.4e5, 4.4e5

# Strain and stress
def eps(v):
    return sym(grad(v))

def sigma(v):
    return 2.0*mu*(eps(v)) + lmbda*tr(eps(v))*Identity(d)

x = SpatialCoordinate(mesh)

#define B(u,x)
def B(u,x):
	return Expression(("0", "x[1]+u1<d0 ? J/(x[1]+u1): J**2/(x[1]+u1)"), \
			u1 = u.sub(1), J = tr(eps(u)), d0 = 0.5, degree = 1)


T = Constant((0.0, 0.0))  # Traction force on the boundary1

# Stored strain energy density
psi = 0.5*inner(sigma(u),eps(u))

# Total potential energy
Pi = psi*dx - dot(B(u,x), u)*dx - dot(T, u)*ds

# Compute first variation of Pi (directional derivative about u in the direction of v)
Pi_u = derivative(Pi, u, v)

# Compute Jacobian of F
Pi_u_u = derivative(Pi_u, u, du)

# Variational problem and solver
problem_u = NonlinearVariationalProblem(Pi_u, u, bcs, Pi_u_u)
solver_u = NonlinearVariationalSolver(problem_u)

# Solver parameters
prm = solver_u.parameters
info(prm, True)

prm['nonlinear_solver'] = 'snes'
prm['snes_solver']['line_search'] = 'basic'
prm['snes_solver']['linear_solver']= 'cg'
prm['snes_solver']['absolute_tolerance'] = 1e-6
prm['snes_solver']['relative_tolerance'] = 1e-6
prm['snes_solver']['maximum_iterations'] = 200


# load steps
load_min = 0. # load multiplier min value
load_max = 1. # load multiplier max value
load_steps = 201 # number of time steps
load_multipliers = np.linspace(load_min,load_max,load_steps)
forces = np.zeros((len(load_multipliers),2))
file_u = File(savedir+"/u.pvd");
file_s = File(savedir+"/stress.pvd");

for (i_t,t) in enumerate(load_multipliers):
	utop.t = t
	solver_u.solve()
	
	W = TensorFunctionSpace(mesh, "Lagrange", 2)
	sigma_w = project(2*mu*sym(grad(u)) + lmbda*tr(grad(u))*Identity(d), W)

	file_u << u;
	file_s << sigma_w
	
	forces[i_t] = np.array([t*0.01,assemble(inner(sigma(u)*e2,e2)*ds(1))])
	np.savetxt(savedir+"/force.txt",forces)

​

I got the error:
In instant.recompile: The module did not compile with command 'make VERBOSE=1', see '/home/fenics/.cache/instant/python2.7/error/dolfin_18092345ac056e1dd526b1f67f9b1976ee5ad336/compile.log'
Traceback (most recent call last):
File "bftest.py", line 84, in <module>
Pi = psi*dx - dot(B(u,x), u)*dx - dot(T, u)*ds
File "bftest.py", line 76, in B
u1 = u.sub(1), J = tr(eps(u)), degree = 1)
File "/usr/local/lib/python2.7/dist-packages/dolfin/functions/expression.py", line 679, in __new__
mpi_comm=kwargs.get("mpi_comm"))
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/expressions.py", line 266, in compile_expressions
mpi_comm=mpi_comm)
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/expressions.py", line 183, in compile_expression_code
mpi_comm=mpi_comm)
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 70, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/lib/python2.7/dist-packages/dolfin/compilemodules/compilemodule.py", line 603, in compile_extension_module
**instant_kwargs)
File "/usr/local/lib/python2.7/dist-packages/instant/build.py", line 577, in build_module
build_system)
File "/usr/local/lib/python2.7/dist-packages/instant/build.py", line 170, in recompile
instant_error(msg % (cmd, compile_log_filename_dest))
File "/usr/local/lib/python2.7/dist-packages/instant/output.py", line 96, in instant_error
raise RuntimeError(text)
RuntimeError: In instant.recompile: The module did not compile with command 'make VERBOSE=1', see '/home/fenics/.cache/instant/python2.7/error/dolfin_18092345ac056e1dd526b1f67f9b1976ee5ad336/compile.log'

I appreciate any help on how to resolve this issue.

Thank you in advance.
Community: FEniCS Project

1 Answer


4
4 months ago by
Stripping down the definition of B(u, x) to
def B(u, x): return Expression(('0', 'J'), J=tr(eps(u)), degree=1)​

gives a more lucid error message:

TypeError: expected default arguments for member variables to be scalars or a scalar GenericFunctions.


I was able to run your code with
d0 = 0.5
def B(u, x):
        u1 = u.sub(1)
        J  = tr(eps(u))
        return as_vector([0.0, conditional(lt(x[1]+u1, d0), J/(x[1]+u1), J**2/(x[1]+u1))])​

but you must use the uflacs representation

parameters['form_compiler']['representation'] = 'uflacs'

see https://bitbucket.org/fenics-project/ffc/issues/107/ffc-fails-with-conditionals-which-are

Thank you, Adam. It helps!
written 4 months ago by Wilhelm  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »