150 lines
3.5 KiB
Python
150 lines
3.5 KiB
Python
# -*- coding: utf-8 -*-
|
|
"""
|
|
Spyder Editor
|
|
|
|
This is a temporary script file.
|
|
"""
|
|
import fenics as fe
|
|
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
import dolfin
|
|
import os
|
|
|
|
# --------------------
|
|
# Functions and classes
|
|
# --------------------
|
|
#def bottom(x, on_boundary):
|
|
# return (on_boundary and fe.near(x[1], 0.0))
|
|
|
|
# Strain function
|
|
def epsilon(u):
|
|
return fe.sym(fe.grad(u))
|
|
|
|
# Stress function
|
|
def sigma(u):
|
|
return lambda_*fe.div(u)*fe.Identity(2) + 2*mu*epsilon(u)
|
|
|
|
|
|
# --------------------
|
|
# Parameters
|
|
# --------------------
|
|
|
|
# Density
|
|
rho = fe.Constant(200.0)
|
|
|
|
# Young's modulus and Poisson's ratio
|
|
E = 135e9
|
|
nu = 0.35
|
|
|
|
# Lame's constants
|
|
lambda_ = E*nu/(1+nu)/(1-2*nu)
|
|
mu = E/2/(1+nu)
|
|
|
|
l_x, l_y = 5.0, 5.0 # Domain dimensions
|
|
n_x, n_y = 20, 20 # Number of elements
|
|
lc = l_x/n_x
|
|
r = 2 # rayon du trou
|
|
|
|
# Load
|
|
f = 600
|
|
F = fe.Constant((0.0, f))
|
|
|
|
# Model type
|
|
model = "plane_strain"
|
|
if model == "plane_stress":
|
|
lambda_ = 2*mu*lambda_/(lambda_+2*mu)
|
|
#%% --------------------
|
|
# Geometry
|
|
# --------------------
|
|
#on adapte la géometrie aux paramètres
|
|
f1 = open("quartvierge.geo",'r')
|
|
toto = f1.read()
|
|
f1.close
|
|
toto = toto.replace("lx =","lx = %f;" % l_x)
|
|
toto = toto.replace("ly =","ly = %f;" % l_y)
|
|
toto = toto.replace("lc =","lc = %f;" % lc)
|
|
toto = toto.replace("r =","r = %f;" % r)
|
|
f2 = open("quart.geo",'w')
|
|
f2.writelines(toto)
|
|
f2.close()
|
|
# Commandes pour générer le maillage en appelant GMSH
|
|
#
|
|
command = 'gmsh -2 quart.geo -format msh2'
|
|
os.system(command)
|
|
#print ("Le maillage est fait!\n")
|
|
command = 'dolfin-convert quart.msh geometry.xml'
|
|
os.system(command)
|
|
#print ("Finished meshing!\n")
|
|
# demonstration du maillage
|
|
mesh = fe.Mesh("geometry.xml")
|
|
plt.figure(1)
|
|
fe.plot(mesh)
|
|
plt.show()
|
|
#%%
|
|
|
|
# Definition des limites
|
|
left = dolfin.CompiledSubDomain("on_boundary && near(x[0],0)")
|
|
right = dolfin.CompiledSubDomain("on_boundary && near(x[0], Lx)", Lx = l_x)
|
|
top = dolfin.CompiledSubDomain("on_boundary && near(x[1], Ly)", Ly = l_y)
|
|
bottom = dolfin.CompiledSubDomain("on_boundary && near(x[1], 0)")
|
|
|
|
boundary_indices = {"left": 0, "right": 1, "top": 2, "bottom": 3}
|
|
boundary_markers = dolfin.MeshFunction("size_t", mesh, dim=1, value=0)
|
|
left.mark(boundary_markers, boundary_indices["left"])
|
|
right.mark(boundary_markers, boundary_indices["right"])
|
|
right.mark(boundary_markers, boundary_indices["right"])
|
|
top.mark(boundary_markers, boundary_indices["top"])
|
|
bottom.mark(boundary_markers, boundary_indices["bottom"])
|
|
|
|
ds = dolfin.ds(domain=mesh, subdomain_data=boundary_markers)
|
|
dx = dolfin.dx(domain=mesh)
|
|
|
|
# --------------------
|
|
# Function spaces
|
|
# --------------------
|
|
V = fe.VectorFunctionSpace(mesh, "CG", 1)
|
|
u_tr = fe.TrialFunction(V)
|
|
u_test = fe.TestFunction(V)
|
|
|
|
# --------------------
|
|
# Boundary conditions
|
|
# --------------------
|
|
bc_l = fe.DirichletBC(V.sub(1), 0.0, left)
|
|
|
|
bc_b = fe.DirichletBC(V.sub(0), 0.0, bottom)
|
|
|
|
#bc = fe.DirichletBC(V, fe.Constant((0.0, 0.0)), bottom)
|
|
|
|
# --------------------
|
|
# Weak form
|
|
# --------------------
|
|
a = fe.inner(sigma(u_tr), epsilon(u_test))*fe.dx
|
|
l = fe.inner(F, u_test)*ds(boundary_indices["top"])
|
|
|
|
# --------------------
|
|
# Solver
|
|
# --------------------
|
|
u = fe.Function(V)
|
|
A_ass, L_ass = fe.assemble_system(a, l)
|
|
|
|
fe.solve(A_ass, u.vector(), L_ass)
|
|
|
|
print("jusqu'ici tout va bien")
|
|
|
|
#%% --------------------
|
|
# Post-process
|
|
# --------------------
|
|
|
|
plt.figure(2)
|
|
fe.plot(u.sub(0), title="Ux")
|
|
plt.show()
|
|
plt.figure(3)
|
|
fe.plot(u.sub(1), title="Uy")
|
|
plt.show()
|
|
plt.figure(4)
|
|
dolfin.plot(u, mode="displacement", title="displacement")
|
|
plt.show()
|
|
plt.figure(5)
|
|
stress = sigma(u)
|
|
fe.plot(stress[1, 1], title="Stress")
|
|
plt.show() |