### (Deleted) Cylinder Mesh and B.C. for an Elastic problem

121
views
1
3 months ago by
Hi Everyone,

I'm currently trying to write with FeniCS a Linear Elastic problem over a Cylindrical domain (Very thin Cylinder, with base radius 400 times the height, here i used 80 to be faster), problems is that i'm very new to this and i don't know many things that are going on.
So far i wrote as follow:

from fenics import *
from dolfin import *
from mshr import *

# Define 3D geometry
Cylinder = Cylinder(dolfin.Point(0, 0, 0.5), dolfin.Point(0, 0, -0.5), 40, 40)

generator = CSGCGALMeshGenerator3D()
generator.parameters["mesh_resolution"] = 15.0;
generator.parameters["cell_size"] = 1.0;

# Mesh generator
mesh = generator.generate(CSGCGALDomain3D(Cylinder))

dolfin.plot(mesh, "3D mesh")

vtkfile = File("cylinder.pvd")
vtkfile << (mesh)

few questions: Is correct to write the Cylinder as a geometry and then mesh it? If yes, i used this Mesh Generator but i don't know practically what's the "mesh resolution" and "cell size".. So i get a bad mesh and overall long to compute, i tried changing the number a few times but still they are pretty bad, if someone knows how to correct please could him/she inform me.
i was also wondering wha't the most efficent way to mesh a 3d domain like this? Is there any way automate the mesh in order to have a good mesh e.g. all the elements with the same size and connectivities all well made?

Thank you everyone,
Matteo

Community: FEniCS Project
3
For more fine control in meshing, use of a dedicated tool can relieve you of some headache:
gmsh: http://gmsh.info/
salome: https://www.salome-platform.org
written 3 months ago by klunkean

I managed to use gmsh to make the geometry and improrted the geometry in fenics through:

dolfin-convert geometry.msh geometry.xml

Now i have a few questions.. i'm having troubles defining the boundary conditions and i think the problem could be the mesh:

Are the mesh keep the same after it is imported? If yes, how is it possible to check ?
Cause it doens't look so (I attached code and pictures)
Moreover i'm trying to define Dirichlet Boundary Condition such as the geometry is fixed along the border loaded only under it self weight . If you have any advice it would be very welcome:

from fenics import *
from dolfin import *
from mshr import *
import numpy as np

#Geometries and Meshing with GMSH

mesh = Mesh("Membrane.xml");
print(mesh) # It return a print that doesn't look that the point (0,0,1) the center of the cylinder is in the right spot

vtkfile = File("cylinder.pvd")
vtkfile << (mesh)

# cd=MeshFunction('size_t', mesh, "Membrane.xml");  This function doesn't work (???)
# fd=MeshFunction('size_t', mesh);

# cd.size() is the total number of finite elements generated by the GMSH mesh.
# cd.size()

#Elasticity

V = VectorFunctionSpace(mesh, "P", 1)

#Constants

r = 400
mu = 1
rho = 1
beta = 1.25
lambda_ = beta
g = 9.81

#Def boundary condition, with tolerance

tol = 1E-14
w_D = Constant((0,0,0))
def clamped_boundary(x, on_boundary):
return on_boundary and np.abs(x[0]**2 + x[1]**2 - r**2) < tol

bc = DirichletBC(V, w_D, clamped_boundary)

# Define Strain

def epsilon(u):
def sigma(u):
return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# print(np.array(x[0]**2 + x[1]**2 -r**2) < tol)

# Define variational problem

u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0, 0, -rho*g))
d = u.geometric_dimension() #SPACE DIMENSION
# print(d)
T = Constant((0,0,0))

#GOVERNING EQUATION

a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution
V = FunctionSpace(mesh, "P", 1)
u_magnitude = sqrt(dot(u, u))
u_magnitude = project(u_magnitude, V)
plot(u_magnitude, title="’Displacement magnitude’")

vtkfile = File("displacement.pvd")
vtkfile << (u_magnitude)

THe geometry is a cylinder with a radius of 400 (Explicit in gmsh) but if i see the print(mesh) as here, the center is clearly moved along positive direction:

written 3 months ago by Matteo