hi all,
im trying to run a loop that will calculate the stress and deformation at different point across my model, but for some reason it is not able to evaluate the calculation in the point
im getting different shape for each point, any idea why? im using the (shape,scheme,degree) option per Dokken’s recommendation, maybe im using it wrong.
im attaching part of my code, hopefully it will be enough
import fenics as fe
from dolfin import *
import matplotlib.pyplot as plt
from ffc.fiatinterface import create_quadrature as cquad
import numpy as np
import meshio
import dolfin
from fenics import *
import sys, os, shutil, math
import matplotlib.pyplot as plt # Plotting Library
F = variable(I + grad(u))
CL = F * F.T
# Deviatoric part
dev_stress = MU * (CL - I)
# Volumetric part
vol_stress = LAMBDA * (J - 1) * I
sigma = (1 / J) * dev_stress + vol_stress
S = FunctionSpace(mesh, "P", 1)
# Example function space
sigma_xx = project(sigma[0, 0], S)
stress123 = Function(S, name='stressi')
dt = 0.1
t, T = 0.0, 100 * dt
TT = Traction()
# creating functions in order to see the stress and deformation over time
V_tensor = TensorFunctionSpace(mesh, "DG", 0)
S = FunctionSpace(mesh, "P", 1)
print(f"V_tensor {V_tensor}")
defGrad = Function(V_tensor, name='F')
stress123 = Function(S, name='stressi')
shape = "interval"
deg =2
scheme = "default"
points1, weights = cquad(shape, deg, scheme)
print(f"weights{weights}")
# creating a file for paraview
stress123.rename("SRESS1", "")
subdomains_function.rename("subdomains", "")
defGrad.rename("Deformation_DRAD", "")
u.rename("Displacement Vector", "")
test_file = fe.XDMFFile("test.xdmf")
test_file.parameters["flush_output"] = True
test_file.parameters["functions_share_mesh"] = True
length = 1
absJ = 1
# Parameters
# Function to check if a point is outside both cylinders
beam_size = 0.2 # Beam side length
cylinder_radius_squared = 0.0182**2
num_points = 20
x_constant = 0.4 # The constant x value
# Function to check if a point is outside both cylinders
def is_outside_cylinders(y, z, y_center1, y_center2, cylinder_radius_squared):
distance1 = (y - y_center1)**2 + (z - 0.1)**2
distance2 = (y - y_center2)**2 + (z - 0.1)**2
return distance1 > cylinder_radius_squared and distance2 > cylinder_radius_squared
# Generate random points
def generate_points(num_points, beam_size, cylinder_radius_squared, x_constant):
points = []
while len(points) < num_points:
# Randomly sample within the beam's boundaries
y = np.random.uniform(0, beam_size)
z = np.random.uniform(0, beam_size)
# Define cylinder centers based on the constant x
y_center1 = 0.05 + 0.1 * x_constant
y_center2 = 0.15 - 0.1 * x_constant
# Check if the point is outside both cylinders
if is_outside_cylinders(y, z, y_center1, y_center2, cylinder_radius_squared):
points.append((x_constant, y, z))
return np.array(points)
# Generate 20 points
points = generate_points(num_points, beam_size, cylinder_radius_squared, x_constant)
# Output points
print("Generated points:")
print(points)
# Function to check if a point is inside both cylinders
def is_inside_cylinders(y, z, y_center1, y_center2, cylinder_radius_squared):
distance1 = (y - y_center1)**2 + (z - 0.1)**2
distance2 = (y - y_center2)**2 + (z - 0.1)**2
return distance1 <= cylinder_radius_squared or distance2 <= cylinder_radius_squared
# Generate random points inside the cylinders
def generate_points_inside_cylinders(num_points, beam_size, cylinder_radius_squared, x_constant):
physical_points2 = []
while len(physical_points2) < num_points:
# Randomly sample within the beam's boundaries
y = np.random.uniform(0, beam_size)
z = np.random.uniform(0, beam_size)
# Define cylinder centers based on the constant x
y_center1 = 0.05 + 0.1 * x_constant
y_center2 = 0.15 - 0.1 * x_constant
# Check if the point is inside either of the cylinders
if is_inside_cylinders(y, z, y_center1, y_center2, cylinder_radius_squared):
physical_points2.append((x_constant, y, z))
return np.array(physical_points2)
# Generate 20 points inside the cylinders
physical_points2 = generate_points_inside_cylinders(num_points, beam_size, cylinder_radius_squared, x_constant)
# Output points
print("Generated points inside cylinders:")
print(physical_points2)
# creating a loop with time dependent in order to apply force over time
while t <= T:
print('time: ', t)
# Increase traction
TT.t = t
tDirBC.time_ = t
tDirBC4.time_ = t
solver.solve()
F = variable(I + grad(u))
J = det(F)
CL = F * F.T
# Deviatoric part
dev_stress = MU * (CL - I)
# Volumetric part
vol_stress = LAMBDA * (J - 1) * I
sigma = (1 / J) * dev_stress + vol_stress
u.rename("u", "displacement")
DF = I + grad(u) # Compute the deformation gradient
defGrad.assign(project(DF, V_tensor)) # Project the deformation gradient
sigma_xx = project(sigma[0, 0], S)
# Extravgacting the xx component of the stress tensor
stress_xx_projected = project(sigma_xx, S) # Projecting this scalar component
stress123.assign(stress_xx_projected) # Assigning the projected stress
print(f"stress_xx_projected {stress_xx_projected}")
stress123.assign(stress_xx_projected) # Assigning the projected stress
# get stretch at a point for plotting
# trying to create avarage stress and deformation
for point in points:
try:
deformation_gradient_values_all11.append(defGrad(point)[0]) # Store the value at the current point
stress_values_all11.append(stress123(point))
stress_values_all111.append(np.dot(stress_values_all11, weights) * absJ)
deformation_gradient_values_all222.append(np.dot(deformation_gradient_values_all11, weights) * absJ)
except Exception as e:
print(f"Could not evaluate at {point}: {e}")
# Loop for physical_points2
for point in physical_points2:
try:
deformation_gradient_values_all4.append(defGrad(point)[0]) # Store the value at the current point
stress_values_all4.append(stress123(point))
stress_values_all5.append(np.dot(stress_values_all4, weights) * absJ)
deformation_gradient_values_all6.append(np.dot(deformation_gradient_values_all4, weights) * absJ)
except Exception as e:
print(f"Could not evaluate at {point}: {e}")
# Print the evaluated values for debugging purposes
print(stress_values_all111)
print(deformation_gradient_values_all222)
print(stress_values_all5)
print(deformation_gradient_values_all6)
deformation_gradient_values_all11 = []
stress_values_all11 = []
deformation_gradient_values_all4 = []
stress_values_all4 = []
# time increment
t += float(dt)
deformation_gradient_values_all333 = np.array(deformation_gradient_values_all222)
stress_values_all222 = np.array(stress_values_all111)
deformation_gradient_values_all7 = np.array(deformation_gradient_values_all6)
stress_values_all6 = np.array(stress_values_all5)
A_tag=0.0002624;
A_0=0.0394752;
stress_body=(((stress_values_all6*A_tag)+(stress_values_all222*A_0))/0.04)
deformation_gradient_body=np.array(((deformation_gradient_values_all7*A_tag)+(deformation_gradient_values_all333*A_0))/0.04)
this is what im getting
Could not evaluate at [0.4 0.01093385 0.13524289]: shapes (1,) and (2,) not aligned: 1 (dim 0) != 2 (dim 0)
Could not evaluate at [0.4 0.16278207 0.16510247]: shapes (3,) and (2,) not aligned: 3 (dim 0) != 2 (dim 0)
Could not evaluate at [0.4 0.07372184 0.19482677]: shapes (4,) and (2,) not aligned: 4 (dim 0) != 2 (dim 0)
Could not evaluate at [0.4 0.19559365 0.0247406 ]: shapes (5,) and (2,) not aligned: 5 (dim 0) != 2 (dim 0)
and so on...
i tried to change the scheme and degree
will appreciate any help