Thank you Dokken for the solution. It works perfectly fine now. I have attached the code at the end for others. Although I have 2 questions.
- I see that I did not understand this correctly. Could you shortly explain (or point me to the right source/documentation) why it has to be like you mentioned? I usually fix, eg: displacement in x direction in the following way. Therefore, I tried it in a similar way. Why does the following way did not work.
# bottom surface (facet tag 3) of cylinder in x direction
facets = facet_tag.indices[facet_tag.values == 3] # fixed in x
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))
- I want to calculate stress/force on the top surface of the cylinder so that I can analyse the stress strain relations in the end. The idea was to calculate the traction on the top surface and multiply it by the area of the circular surface. For this, I added the code mentioned below in the time loop and it calculates it. However I am facing issues to save this in a file. I usually use fem. Expression() to evaluate and save data using VTXWriter. But I see that Expressions won’t work for integral measures. What would be the right way to save this data in an external file?
#---------stress on top surface-------#
traction = ufl.inner(df, ufl.dot(P,normal_vec))
Top_force = fem.assemble_scalar(fem.form(traction * ds(2)))
Top_stress = Top_force / (3.14 * pow(0.105,2))
print('traction =',Top_stress)
Here is the code that works fine now. Instead of NeoHookean, I have used a model called Yeoh to represent elasticity as I have the material parameters available for the compression test I am trying to simulate. It is a poroelastic problem.
from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl
import dolfinx
print(f"DOLFINx version: {dolfinx.__version__} based on GIT commit: {dolfinx.git_commit_hash} of https://github.com/FEniCS/dolfinx/")
from mpi4py import MPI
from dolfinx import fem
from petsc4py.PETSc import ScalarType
from petsc4py import PETSc #solvers
from pathlib import Path
from dolfinx.io import gmshio
import meshio
#------------------------------meshio-------------------------------#
folder = "yeoh_bi"
#---------------------- 2nd order tri elements (set order 1 in gmsh)-------------------#
proc = MPI.COMM_WORLD.rank
# Extract data and returns meshio mesh including physical markers
def create_mesh(mesh, cell_type, prune_z=False):
cells = mesh.get_cells_type(cell_type)
cell_data = mesh.get_cell_data("gmsh:physical", cell_type)
points = mesh.points[:, :2] if prune_z else mesh.points
out_mesh = meshio.Mesh(points=points, cells={cell_type: cells}, cell_data={"name_to_read": [cell_data.astype(np.int32)]})
return out_mesh
if proc == 0:
# Read in mesh
msh = meshio.read(Path("geometry")/"cylinder_21.msh") #"meshmixer_freecad_expand_v3.msh") #"meshmixer_freecad_v3.msh")
print(msh)
# Create and save one file for the mesh, and one file for the facets
triangle_mesh = create_mesh(msh, "tetra")
line_mesh = create_mesh(msh, "triangle")
meshio.write(Path(folder)/"mesh.xdmf", triangle_mesh)
meshio.write(Path(folder)/"mt.xdmf", line_mesh)
MPI.COMM_WORLD.barrier()
#Read the data in DOLFINx using XDMFFile.read_mesh and XDMFFile.read_meshtags.
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, Path(folder)/"mesh.xdmf", "r") as xdmf:
domain = xdmf.read_mesh(name="Grid") # cell tags
ct = xdmf.read_meshtags(domain, name="Grid")
domain.topology.create_connectivity(domain.topology.dim, domain.topology.dim - 1)
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, Path(folder)/"mt.xdmf", "r") as xdmf:
facet_tag = xdmf.read_meshtags(domain, name="Grid") # facet tags
fdim = domain.topology.dim - 1
dim = domain.topology.dim
#-----------------------Function Spaces and Mixed Space for displacement and pressure-------------------#
u_el = ufl.VectorElement("CG", domain.ufl_cell(), 2) # u quad element
p_el = ufl.FiniteElement("CG", domain.ufl_cell(), 1) # p linear element
element = ufl.MixedElement([u_el, p_el]) # Mixed Element
V = fem.FunctionSpace(domain, element) # Mixed Function space for all dofs
up = fem.Function(V) # total dofs in the problem/ func for all dofs
u_sp = fem.FunctionSpace(domain, u_el) # u Function space
p_sp = fem.FunctionSpace(domain, p_el) # p Function space
Vu, up_to_u = V.sub(0).collapse() # dofs in u func subspace
u_n = fem.Function(Vu) # dofs in u func
Vp, up_to_p = V.sub(1).collapse() # dofs in u func subspace
p_n = fem.Function(Vp) # dofs in u func
u, p = ufl.split(up)
print("all dofs:", len(u))
(del_u, del_p) = ufl.TestFunctions(V)
#---------------------------temporal parameters---------------------------#
# Define temporal parameters
t = 0 # Start time
T = 100 #80000
num_steps = 100
dt = T / num_steps # time step size
#---------------------------Initial Conditions---------------------------#
# Initialize history values
up.x.array[:] = 0.0
u_n.x.array[:] = 0.0
#---------------------------Boundary Conditions---------------------------#
def bc_vec_zero(x):
return np.zeros((dim, x.shape[1]), dtype=PETSc.ScalarType)
def bc_scalar_zero(x):
return np.zeros((1, x.shape[1]), dtype=PETSc.ScalarType)
list_bc = []
# Displacement Boundary
u_bc = fem.Function(u_sp)
u_bc.interpolate(bc_vec_zero)
#-----------------Fixed Bottom (3)---------------------
facets = facet_tag.indices[facet_tag.values == 3] # fixed in x
left_dofs = fem.locate_dofs_topological((V.sub(0).sub(0), u_sp.sub(0)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, left_dofs, V.sub(0)))
facets = facet_tag.indices[facet_tag.values == 3] # fixed in y
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(1), u_sp.sub(1)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))
facets = facet_tag.indices[facet_tag.values == 3] # fixed in z
bottom_dofs = fem.locate_dofs_topological((V.sub(0).sub(2), u_sp.sub(2)), fdim, facets)
list_bc.append(fem.dirichletbc(u_bc, bottom_dofs, V.sub(0)))
#--------- Compressive Strain on Top Surface (2)------------#
displacement_rate = 0.001 # 1mm/s per second
def axial_strain(t):
if t <= 3.0:
return -t * displacement_rate
else:
return -0.003
def u_exact(x):
return np.arry([0, 0, -0.001], (dim, x.shape[1]), dtype=PETSc.ScalarType)
u_bc_2 = fem. Constant(domain, ScalarType(axial_strain(t)))
top_facets = facet_tag.indices[facet_tag.values == 2]
top_dofs = fem.locate_dofs_topological(V.sub(0).sub(2), fdim, top_facets)
list_bc.append(fem.dirichletbc(u_bc_2, top_dofs, V.sub(0).sub(2)))
#----Normal Vector--------#
normal_vec = ufl.FacetNormal(domain)
#------------------------------Parameters-----------------------------------#
# Parameters for hyperelastic Yeoh model (could be neohookean as well)
c1 = 0.29 #material parameters
c2 = 1.01 #material parameter
d = 1.60 #material parameter
e= 0.90 #void ratio
k = 1.1e-10 #permeability
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
F = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)
b = ufl.variable(F * F.T)
# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
# Time Derivatives
vel_s = (u - u_n) / dt # Solid Velocity
# Gradients
vel_grad = 0.5 * (ufl.grad(vel_s) + ufl.transpose(ufl.grad(vel_s)))
grad_p = ufl.grad(p)
#Void ratio
nF = e/(1+e)
#--------------- Darcy Law --------------#
dar = k * grad_p
w_FS = dar/nF
#-----Yeoh stress (could be Neo Hookean as well)--------#
P_yeoh_b = ((c1) + (c2 * 2 * (Ic - 3))) * (1/J)
P_yeoh_i = ((1/d) * (J - 1))
P_yeoh = (P_yeoh_b * b) + (P_yeoh_i * I)
#------Total Stress----#
P = P_yeoh - p*I
# Integration measures
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)
#-----Weak Forms--------#
#Momentum Balance
Mom_B = ufl.inner(P, ufl.grad(del_u)) * dx # Solid Stress part
#Volume Balance
Vol_M = (ufl.div(vel_s) * del_p) * dx # div velocity
Vol_M += - ufl.inner(dar, ufl.grad(del_p)) * dx # darcy
weak_form = Mom_B + Vol_M
#-------------------------------Solver----------------------------#
# Initialise non-linear problem
problem = NonlinearProblem(weak_form, up, list_bc)
# Initialise newton solver
solver = NewtonSolver(domain.comm, problem)
# Set Newton solver options
solver.atol = 1e-8 #1e-8
solver.rtol = 1e-8 #1e-8
solver.convergence_criterion = "incremental"
solver.max_it = 50
solver.report = True
# Configure mumps
ksp = solver.krylov_solver
opts = PETSc.Options()
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()
#-----------------------------Solution Settings-----------------------------#
# Create solution variables
Vu_sol, up_to_u_sol = V.sub(0).collapse() # dofs in u func subspace
u_sol = fem.Function(Vu_sol)
Vp_sol, up_to_p_sol = V.sub(1).collapse() # dofs in u func subspace
p_sol = fem.Function(Vp_sol)
V_vec = fem.VectorFunctionSpace(domain,("CG", 2)) #
grad_p_o = fem.Function(V_vec) #
vel_s_o = fem.Function(u_sp)
V_tensor = fem.FunctionSpace(domain,("CG", 2, (3,3))) # 2 is for 2nd order elements. (3,3) is size of stress tensor
stress_o = fem.Function(V_tensor)
vel_grad_o = fem.Function(V_tensor)
# Force top surface space and function
force_FS = fem.VectorFunctionSpace(domain,("CG", 1)) #fem.VectorFunctionSpace(domain, "Real", 0)
df = ufl.TestFunction(force_FS)
# Expressions for output
w_FS_out_exp = fem.Expression(w_FS, u_sp.element.interpolation_points())
grad_p_out_exp = fem.Expression(grad_p, u_sp.element.interpolation_points()) #
stress_out_exp = fem.Expression(P, u_sp.element.interpolation_points())
vel_s_out_exp = fem.Expression(vel_s, u_sp.element.interpolation_points())
# Rename for Paraview file
u_sol.name = "displacement"
p_sol.name = "pressure"
grad_p_o.name = "grad_p"
vel_s_o.name = "vel_s"
stress_o.name = "stress"
vel_grad_o.name = "vel_grad"
# Writing data to files
file1 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"disp.bp", u_sol, "bp4")
file2 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"press.bp", p_sol, "bp4")
file3 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"grad_p.bp", grad_p_o, "bp4") #
file4 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"vel_s.bp", vel_s_o, "bp4")
file5 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"stress.bp", stress_o, "bp4") #
file6 = dolfinx.io.VTXWriter(domain.comm, Path(folder)/"vel_grad.bp", vel_grad_o, "bp4") #
#---------------------------Solve problem-------------------------------------#
duration_solve = 0.0
for n in range(num_steps):
# Update load/displacement
u_bc_2.value = axial_strain(t)
print('disp =',u_bc_2.value)
# Solve newton-steps
log.set_log_level(log.LogLevel.INFO)
duration_solve -= MPI.Wtime()
niter, converged = solver.solve(up)
duration_solve += MPI.Wtime()
PETSc.Sys.Print(
"Phys. Time {:.4f}, Calc. Time {:.4f}, Num. Iter. {}".format(
t, duration_solve, niter
)
)
u_n.x.array[:] = up.x.array[up_to_u]
u_sl = up.sub(0).collapse()
p_sl = up.sub(1).collapse()
u_sl.x.scatter_forward()
p_sl.x.scatter_forward()
u_sol.interpolate(u_sl)
p_sol.interpolate(p_sl)
grad_p_o.interpolate(grad_p_out_exp) #
vel_s_o.interpolate(vel_s_out_exp)
stress_o.interpolate(stress_out_exp) #
#---------stress on top surface-------#
traction = ufl.inner(df, ufl.dot(P,normal_vec))
Top_force = fem.assemble_scalar(fem.form(traction * ds(2)))
Top_stress = Top_force / (3.14 * pow(0.105,2))
print('traction =',Top_stress)
file1.write(t)
file2.write(t)
file3.write(t)
file4.write(t)
file5.write(t)
file6.write(t) #
#Update time
t += dt
file1.close()
file2.close()
file3.close()
file4.close()
file5.close()
file6.close() #
The msh file usedcan be found on this link: cylinder_21.msh - Google Drive
