Hello everybody,
for computing hyperelastic buckling problems I have set up the following code:
from __future__ import print_function
from mshr import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from dolfin import *
import matplotlib.pyplot as plt
import pandas as pd
if not has_petsc():
print("DOLFIN must be compiled at least with PETSc 3.6 to run this demo.")
exit(0)
parameters["linear_algebra_backend"] = "PETSc"
#========================================================Geometriewerte===================================================================
Ra_EBR=86.2/2000
Ri_EBR=75/2000
Ri_OR=75/2000
h_EBR=11/1000
d=8/1000
#======================================================2D-domain and mesh=================================================================
domain= Circle(Point(Ri_OR+d/2, h_EBR/2), d/2)
# Create mesh
mesh= generate_mesh(domain, 30)
mesh= refine(mesh)
# Create function space (aus ufl-UnifiedFormLanguage)
V= VectorFunctionSpace(mesh, "Lagrange", 1) #Legt ein Muster für einen 2x1-Vector auf jedem Punkt in mesh an
# Create solution, trial and test functions
u, du, v = Function(V), TrialFunction(V), TestFunction(V)
x = SpatialCoordinate(mesh)
#========================================================Definition des Verzerrungsverktors===============================================
# Definition of eps as a function -ufl
def eps(u):
return sym(as_tensor([[u[0].dx(0) , 0 , 0.5*(u[1].dx(0)+u[0].dx(1)) ],
[0 , u[0]/x[0] , 0 ],
[0.5*(u[1].dx(0)+u[0].dx(1)) , 0 , u[1].dx(1) ]]))
#=========================================================Definition der Randbedingungen==================================================
# Box constrains
constraint_l = Expression(("x_min-x[0]", "y_min-x[1]"), x_min=Ri_EBR, y_min=0, degree=2)
constraint_u = Expression(("x_max-x[0]", "y_max-x[1]"), x_max=Ra_EBR, y_max=h_EBR, degree=2)
u_min = interpolate(constraint_l, V)
u_max = interpolate(constraint_u, V)
# Surface force
f = Constant((0, 0))
#=============================================================Material Model==============================================================
E, nu = 8* 10**6, 0.4
mu = Constant(E/(2.0*(1.0+nu)))
lmbda = Constant(E*nu/((1.0+nu)*(1.0-2.0*nu)))
# Compressible neo-Hookean model
I = Identity(3)
F = I + eps(u)
C = F.T*F
Ic = tr(C)
J = det(F)
psi = (mu/2)*(Ic-3)-mu*ln(J)+(lmbda/2)*(ln(J))**2
#======================================================= Variational formulation==========================================================
elastic_energy = psi*dx - dot(f, u)*ds
grad_elastic_energy = derivative(elastic_energy, u, v)
H_elastic_energy = derivative(grad_elastic_energy, u, du)
#===============================================================Solver====================================================================
# Define the minimisation problem by using OptimisationProblem class -dolfin
class BucklingProblem(OptimisationProblem):
def __init__(self):
OptimisationProblem.__init__(self)
# Objective function
def f(self, x):
u.vector()[:] = x
return assemble(elastic_energy)
# Gradient of the objective function
def F(self, b, x):
u.vector()[:] = x
assemble(grad_elastic_energy, tensor=b)
# Hessian of the objective function
def J(self, A, x):
u.vector()[:] = x
assemble(H_elastic_energy, tensor=A)
# Create the PETScTAOSolver
solver = PETScTAOSolver()
# Set some parameters
#solver.parameters["linear_solver"] = "nash"
solver.parameters["method"] = "bqpip"
solver.parameters["monitor_convergence"] = True
solver.parameters["gradient_absolute_tol"] = 1.0e-1
solver.parameters["gradient_relative_tol"] = 1.0e-06
solver.parameters["report"] = False
solver.parameters["maximum_iterations"] = 1000
# Parse (PETSc) parameters
parameters.parse()
# Solve the problem
solver.solve(BucklingProblem(), u.vector(), u_min.vector(), u_max.vector())
#===================================================Transformation des Ergebnisvektors====================================================
num_dofs_per_component = int(V.dim()/V.num_sub_spaces()) #Menge der Knotenpunkte
num_sub_spaces = V.num_sub_spaces() #Unterräume...hier gleichbedeutend mit Koordinaten
uv = np.zeros((num_sub_spaces, num_dofs_per_component)) #Aufbau einer Matrix der Größe Koordinaten X Knotenpunkte
for i in range(num_sub_spaces):
uv[i] = u.sub(i, deepcopy=True).vector().get_local()
rz = V.sub(0).collapse().tabulate_dof_coordinates()
uv = uv.T
ur = np.zeros(num_dofs_per_component)
uz = np.zeros(num_dofs_per_component)
r = np.zeros(num_dofs_per_component)
z = np.zeros(num_dofs_per_component)
i=0
for coord, vec in zip(rz, uv):
ur[i] = vec[0]
uz[i] = vec[1]
r[i] = coord[0]
z[i] = coord[1]
i = i+1
r_end=r+ur
z_end=z+uz
get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x >= y]
ind = get_indexes(Ri_EBR+(Ra_EBR-Ri_EBR)/1000,r_end)
z_i0 = z[ind]
r_i0 = r[ind]
z_i = z_end[ind]
r_i = r_end[ind]
get_indexes = lambda x, xs: [i for (y, i) in zip(xs, range(len(xs))) if x <= y]
ind = get_indexes(Ra_EBR-(Ra_EBR-Ri_EBR)/1000,r_end)
z_a0 = z[ind]
r_a0 = r[ind]
z_a = z_end[ind]
r_a = r_end[ind]
#========================================================Berechnen der mechanischen Spannung==============================================
W = TensorFunctionSpace(mesh, 'DG', 1, shape=(3,3)) #analog V
sig_ = 2*mu*sym(eps(u)) + (lmbda*tr(eps(u)))*Identity(3)
sig_w = project(sig_, W)
[s00, s01, s02, s10, s11, s12, s20, s21, s22] = sig_w.split(True)
sig_i=np.zeros(len(r_i))
for i in range (len(r_i)):
sig_i[i]=s00(r_i0[i],z_i0[i])
sig_a=np.zeros(len(r_a))
for i in range (len(r_a)):
sig_a[i]=s00(r_a0[i],z_a0[i])
#===============================================Plotten der Erebnisse/Speichern der Plots=================================================
EBRx= np.array([Ri_EBR, Ri_EBR, Ra_EBR, Ra_EBR, Ri_EBR])
EBRy= np.array([0, h_EBR, h_EBR, 0, 0])
plt.figure()
plot(mesh, linewidth=0.8)
plt.plot(EBRx, EBRy, linewidth=1, color='k', linestyle='-')
plt.axis('equal')
plt.grid(color='k', linestyle='-', linewidth=0.2)
plt.ylabel('x [m]')
plt.xlabel('r [m]')
plt.title('Vernetzung')
plt.savefig('Vernetzung.png')
#plt.show()
EBRx =EBRx *1000
EBRy =EBRy *1000
plt.figure()
plt.plot(Ri_EBR*1000+(sig_i/max(abs(sig_a)))*500*(Ra_EBR-Ri_EBR), z_i*1000, '.') #*10**(-6)
plt.plot(Ra_EBR*1000-(sig_a/max(abs(sig_a)))*500*(Ra_EBR-Ri_EBR), z_a*1000, '.')
plt.plot(EBRx, EBRy, linewidth=1, color='k', linestyle='dotted')
epsilon=plt.scatter((r+ur)*1000,(z+uz)*1000,c=ur*1000, s=0.5)
plt.arrow(Ra_EBR*1000, h_EBR*1000, 600*(Ra_EBR-Ri_EBR), 0, width = 0.02)
plt.text(Ra_EBR*1000+600*(Ra_EBR-Ri_EBR), h_EBR*1050, 'p/p_max', horizontalalignment='center', verticalalignment='center')
plt.text(Ra_EBR*1000+500*(Ra_EBR-Ri_EBR), h_EBR*1050, '1', horizontalalignment='center', verticalalignment='center')
plt.text(Ra_EBR*1000+250*(Ra_EBR-Ri_EBR), h_EBR*1050, '1/2', horizontalalignment='center', verticalalignment='center')
plt.text(Ra_EBR*1000, h_EBR*1050, '0', horizontalalignment='center', verticalalignment='center')
plt.plot([Ra_EBR*1000+500*(Ra_EBR-Ri_EBR), Ra_EBR*1000+500*(Ra_EBR-Ri_EBR)], [h_EBR*1000, h_EBR*0], linewidth=0.5, color='k', linestyle='-')
plt.plot([Ra_EBR*1000+250*(Ra_EBR-Ri_EBR), Ra_EBR*1000+250*(Ra_EBR-Ri_EBR)], [h_EBR*1000, h_EBR*0], linewidth=0.5, color='k', linestyle='-')
plt.plot([Ra_EBR*1000, Ra_EBR*1000], [h_EBR*1000, 0], linewidth=0.7, color='k', linestyle='-')
plt.arrow(Ri_EBR*1000, h_EBR*1000, -600*(Ra_EBR-Ri_EBR), 0, width = 0.02)
plt.text(Ri_EBR*1000-600*(Ra_EBR-Ri_EBR), h_EBR*1050, 'p/p_max', horizontalalignment='center', verticalalignment='center')
plt.text(Ri_EBR*1000-500*(Ra_EBR-Ri_EBR), h_EBR*1050, '1', horizontalalignment='center', verticalalignment='center')
plt.text(Ri_EBR*1000-250*(Ra_EBR-Ri_EBR), h_EBR*1050, '1/2', horizontalalignment='center', verticalalignment='center')
plt.text(Ri_EBR*1000, h_EBR*1050, '0', horizontalalignment='center', verticalalignment='center')
plt.plot([Ri_EBR*1000-500*(Ra_EBR-Ri_EBR), Ri_EBR*1000-500*(Ra_EBR-Ri_EBR)], [h_EBR*1000, h_EBR*0], linewidth=0.5, color='k', linestyle='-')
plt.plot([Ri_EBR*1000-250*(Ra_EBR-Ri_EBR), Ri_EBR*1000-250*(Ra_EBR-Ri_EBR)], [h_EBR*1000, h_EBR*0], linewidth=0.5, color='k', linestyle='-')
plt.plot([Ri_EBR*1000, Ri_EBR*1000], [h_EBR*1000, 0], linewidth=0.7, color='k', linestyle='-')
plt.axis('equal')
plt.grid(b=None, which='major', axis='both')
plt.ylabel('x [mm]')
plt.xlabel('r [mm]')
cbar=plt.colorbar(epsilon)
cbar.set_label('Verschiebung [mm]', rotation=270)
plt.title('Verschiebung')
plt.savefig('Verschiebung.png')
plt.show()
(in fact, the dolfin example by Tianyi Li for hyperelastic beams and some other examples did the main work. 
)
which works completely fine, but now I want to use an polygonal domain for the goal domain.
What I tried is to define the boundary conditions as functions.
def constraint_u(self, x):
if x[0] > Ra_EBR-1.5/1000:
return Expression(("x_max-x[0]", "y_max-x[1]"), x_max=Ra_EBR, y_max=h_EBR, degree=2)
else:
return Expression(("x_max-x[0]", "y_max-x[1]"), x_max=Ra_EBR, y_max=h_EBR-6.8/1000, degree=2)
But beside the fact, that this is very complicated to set up for different domain geometries, it doesn’t work, because i don’t know how to involve this function into the solver.
Does anyone have an advice how to realize it? Or, what would be even better, how to involve an defined polygonal domain instead of u_min and u_max. A domain defined like
domain_vertices = [Point(Ri_EBR,0),
Point(Ra_EBR,0),
Point(Ra_EBR,h_EBR*3/4),
Point(Ra_EBR-(Ra_EBR-Ri_EBR)/2,h_EBR*3/4),
Point(Ra_EBR-(Ra_EBR-Ri_EBR)/2,h_EBR),
Point(Ri_EBR,h_EBR),
Point(Ri_EBR,0)]
domain= Polygon(domain_vertices)
Thank you very much.