### How to interpolate a initial condition for transient in Mixed formulation and static displacement

178

views

0

I have got o error to interpolate initial condition in my mixed formulation. I have a (0. , 0. ) displacement in my Elastic phasephield model. Could someone help me to understand how to interpolate initial condition for displacement ?

my code is:

from fenics import*

n = 20

mesh = UnitSquareMesh(n,n)

V = VectorElement("P", mesh.ufl_cell(), 2)

CG = FiniteElement("P", mesh.ufl_cell(), 1)

W = FunctionSpace(mesh, MixedElement([V, CG]))

#V0=FunctionSpace(mesh, "CG",1)

u, phi_u = TrialFunctions(W)

v, phi_v = TestFunctions(W)

#u_phi_u = TrialFunction(W)

#u, phi_u = split(u_phi_u)

#phi0_u=Constant(0.0)

#phi0_u = project(phi0_u, CG)

E = 1.E+0

nu = 0.3

b = Constant((0., 1.E-5))

s= Constant((1.E-5, 0.))

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

class inferior(SubDomain):

def inside(self,x,on_boundary):

tol = 1E-14

return abs(x[1]) < tol and on_boundary

class topo(SubDomain):

def inside(self,x,on_boundary):

tol = 1E-14

return abs(x[1]-1.) < tol and on_boundary

class esquerda(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[0]) < tol

class direita(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[0] - 1.) < tol

contorno = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

ds = Measure("ds", domain=mesh, subdomain_data=contorno)

contorno.set_all(0)

Direita = direita()

Direita.mark(contorno, 1)

Esquerda = esquerda()

Esquerda.mark(contorno, 2)

Topo = topo()

Inferior = inferior()

inferior_x = DirichletBC(W.sub(0),Constant((0.0,0.0)),Inferior)

superior_x = DirichletBC(W.sub(0),Constant((0.0,1.E-5)),Topo)

gf=1.0E-3

l=1.E-0

#BCS = [inferior_x, inferior_y, superior_x]

BCS = [inferior_x, superior_x]

#BCS = [inferior_x]

def Pii(phi):

return gf/l*phi-2.0*(1-phi)*(lmbda*tr(epsilon(u))**2+mu*(epsilon(u)**2))

def epsilon(v):

return 0.5*(grad(v) + grad(v).T)

def sigma(u):

return lmbda*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)

def factor(phi):

return ((1-phi)**2)

dt = 1.E-1

t = float(dt)

Tf=10*dt

beta=1.E-3

w = Function(W)

w0=Function(W)

u, phi_u = split(w)

u0, phi0_u = split(w0)

F = +inner(factor(phi_u)*sigma(u),grad(v))*dx \

- dot(b,v)*dx \

+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \

+ inner(gf*l*grad(phi_u),grad(phi_v))*dx \

+ dot(Pii(phi_u),phi_v)*dx \

+beta*(1.0/dt)*inner(phi_u-phi0_u, phi_v)*dx

R = action(F, w)

DR = derivative(R, w)

problem = NonlinearVariationalProblem(R, w,BCS, DR)

solver = NonlinearVariationalSolver(problem)

contador=0

while t <= Tf:

interactive()

my code is:

from fenics import*

n = 20

mesh = UnitSquareMesh(n,n)

V = VectorElement("P", mesh.ufl_cell(), 2)

CG = FiniteElement("P", mesh.ufl_cell(), 1)

W = FunctionSpace(mesh, MixedElement([V, CG]))

#V0=FunctionSpace(mesh, "CG",1)

u, phi_u = TrialFunctions(W)

v, phi_v = TestFunctions(W)

#u_phi_u = TrialFunction(W)

#u, phi_u = split(u_phi_u)

#phi0_u=Constant(0.0)

#phi0_u = project(phi0_u, CG)

E = 1.E+0

nu = 0.3

b = Constant((0., 1.E-5))

s= Constant((1.E-5, 0.))

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

class inferior(SubDomain):

def inside(self,x,on_boundary):

tol = 1E-14

return abs(x[1]) < tol and on_boundary

class topo(SubDomain):

def inside(self,x,on_boundary):

tol = 1E-14

return abs(x[1]-1.) < tol and on_boundary

class esquerda(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[0]) < tol

class direita(SubDomain):

def inside(self, x, on_boundary):

tol = 1E-14

return on_boundary and abs(x[0] - 1.) < tol

contorno = MeshFunction("size_t", mesh, mesh.topology().dim()-1)

ds = Measure("ds", domain=mesh, subdomain_data=contorno)

contorno.set_all(0)

Direita = direita()

Direita.mark(contorno, 1)

Esquerda = esquerda()

Esquerda.mark(contorno, 2)

Topo = topo()

Inferior = inferior()

inferior_x = DirichletBC(W.sub(0),Constant((0.0,0.0)),Inferior)

superior_x = DirichletBC(W.sub(0),Constant((0.0,1.E-5)),Topo)

gf=1.0E-3

l=1.E-0

#BCS = [inferior_x, inferior_y, superior_x]

BCS = [inferior_x, superior_x]

#BCS = [inferior_x]

def Pii(phi):

return gf/l*phi-2.0*(1-phi)*(lmbda*tr(epsilon(u))**2+mu*(epsilon(u)**2))

def epsilon(v):

return 0.5*(grad(v) + grad(v).T)

def sigma(u):

return lmbda*tr(epsilon(u))*Identity(2) + 2.0*mu*epsilon(u)

def factor(phi):

return ((1-phi)**2)

dt = 1.E-1

t = float(dt)

Tf=10*dt

beta=1.E-3

w = Function(W)

w0=Function(W)

u, phi_u = split(w)

u0, phi0_u = split(w0)

F = +inner(factor(phi_u)*sigma(u),grad(v))*dx \

- dot(b,v)*dx \

+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \

+ inner(gf*l*grad(phi_u),grad(phi_v))*dx \

+ dot(Pii(phi_u),phi_v)*dx \

+beta*(1.0/dt)*inner(phi_u-phi0_u, phi_v)*dx

R = action(F, w)

DR = derivative(R, w)

problem = NonlinearVariationalProblem(R, w,BCS, DR)

solver = NonlinearVariationalSolver(problem)

contador=0

while t <= Tf:

solver.solve()

# (u, phi_u) = w.split()

phi0_u.assign(phi_u)

t += float(dt)

contador+=1

plot(u, key="u", title='Solution at t = %g contador = %g' % (t, contador))

plot(phi_u, key="phi_u", title='Solution at t = %g contador = %g' % (t, contador))

interactive()

Community: FEniCS Project

### 2 Answers

2

Check out the paragraph

*in the FEniCS Tutorial, p. 80. (https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf)***Setting initial conditions for mixed systems**0

I solved it:

you should add w0.assign(w)

w = Function(W)

w0=Function(W)

u, phi_u = split(w)

u0, phi0_u = split(w0)

F = +inner(factor(phi_u)*sigma(u),grad(v))*dx \

- dot(b,v)*dx \

+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \

+ inner(gf*l*grad(phi_u),grad(phi_v))*dx \

+ dot(Pii(phi_u),phi_v)*dx \

+beta*(1.0/dt)*inner(phi_u-phi0_u, phi_v)*dx

I hope help someone.

Thanks all!!

you should add w0.assign(w)

w = Function(W)

w0=Function(W)

u, phi_u = split(w)

u0, phi0_u = split(w0)

F = +inner(factor(phi_u)*sigma(u),grad(v))*dx \

- dot(b,v)*dx \

+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \

+ inner(gf*l*grad(phi_u),grad(phi_v))*dx \

+ dot(Pii(phi_u),phi_v)*dx \

+beta*(1.0/dt)*inner(phi_u-phi0_u, phi_v)*dx

I hope help someone.

Thanks all!!

Plot(u, mode='displacement' ) or see demo_plot.py

written
5 weeks ago by
hsk

Please login to add an answer/comment or follow this question.

according to your code if you wants to separate vector element and scalar element, for Plotting

can you help, how to plot vector element separately, its a silly question, but it would help alot , i am a beginner

Thanks