version autonome

quartrou.py

@@ -0,0 +1,150 @@
+# -*- 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()

quartvierge.geo

@@ -0,0 +1,22 @@
+lx =
+ly =
+r =
+lc = 
+Point (1) = {0 ,0 ,0, lc}; // centre du trou
+Point (2) = {0 ,r ,0 ,lc};
+Point (3) = {0 ,ly ,0 ,lc};
+Point (4) = {lx ,ly ,0 ,lc};
+Point (5) = {lx ,0 ,0 ,lc};
+Point (6) = {r ,0 ,0 ,lc};
+Line (1) = {5 ,6}; // bas
+Line (2) = {4 ,5}; // droite
+Line (3) = {3 ,4}; // haut
+Line (4) = {2 ,3}; // gauche
+Circle (5) = {6 ,1 ,2};
+Line Loop (6) = {1 ,2 ,3 ,4 ,5};
+Plane Surface(7)={6};
+Physical Line ("bas") = {5 ,6}; // bas
+Physical Line ("droite") = {4 ,5}; // droite
+Physical Line ("haut") = {3 ,4}; // haut
+Physical Line ("gauche") = {2 ,3}; // gauche
+Physical Surface(1) ={7};