0

I'm new to FEniCS and optimization. Trying to run my code I receive the Error UFLException: Can't add expressions with different shapes. The Code:

from dolfin import *
from dolfin_adjoint import *
mesh = Mesh("beam.xml")         # beam x[0] 0 - 5, x[1] 0 - 1, x[2] 0 - 1
V = VectorFunctionSpace(mesh, "CG", 1)
y = Function(V, name="State")   # Trialfunction State
d = y.geometric_dimension()     # space dimension
u = Function(V, name="Control") # Trialfunction Control
v = TestFunction(V)             # Testfunction
x = SpatialCoordinate(mesh)

alpha = Constant(10e-6)
beta = Constant(0.2)

rho = Constant(7800) # Density of steel [kg/m³]
E = Constant(210000)  # Young's modulus steel [MN/m²]
nu = Constant(0.3)  # Poisson's ratio steel

lmbda = E*nu/(1+nu)/(1-2*nu) # Lame's constant
mu = E/(2*(1+nu)) # Shear modulus

G = Constant(pi**2) #gravity

F_z = Constant(rho*G) #weight force
F = Constant((0,0,-F_z))
g_z = Constant(2.9575e5)
g = Constant((0,0,g_z))

y_Desired = Constant((0,0,0.0))

top = AutoSubDomain(lambda x: near(x[1], 1))

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0) 

top.mark(boundaries, 1)
ds_N = ds(subdomain_data=boundaries) # use in ds_N(1) as in inner(g,v)*ds_N(1)

bottom = AutoSubDomain(lambda x: near(x[1], 0))

boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0) 

bottom.mark(boundaries, 1)
ds_0 = ds(subdomain_data=boundaries) # use in ds_0(1) as in beta(abs(u))*ds_O(1)


def bottom(x, on_boundary):
    return(on_boundary and near (x[1], 0.0))

bc_bottom = DirichletBC(V, u, bottom)


def left(x, on_boundary):
    return(on_boundary and near (x[0], 0.0))
bc_left = DirichletBC(V, Constant((0, 0, 0)), left)


def right(x, on_boundary):
    return(on_boundary and near (x[0], 5.0))
bc_right = DirichletBC(V, Constant((0, 0, 0)), right)

bcs = [bc_left, bc_right, bc_bottom]

def epsilon(y):
    return 0.5*(nabla_grad(y) + nabla_grad(y).T)


def sigma (y):
    return lmbda*div(y)*Identity(d) + 2*mu*epsilon(y)

S = inner(sigma(y),epsilon(v))*dx - inner((F+u),v)*dx - inner(g,v)*ds_N(1)

solve(S == 0, y, bcs)

J = assemble((0.5*inner(y- y_Desired, y - y_Desired))*dx + alpha / 2 * inner((F+u), v)* dx + beta*abs(u)*ds_0)

I tried to make y_Desired also 3D but then I receive UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (3,) and free indices with labels ().

Can anyone help me on this?

DocGoccia
  • 1
  • 1

0 Answers0