An easy to use immersed boundary method in 2D, with full implementations in MATLAB and Python that contains over 75 built-in examples, including multiple options for fiber-structure models and advection-diffusion, Boussinesq approximations, and/or artificial forcing.
from the please_find_lagrangnian_force_on_eulerian_grid.py, from line 821-864 is below
id_1 = pts_1[ii] # 1ST Node index
id_2 = pts_2[ii] # (MIDDLE) 2nd Node index -> index that gets force applied to it!
id_3 = pts_3[ii] # 3RD Node index
k_Beam = K_Vec[ii] # Beam stiffness of i-th spring
C = C_Vec[ii] # Curvature of the beam between these three nodes
Xp = xLag[id_1] # xPt of 1ST Node Pt. in beam
Xq = xLag[id_2] # xPt of 2ND (MIDDLE) Node Pt. in beam
Xr = xLag[id_3] # xPt of 3RD Node Pt. in beam
Yp = yLag[id_1] # yPt of 1ST Node Pt. in beam
Yq = yLag[id_2] # yPt of 2ND (MIDDLE) Node Pt. in beam
Yr = yLag[id_3] # yPt of 3RD Node Pt. in beam
# Checks if Lag. Pts. have passed through the boundary and translates appropriately
Xp,Xq,Xr = check_If_Beam_Points_Pass_Through_Boundary(ds,Lx,Xp,Xq,Xr)
Yp,Yq,Yr = check_If_Beam_Points_Pass_Through_Boundary(ds,Ly,Yp,Yq,Yr)
# Compute Cross-Product
cross_prod = (Xr-Xq)*(Yq-Yp) - (Yr-Yq)*(Xq-Xp)
# FORCES FOR LEFT NODE
bF_x_L = -k_Beam * ( cross_prod - C ) * ( Yr-Yq )
bF_y_L = k_Beam * ( cross_prod - C ) * ( Xr-Xq )
# FORCES FOR MIDDLE NODE
bF_x_M = k_Beam * ( cross_prod - C ) * ( (Yq-Yp) + (Yr-Yq) )
bF_y_M = -k_Beam * ( cross_prod - C ) * ( (Xr-Xq) + (Xq-Xp) )
# FORCES FOR RIGHT NODE
bF_x_R = -k_Beam * ( cross_prod - C ) * ( Yq-Yp )
bF_y_R = k_Beam * ( cross_prod - C ) * ( Xq-Xp )
fx[id_1] -= bF_x_L # Sum total forces for left node,
# in x-direction (this is LEFT node for this beam)
fy[id_1] -= bF_y_L # Sum total forces for left node,
# in y-direction (this is LEFT node for this beam)
fx[id_2] += bF_x_M # Sum total forces for middle node,
# in x-direction (this is MIDDLE node for this beam)
fy[id_2] += bF_y_M # Sum total forces for middle node,
# in y-direction (this is MIDDLE node for this beam)
fx[id_3] -= bF_x_R # Sum total forces for right node,
# in x-direction (this is RIGHT node for this beam)
fy[id_3] -= bF_y_R # Sum total forces for right node,
# in y-direction (this is RIGHT node for this beam)
it can be seen that bF_x_L + bF_x_M + bF_x_R=0, but it then adds the force to node as -bF_x_L, bF_x_M, -bF_x_R ?
The TOTAL force to add on the three nodes should be zero, but it doesn't here, is it a bug?
from the please_find_lagrangnian_force_on_eulerian_grid.py, from line 821-864 is below id_1 = pts_1[ii] # 1ST Node index id_2 = pts_2[ii] # (MIDDLE) 2nd Node index -> index that gets force applied to it! id_3 = pts_3[ii] # 3RD Node index k_Beam = K_Vec[ii] # Beam stiffness of i-th spring C = C_Vec[ii] # Curvature of the beam between these three nodes
it can be seen that bF_x_L + bF_x_M + bF_x_R=0, but it then adds the force to node as -bF_x_L, bF_x_M, -bF_x_R ? The TOTAL force to add on the three nodes should be zero, but it doesn't here, is it a bug?