Dear FEniCS community
I’m a new user of the FEniCS library. I would really appreciate your help.
I’m trying to implement the variational formulation reported in this article:A structural model for the viscoelastic behavior of arterial walls: Continuum formulation and finite element analysis - ScienceDirect
The variational reported in the article is for a hyperelastic fiber-reinforced material with viscoelasticity
In this first version of the code, I only try to solve the hyperelastic problem.
I looked for a solution in the several topics open in the blog, which helped me a lot in setting the boundary condition, but the problem with the variational formulation still remains.
Here I reported my code:
from dolfinx import *
from matplotlib import pyplot as plt
import meshio
from dolfinx.io import gmshio
import pyvista
import numpy as np
from mpi4py import MPI
import math
import os
import sys
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import pyvista
from dolfinx import mesh, fem, plot, default_scalar_type, log
from basix.ufl import element, mixed_element
import ufl
#Create the mesh
length = 10
width = 4
depth = 2
num_ele_along_depth = 4
size_ele = depth/num_ele_along_depth
print(f"Dimensione elementi : {size_ele}")
n = [int(length/size_ele), int(width/size_ele), int(depth/size_ele)]
print(f"Numero elementi lungo x: {n[0]} \nNumero elementi lungo y: {n[1]} \nNumero elementi pungo z : {n[2]}")
domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([length, width, depth])], n, cell_type=mesh.CellType.hexahedron)
print("Ho creato il mio dominio trave 3D")
print(domain)
print(f"Numero di vertici: {domain.topology.index_map(0).size_local}")
print(f"Numero di elementi: {domain.topology.index_map(domain.topology.dim).size_local}")
print(f"Dimensione topologica: {domain.topology.dim}") # Dovrebbe stampare 3
print(f"Dimensione geometrica: {domain.geometry.dim}") # Dovrebbe stampare 3
#Visualization of domain
import pyvista as pv
from dolfinx import plot
topology, cell_types, geometry = plot.vtk_mesh(domain)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)
plotter = pv.Plotter()
plotter.add_mesh(grid, show_edges=True)
plotter.add_text("Definizione della mesh della trave 3D", position="upper_edge", font_size=12, color="blue")
plotter.view_isometric()
plotter.show_axes()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("Image_domain_3D.png")
#Create mixed element
Ue = element("CG", domain.ufl_cell().cellname(), 2, shape=(3,)) #Displacement
Pe = element("CG", domain.ufl_cell().cellname(), 1) #Pressure
Je = element("CG", domain.ufl_cell().cellname(), 1) #Volumetric strain
Me = mixed_element([Ue, Pe, Je])
print("Ho creato gli elementi mixati")
#Create function space
U = fem.functionspace(domain, Ue)
P = fem.functionspace(domain, Pe)
J = fem.functionspace(domain, Je)
print(f"Function spaces created: {U}, {P}, {J}")
M = fem.functionspace(domain, Me)
print(f"Function space created: {M}")
#Definizione delle condizioni a contorno
tdim = domain.topology.dim
fdim = domain.topology.dim - 1
def boundary_fixed(x):
return np.isclose(x[0], 0)
fixed_facets = mesh.locate_entities_boundary(domain, fdim, boundary_fixed)
marked_facets = np.hstack(fixed_facets)
marked_values = np.hstack(np.full_like(fixed_facets, 1))
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])
u_fixed = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
fixed_dofs = fem.locate_dofs_topological(U, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_fixed, fixed_dofs, U)]
domain.topology.create_connectivity(fdim, 0) # Ensure facet-to-node connectivity
facet_to_vertex = domain.topology.connectivity(fdim, 0)
# Extract all unique node indices from marked facets
fixed_vertices = np.unique(np.hstack([facet_to_vertex.links(facet) for facet in facet_tag.find(1)]))
#Visualization of boundary conditions
import pyvista as pv
# Convert mesh and boundary facets to PyVista format plot.vtk_mesh(domain)
topology, cell_types, geometry = plot.vtk_mesh(domain, domain.topology.dim)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)
# Find boundary facets
bndry_topology, bndry_cell_types, bndry_geometry = plot.vtk_mesh(domain, fdim)
boundary_grid = pv.UnstructuredGrid(bndry_topology, bndry_cell_types, bndry_geometry)
# Extract only fixed boundary facets
fixed_coordinates = domain.geometry.x[fixed_vertices]
#Per verificare che anche quelli fissi siano corretti
dof_coordinates = U.tabulate_dof_coordinates()[fixed_dofs]
# Plot mesh and highlight fixed nodes
plotter = pv.Plotter()
plotter.add_mesh(grid, style="wireframe", color="white") # Plot full mesh
plotter.add_mesh(fixed_coordinates, color="red", render_points_as_spheres=True, point_size=10)
# Replace fixed_coordinates with dof_coordinates in visualization
plotter.add_mesh(dof_coordinates, color="green", render_points_as_spheres=True, point_size=10)
plotter.add_text("Dirichlet DOFs (Green)", position="upper_edge", font_size=12, color="green")
plotter.add_text("Boundary Conditions (Fixed Nodes in Red)", position="upper_edge", font_size=12, color="blue")
print("Visualizzazione pronta")
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("Image_Boundary_Condition.png")
#Traction vector
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
#Body force vector
B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
##################################SECOND VISUALIZATION
###############################################
###########VERSIONE 1
# # Define functions
# u = fem.Function(U)
# p = fem.Function(P)
# j = fem.Function(J)
# m = fem.Function(M)
# m.x.array[:] = 0.0 #This will help Newton’s solver start from a reasonable initial state.
# print("Functions defined")
# δu, δp, δj = δw = ufl.TestFunctions(M)
#####################################################################FINE
#############VERSIONE 2
m = fem.Function(M)
p, j = ufl.split(m)[1], ufl.split(m)[2],
u = ufl.split(m)[0]
m.x.array[:] = 0.0
δu, δp, δj = δw = ufl.TestFunctions(M)
#####################################################################FINE
assert M.num_sub_spaces == 3, f"Error: Expected 3 subspaces in M, but got {M.num_sub_spaces}"
print(f"Number of subspaces in M: {M.num_sub_spaces}")
print(f"Displacement space shape: {M.sub(0).element.space_dimension}")
print(f"Pressure space shape: {M.sub(1).element.space_dimension}")
print(f"Volumetric strain space shape: {M.sub(2).element.space_dimension}")
####### VARIATIONAL FORMULATION #######
#######################################
# Spatial dimension
d = len(u)
# Identity tensor
I = ufl.variable(ufl.Identity(d))
# Deformation gradient
FF = ufl.variable(I + ufl.grad(u))
# Right Cauchy-Green tensor
C = ufl.variable(ufl.transpose(FF) * FF)
# Invariants of deformation tensors
Ic1 = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(FF))
# distortional part of C
C_dist = ufl.variable(J**(-2/3)*C)
Ic1_dist = ufl.variable(ufl.tr(C_dist))
#Fiber definition and orientation
a0 = fem.Constant(domain, default_scalar_type((1, 0, 0)))
Ic4 = ufl.variable(ufl.inner(a0, C*a0))
Ic4_dist = ufl.variable(J**(-2/3)*Ic4)
A0=ufl.variable(ufl.outer(a0,a0))
# Material parameters
C_m=default_scalar_type(100)
K1=default_scalar_type(100)
K2=default_scalar_type(10)
K11 = fem.Constant(domain,K1)
K22 = fem.Constant(domain,K2)
# Stored strain energy density
#psi_volumetrico_matrice=
psi_isochorico_matrice= (C_m/2)*(Ic1_dist -3)
psi_isochorico_fibre= K11/(2*K22)*(ufl.exp(K2*(Ic4_dist-1)**2)-1)
# Second Piola-Kirchhoff stress tensor
SPK_volumetrico_matrice = J*p*ufl.inv(C)
SPK_isochorico_matrice = 2 * J**(-2/3) * ufl.diff(psi_isochorico_matrice,Ic1_dist) * (I-(1/3)*ufl.inner(I,C)*ufl.inv(C))
SPK_isochorico_fibre = J**(-2/3) * ufl.diff(psi_isochorico_fibre,Ic4_dist) * (A0-(1/3)*ufl.inner(A0,C)*ufl.inv(C))
SPK_tot = SPK_volumetrico_matrice + SPK_isochorico_matrice + SPK_isochorico_fibre
# First Piola-Kirchhoff stress tensor
P = ufl.inv(FF)*SPK_tot
def Jacobian_F(u):
d = len(u)
# Identity tensor
I = ufl.Identity(d)
# Deformation gradient
FF = I + ufl.grad(u)
return ufl.variable(ufl.det(FF))
metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain, metadata=metadata)
#Per vito che leggerà nel futuro
#Il traction vector non è applicato su nessuna superficie
dx = ufl.Measure("dx", domain=domain, metadata=metadata) # Elementarily volume
# Define form F (we want to find u such that F(u) = 0)
def E(r):
d = len(r)
# Identity tensor
I = ufl.Identity(d)
# Deformation gradient
FF = ufl.variable(I + ufl.grad(r))
# Right Cauchy-Green tensor
C = ufl.variable(ufl.transpose(FF) * FF)
return ufl.variable(0.5*(C - I))
bulk_modulus = fem.Constant(domain, default_scalar_type(1000))
F1 = ufl.inner(E(δu) , SPK_volumetrico_matrice) * dx + ufl.inner(E(δu) , SPK_isochorico_matrice) * dx + ufl.inner(E(δu) , SPK_isochorico_fibre) * dx - ufl.inner(δu, B) * dx - ufl.inner(δu, T) * ds
F2 = (Jacobian_F(u)-j) *δp *dx
#U(JJ_trial) = bulk_modulus/2 * (JJ_trial - 1)**2
F3 = (bulk_modulus*(j - 1)-p)*δj*dx
F = F1 +F2 +F3
Is the variational formulation that I whore correct?
Any help is very appreciated