Error in bilinear Form

[Akshay]

I was solving the Lame’s problem by exploiting symmetry. Domain is semi-circular pipe with some thickness. I have written the code but not able to figure out the error.

Code

from fenics import *
import numpy as np
import math
parameters['form_compiler']['cpp_optimize'] = True
parameters['form_compiler']['optimize'] = True
from mshr import*
set_log_active(False)
import matplotlib.pyplot as plt

mesh = Mesh("SE.xml")
plot(mesh)
plt.xlabel("X")
plt.ylabel("Y")
plt.show()


'''Printing Mesh parameters'''

print('calculating the Length Scale Parameter and Total no of Elements')
print ('Length Scale Parameter =',2*mesh.hmin())
print ('Total no of Elements = ', mesh.num_cells())

ndim = mesh.geometry().dim()

'''--------------------------------------------------------------------
#                       Defining Functionspaces
#--------------------------------------------------------------------'''

V = FunctionSpace(mesh,"CG",1)
W = VectorFunctionSpace(mesh,"CG",1)
WW = FunctionSpace(mesh,"DG",0)

p , q = TrialFunction(V) , TestFunction(V)
u , v = TrialFunction(W) , TestFunction(W)

'''--------------------------------------------------------------------
#                       Material Properties
#--------------------------------------------------------------------'''

E = 210e3
nu = 0.31
lmbda , mu = Constant(E*nu/((1.0 + nu )*(1.0-2.0*nu))) , Constant(E/(2*(1+nu)))

'''--------------------------------------------------------------------
#                       Constitutive Relations
#--------------------------------------------------------------------'''

def epsilon(u):
    return sym(grad(u))

def sigma(u):
    return 2*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(ndim)

'''Strain Energy Functional and History Field Variable'''

def psi(u):
    return 0.5*(lmbda + mu)*(0.5*(tr(epsilon(u)) + abs(tr(epsilon(u)))))**2 +             mu*inner(dev(epsilon(u)),dev(epsilon(u)))
            
'''--------------------------------------------------------------------
#                       Subdividing Boundaries
#--------------------------------------------------------------------'''

top = CompiledSubDomain("near(x[0],0) and on_boundary")
bottom = CompiledSubDomain("near(x[1],0) and on_boundary")
Inner = CompiledSubDomain("near(((x[0]*x[0]-x0*x0) + (x[1]*x[1] - x1*x1) - r*r),0,0.1)",x0=0,x1=0,r=1)


'''--------------------------------------------------------------------
#                       Boundaries Conditions
#--------------------------------------------------------------------'''

bcbot = DirichletBC(W.sub(1),Constant(0),bottom)
bctop = DirichletBC(W.sub(0),Constant(0),top)


bc_u = [bcbot,bctop]

'''--------------------------------------------------------------------
#                       Numbering the Boundaries
#--------------------------------------------------------------------'''

Boundaries = MeshFunction('size_t',mesh,ndim-1)
Boundaries.set_all(0)
top.mark(Boundaries,1)
bottom.mark(Boundaries,2)
Inner.mark(Boundaries,3)
ds = Measure('ds')(subdomain_data = Boundaries)
n = FacetNormal(mesh)


'''--------------------------------------------------------------------
#                       Variational Formulation
#--------------------------------------------------------------------'''

u = Function(W)
f = Constant((1,1))
E = inner(sigma(u),epsilon(v))*dx - inner(f,v)*ds(3)
a = lhs(E)
l = rhs(E)
solve(a==l,u,bc_u)


# displacement plot
plot(u)
plt.colorbar(plot(u))
plt.show()

print ('Simulation done with no error')    



assemble(psi(unew)*dx)

Domain is as indicated :

In[ ]:

[dokken]

Please format your code by encapsulating each cell with ```, such that it is more readable.

[dokken]

@Akshay

u = Function(W)

This is your problem, you are refining your TrialFunction before using it in the variational form.
As your are solving this with a==l, u has to be a trialfunction

[Akshay]

So what is solution??What should I do with u = Function(W)

[dokken]

Move it after the definition of the bilinear form:

u,v = Trialfunction(W), TestFunction(W)

E = inner(sigma(u),epsilon(v))*dx - inner(f,v)*ds(3)
a = lhs(E)
l = rhs(E)
uh =Function(W)
solve(a==l,uh, bc_u)

Now your solution is stored in uh

[Akshay]

Thanks dokken…Kudos
It worked