Direkt zum Inhalt | Direkt zur Navigation

UniBwM » LRT » LRT 1 » Prof. Gerdts » Teaching » Finite Elemente » Uebung7_A29_lsg

Uebung7_A29_lsg

from dolfin import *

#Gitter einlesen 
mesh=Mesh("bookshelf.xml");

#Def. FE-Raum
V=VectorFunctionSpace(mesh, 'CG', 10)

#Gitter auf Randfacetten getrennt ansprechen
boundary_parts=MeshFunction("uint", mesh, 1)

#Def. Materialdaten etc.
rho=0.0 #g/cm^3
g=9.81*(0.01)**2 #N/cm^2
E=1.0E6 #50.0E6*(0.01**2) #in N/cm^2
nu=0.2
G=2.0E2 #Druck in N/cm #11 kg

#Def. Dirichlet-Rand links
class Gamma_l(SubDomain):
   def inside(self, x, on_boundary):
       tol = 1E-14
       return on_boundary and x[0] < tol

GammaL=Gamma_l()
GammaL.mark(boundary_parts, 1)

u_l=Constant((0.0, 0.0))
bcd=DirichletBC(V,u_l,GammaL)

#Def. Neumann-Rand oben
class Gamma_o(SubDomain):
   def inside(self, x, on_boundary):
       tol = 1E-14
       return on_boundary and x[1] > 11.0-tol

GammaO=Gamma_o()
GammaO.mark(boundary_parts, 0)

#Variationsproblem
lambdaC=Constant(E*nu/((1.0 + nu)*(1.0 - 2.0*nu)))
muC=Constant(E/4.0/(1.0+nu))

u=TrialFunction(V)
v=TestFunction(V)

f=Constant((0.0,-rho*g))
g=Constant((0.0,-G))

def sigma(v):
    return lambdaC*tr(sym(grad(v)))*Identity(v.cell().d)+2.0*muC*sym(grad(v))

a=inner(sigma(u), sym(grad(v)))*dx
L=dot(f,v)*dx+dot(g,v)*ds(0)

#LGS loesen

u=Function(V)
solve(a == L, u, bcd)

#Visualisierung
plot(mesh)
plot(u, mode="displacement")
#file=File("Displ.pvd")
#file << u
interactive()