Dear,
I have two doubts:
Dear,
I have two doubts:
1st question - According to the example below, I have two volumes and I would like to consider two different types of mtarial in each volume. The model considers the individual volumes in gmsh using BooleanFragments. However, how do I consider each volume to be a material in the code below? I’m sorry, I did some research on discursse but I still can’t make this identification of two materials in the code.
2nd question - Is it possible in fenics to consider this model as a mechanism? Consider the contact between the two volumes? or consider some kind of rotation between them?
gsmh code:
//+
SetFactory("OpenCASCADE");
Box(1) = {0, 0, 0, 3, 1, 1};
//+
Cylinder(2) = {1, -0.2, 0.5, 0, 2.2, 0, 0.3, 2*Pi};
//+
BooleanFragments{ Volume{1}; Delete; }{ Volume{2}; Delete; }
//+
Physical Surface(20) = {11};
//+
Physical Surface(21) = {6};
//+
Physical Volume(22) = {1};
//+
Physical Volume(23) = {4, 2, 3};
Example code:
from fenics import *
import meshio
# Load
F = 2e2
# Define 3D geometry
mesh = Mesh()
with XDMFFile("mecanismo_mesh.xdmf") as infile:
infile.read(mesh)
mvc2 = MeshValueCollection("size_t", mesh, 2)
with XDMFFile("mecanismo_facet.xdmf") as infile:
infile.read(mvc2, "name_to_read")
mf = cpp.mesh.MeshFunctionSizet(mesh, mvc2)
mvc3 = MeshValueCollection("size_t", mesh, 3)
with XDMFFile("mecanismo_mesh.xdmf") as infile:
infile.read(mvc3, "name_to_read")
cf = cpp.mesh.MeshFunctionSizet(mesh, mvc3)
dx1 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=22)
dx2 = Measure("dx", domain=mesh, subdomain_data=cf, subdomain_id=23)
# Define function space and base functions
V = VectorFunctionSpace(mesh, 'Lagrange', degree=1)
# Boundary Condition
bc = DirichletBC(V, Constant((0.0, 0.0, 0.0)), mf, 20)
BCs = [bc]
# Load
ds = Measure("ds", domain=mesh, subdomain_data=mf, subdomain_id=21)
t = Expression(('q_x', 'q_y', 'q_z'), degree = 1, q_x = 0.0, q_y = 0.0, q_z = -F)
## Constitutive relation
def eps(v):
return sym(grad(v))
E = Constant(210e3)
nu = Constant(0.3)
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(v):
dim = v.geometric_dimension()
return lmbda*tr(eps(v))*Identity(dim) + 2.0*mu*eps(v)
## Variational formulation
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
L = dot(t, u_)*ds
A = assemble(a)
B = assemble(L)
[bc.apply(A, B) for bc in BCs]
u = Function(V, name="Displacement")
solve(A ,u.vector() ,B)
#file_results
file_results = XDMFFile("u_new_3D.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
