hi all,
im trying to create a rod with two diffenent subdomains, each for each module young.
when im trying to apply the wrap by vector in para view i see that the vectors are not consistent , any idea why? attaching relevant code snippets
length, width, height = 1.0, 0.2, 0.2 # Dimensions of the rod
num_elements = 12 # Number of elements along each dimension
E1, E2 = 1E8, 2E6 # Young's moduli for the two halves of the rod
# Create a 3D mesh for the rod
mesh = BoxMesh(Point(0, 0, 0), Point(length, width, height), num_elements, num_elements, num_elements)
# Define boundary
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= length/2 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > length/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)
DG = FunctionSpace(mesh, "DG", 0)
materials_function = Function(DG)
E_values = [E1, E2]
help = np.asarray(materials.array(), dtype=np.int32)
def assign_youngs_modulus(cell):
if materials[cell.index()] == 1:
return E1 # Young's modulus for the left half
else:
return E2 # Young's modulus for the right half
# Assign Young's modulus to each cell in the mesh
E_values = [assign_youngs_modulus(cell) for cell in cells(mesh)]
materials_function.vector()[:] = E_values
# Boundary conditions
class Boundaryleft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, DOLFIN_EPS)
class Boundaryright(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, DOLFIN_EPS)
V = VectorFunctionSpace(mesh, 'Lagrange', 2)
# ADDED
facets = MeshFunction("size_t", mesh, 2)
#facets.set_all(0)
Boundaryleft().mark(facets, 1)
Boundaryright().mark(facets, 2)
bc1 = DirichletBC(V.sub(0), Constant(0), facets, 1)
bc2 = DirichletBC(V.sub(0), Constant(1), facets, 2)
bc = [bc1, bc2]
then im using class to define the lamada and mu for each side of the rod and of course the sigma and epsilon
and then the solution
and this is my xdmf definition
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_solution, 0.0)
beam_deflection_file.write(u_magnitude, 0.0)
beam_deflection_file.write(sigma_11, 0.0)
beam_deflection_file.write(materials_function, 0.0)
beam_deflection_file.write(von_Mises_stress, 0.0)
this is the pic of how i see it
maybe its a normal behavior, i dont know