HI ALL,
im setting a comined boundary conditions on a rode(direchlet and neumann)
for some reason i only set direchlet condition on x direction via para view it seems to be moving from the initial location and not staying in the same location , any idea why?
im using the code of the documantion of the hyperelastic material
when im setting the direclet condition on all on the directions its being clamped and not moving
attaching also screenshots of what i mean
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
import os
# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
E1, E2, nu = 1e5, 6e5, 0.35
# Create mesh and define function space
length, width, height = 1.0, 0.2, 0.2 # Dimensions of the rod
num_elements = 10 # Number of elements along each dimension
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, 4, 4)
#mesh = UnitCubeMesh(24, 16, 16)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= 1/2 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > 1/2 - DOLFIN_EPS
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
# Boundary conditions
boundary_parts = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
left = AutoSubDomain(lambda x: near(x[0], 0))
right = AutoSubDomain(lambda x: near(x[0], 1))
left.mark(boundary_parts, 1)
right.mark(boundary_parts, 2)
bc1 = DirichletBC(V.sub(0), Constant(0), boundary_parts, 1)
#bc2 = DirichletBC(V.sub(1), Constant(0), boundary_parts, 1)
#bc3 = DirichletBC(V.sub(2), Constant(0), boundary_parts, 1)
bcs = [bc1]
neumann = Expression(("0.1*x[0]"), degree=1)
ds = Measure('ds', domain=mesh, subdomain_data=boundary_parts)
class K(fe.UserExpression):
def __init__(self, materials, k_0, k_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return ()
E_global = K(materials, E1, E2, degree=0)
DG = FunctionSpace(mesh, "DG", 2)
materials_function = Function(DG)
E_values = [E1, E2]
materials_function = project(E_global, DG)
MU1, LAMBDA1 = Constant(E1/(2*(1 + nu))), Constant(E1*nu/((1 + nu)*(1 - 2*nu)))
MU2, LAMBDA2 = Constant(E2/(2*(1 + nu))), Constant(E2*nu/((1 + nu)*(1 - 2*nu)))
class MUMU(fe.UserExpression):
def __init__(self, materials, MU_0, MU_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.MU_0 = MU_0
self.MU_1 = MU_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.MU_0
else:
values[0] = self.MU_1
def value_shape(self):
return ()
MU = MUMU(materials, MU1, MU2, degree=0)
class lamabada(fe.UserExpression):
def __init__(self, materials, lam_0, lam_1, **kwargs):
super().__init__(**kwargs)
self.materials = materials
self.lam_0 = lam_0
self.lam_1 = lam_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 1:
values[0] = self.lam_0
else:
values[0] = self.lam_1
def value_shape(self):
return ()
LAMBDA = lamabada(materials, LAMBDA1, LAMBDA2, degree=0)
# Define functions
du = TrialFunction(V) # Incremental displacement
v = TestFunction(V) # Test function
u = Function(V) # Displacement from previous iteration
B = Constant((0.0, 0.0, 0.0)) # Body force per unit volume
T = Constant((0.0, 0.0, 0.0)) # Traction force on the boundary
#u.vector()[:] = 0.0
# Define an initial guess as an Expression
initial_guess = Expression(("x[0]", "0", "0"), degree=1)
# Project or interpolate the initial guess onto u
u.interpolate(initial_guess)
# Kinematics
d = len(u)
I = Identity(d) # Identity tensor
F = I + grad(u) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
CL = F*F.T
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# CS= MU*(C-I) # cauchy stress
# S= J*CS*inv(F).T
# Stored strain energy density (compressible neo-Hookean model)
psi = (MU/2)*(Ic - 3) - MU*ln(J) + (LAMBDA/2)*(ln(J))**2
n = FacetNormal(mesh)
# Total potential energy
Pi = psi*dx - dot(B, u)*dx - neumann * dot(n, u)*ds(2)
# Compute first variation of Pi (directional derivative about u in the direction of v)
F = derivative(Pi, u, v)
# Compute Jacobian of F
Je = derivative(F, u, du)
# Solve variational problem
solver_parameters = {'newton_solver': {'maximum_iterations': 100,
'absolute_tolerance': 1E-8,
'relative_tolerance': 1E-7}}
solve(F == 0, u, bcs, J=Je, solver_parameters=solver_parameters)
# Deviatoric part
dev_stress = MU * (CL - I)
# Volumetric part
vol_stress = LAMBDA * ln(J) * I
sigma = (1/J) * (dev_stress + vol_stress)
S = FunctionSpace(mesh, "P", 1) # Example function space
sigma_xx = project(sigma[0, 0], S)
points_of_interest = [(0.1, 0.1, 0.1), (0.5, 0.1, 0.1), (0.9, 0.1, 0.1)] # Example points
for point in points_of_interest:
try:
sigma_xx_value = sigma_xx(point)
print(f"sigma_xx at {point}: {sigma_xx_value}")
except Exception as e:
print(f"Could not evaluate sigma_xx at {point}: {e}")
# Save solution in VTK format
file = File("displacement.pvd");
file << u;
File("Materials.pvd")<< materials
vtkfile_subdomains = File('subdomains.pvd')
vtkfile_subdomains << materials
# Time-stepping loop
q = Constant(0.5)
dt = 0.25
dt = Constant(dt)
bF_magnitude = Constant(0.0)
t = 0
while t <= 5:
t += float(dt)
info("Time: {}".format(t))
# Increase traction
bF_magnitude.assign(100.0 * t)
# Prepare to solve and solve
solve(F == 0, u, bcs, J=Je, solver_parameters=solver_parameters)
# Close files
materials_function.rename("materials", "")
u.rename("Displacement Vector", "")
beam_deflection_file = fe.XDMFFile("beam_deflection.xdmf")
beam_deflection_file.parameters["flush_output"] = True
beam_deflection_file.parameters["functions_share_mesh"] = True
beam_deflection_file.write(u, 0.0)
beam_deflection_file.write(materials_function, 0.0)
# Plot solution
plot(u)
plt.show()
# Close the XDMF files
beam_deflection_file.close()
