Hello everyone. I have a small doubt, in coupling of mechanical and phase field equation, When I tried to view this stable structure in paraview, for checking the implementation of boundary condition, unfortunately it is showing deformation in negative direction. Also the applied traction is also not visible. I hereby attach my code and paraview outcome,
import os
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import numpy as np
import ufl
from basix.ufl import element
from dolfinx import default_real_type, log, io
from dolfinx.fem import Function, functionspace,assemble_scalar
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.io import XDMFFile
from dolfinx.nls.petsc import NewtonSolver
from ufl import ds, dx, grad, inner
from dolfinx import fem, default_scalar_type,mesh
#Creating mesh
msh = mesh.create_rectangle (comm=MPI.COMM_WORLD, points=((0.0, 0.0), (6.0,1.0)), n=(384,64), cell_type=mesh.CellType.quadrilateral)
#Create function space
V = ufl.VectorElement("Lagrange", msh.ufl_cell(),1) #For deformation map
P1 = ufl.FiniteElement("Lagrange", msh.ufl_cell(),1)
ME = fem.FunctionSpace(msh, ufl.MixedElement([V, P1, P1]))
#Creating Test functions and function
a, q, v= ufl.TestFunctions(ME)
u = Function(ME) # current solution
u0 = Function(ME) # solution from previous converged step
# Split mixed functions
d, p, mu= ufl.split(u)
d0, p0, mu0= ufl.split(u0)
# Interpolate initial condition for the concentration
const_val = 0.5
u.sub(1).interpolate(lambda x: np.full((x.shape[1],), const_val))
u.x.scatter_forward()
#Applying boundary condition to beam
def left_corner(x):
return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],0.0))
ME0, _ = ME.sub(0).collapse()
bottom_value1 = fem.Function(ME0)
bottom_value1.x.array[:]= 0
left_corner_vertex = mesh.locate_entities_boundary(msh,0,left_corner)
dofs = fem.locate_dofs_topological((ME.sub(0),ME0),0,left_corner_vertex)
bc_left_corner = fem.dirichletbc(bottom_value1,dofs,ME.sub(0))
def right_corner(x):
return np.logical_and(np.isclose(x[1],0.0),np.isclose(x[0],6.0))
ME0y, _ = ME.sub(0).sub(1).collapse()
bottom_value2 = fem.Function(ME0y)
bottom_value2.x.array[:]= 0
right_corner_vertex = mesh.locate_entities_boundary(msh,0,right_corner)
right_dof = fem.locate_dofs_topological((ME.sub(0).sub(1),ME0y),0,right_corner_vertex)
bc_right_corner = fem.dirichletbc(bottom_value2, right_dof, ME.sub(0).sub(1))
bcs = [bc_left_corner, bc_right_corner]
#Defining traction and body forces
Traction=fem.Constant(msh,default_scalar_type((0,0)))
B=fem.Constant(msh,default_scalar_type((0,0)))
#Spatial dimension
dim=len(d)
#Identity tensor
I=ufl.variable(ufl.Identity(dim))
#Deformation gradient
F=ufl.variable(I+grad(d))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Elasticity parameters
E =default_scalar_type(74000)
nu = default_scalar_type(0.33)
G = fem.Constant(msh, E / (2 * (1 + nu)))
K = fem.Constant(msh, E / (3 * (1 - 2 * nu)))
void_c = default_scalar_type(0.001)
#Defining parameters
p = ufl.variable(p)
epsilon = default_scalar_type(1/64) # interface thickness
areainteng = default_scalar_type(0.3) #Area specific interface energy
dt = 5.0e-06 # time step
k = default_scalar_type(0.05) #Dissipation term
vp = ufl.variable((3 * p**2) - (2 * p**3))
void_c = default_scalar_type(0.001)
#Bulk energy expression
Psi = (G/2)*(((J**2/3) * Ic)-3) + (K/2)*(((J**2)/2)-(1/2)- ufl.ln(J))
Tot_Psi = ((1-vp)*void_c+ vp)*Psi
P = ufl.diff(Tot_Psi, F)
Int_Psi = (6*(1/epsilon)* areainteng)*((p**2)*(1-p)**2)
# dintdp = ufl.diff (Int_Psi, p)
# dpsidp = ufl.diff(Tot_Psi, p)
# dfdc = dintdp - dpsidp
psi_tot = Int_Psi + Tot_Psi
dfdc = ufl.diff (psi_tot, p)
ds = ufl.Measure("ds", domain=msh)
dx = ufl.Measure("dx", domain=msh)
#theta =1
#mu_mid = (1.0 - theta) * mu0 + theta * mu
#Weak form equations
F0 = inner(P, grad(a)) * dx - inner(B, a) * dx - inner(Traction, a) * ds
F1 = inner(mu,q)*dx - inner(dfdc,q)*dx - (3*(1/epsilon)*areainteng*inner(grad(p),grad(q)))*dx
F2 = inner(p,v)*dx - inner(p0,v)*dx + (dt*k*inner(grad(mu),grad(v)))*dx
NF = F0+F1+F2
#Solving Nonlinear problem
problem=NonlinearProblem(NF,u,bcs)
#Solver configurations
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.max_it = 20 # Adjust the maximum number of iterations
solver.atol = 1e-8 # Adjust the absolute tolerance
solver.rtol = 1e-6 # Adjust the relative tolerance
# Output file
file = XDMFFile(MPI.COMM_WORLD, "bitte3.xdmf", "w")
file.write_mesh(msh)
# We can customize the linear solver used inside the NewtonSolver by
# modifying the PETSc options
ksp = solver.krylov_solver
opts = PETSc.Options() # type: ignore
option_prefix = ksp.getOptionsPrefix()
opts[f"{option_prefix}ksp_type"] = "preonly"
opts[f"{option_prefix}pc_type"] = "lu"
opts[f"{option_prefix}pc_factor_mat_solver_type"] = "mumps"
ksp.setFromOptions()
log.set_log_level(log.LogLevel.INFO)
# Step in time
t = 0.0
# Reduce run time if on test (CI) server
if "CI" in os.environ.keys() or "GITHUB_ACTIONS" in os.environ.keys():
T = 1 * dt
else:
T = 5 * dt
V0, dofs1 = ME.sub(0).collapse()
val2 = dolfinx.fem.Function(V0)
V1, dofs2 = ME.sub(1).collapse()
val3 = dolfinx.fem.Function(V1)
d = u.sub(0)
u0.x.array[:] = u.x.array
while (t < T):
t += dt
r = solver.solve(u)
print(f"Step {int(t/dt)}: num iterations: {r[0]}")
u0.x.array[:] = u.x.array
val2.x.array[:] = u.x.array[dofs1]
val2.x.scatter_forward()
val3.x.array[:] = u.x.array[dofs2]
val3.x.scatter_forward()
file.write_function(d, t)
file.close()
Could some one please confirm me, whether these type of deformation behaviour encountered while coupling or there is error in my coding. In my opinion, I wrote all the equations correct and made boundary conditions of the left bottom corner fixed and right bottom corner constrained only in y direction, traction being applied at top middle. Unfortunately the output of the deformed and undeformed structure with the displacement values were with negative side deformation and applied traction on top middle surface is not visible,
Could some one please advice me over this regard. Thanks for your time and considerations. It will be very much helpful.