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

277
views
0
8 months ago by
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):

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)

- dot(b,v)*dx \
+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \
+ 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)

while t <= Tf:
solver.solve()
# (u, phi_u) = w.split()
phi0_u.assign(phi_u)
t += float(dt)
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
8 months ago by
Check out the paragraph Setting initial conditions for mixed systems in the FEniCS Tutorial, p. 80. (https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf)
0
8 months ago by
I solved it:

w = Function(W)
w0=Function(W)

u, phi_u = split(w)
u0, phi0_u = split(w0)

- dot(b,v)*dx \
+ dot(s,v)*ds(1) - dot(s,v)*ds(2) \
+ 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!!

Hello, your code Helped me, to solve my  issue ,
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
written 4 months ago by Md Ayub Ansari
Plot(u, mode='displacement' ) or see demo_plot.py
written 4 months ago by hirshikesh