- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
I asked GOOGLE Ai for a Fortran program, actually three, just to see if it was easy hard or just plain dumb.
It is just plain dumb.
This is a 3D plastic beam analysis program, it works until you ask it to include the Pardiso solver and then it tanks. I know how to do it with an older program just steal the Fortran code, but I was wondering about AI.
Steve: Can you make this work without relying on @mecej4 rather brilliant solution. I cannot not. I wish I could ask Jim Dempsey or mecej4 as well but they have disappeared from the forum.
! =====================================================================
! 3D Space Frame Plastic Solver with Global Sparse Matrix Assembly,
! CSR Pointer Array Mapping, and an Intel MKL PARDISO Solver Engine
! =====================================================================
module global_fea_mod
use mkl_pardiso
implicit none
private
public :: element_type, setup_element, assemble_sparse_triplets, solve_sparse_pardiso
! System topology sizing limits
integer, parameter, public :: max_nodes = 3
integer, parameter, public :: dof_per_node = 6
integer, parameter, public :: total_g_dof = max_nodes * dof_per_node
integer, parameter, public :: max_triplets = 12 * 12 * 2 ! Max non-zeros across 2 elements
type :: element_type
integer :: id
integer :: node_i, node_j ! Node connectivity indices
integer :: LM(12) ! Location Matrix mapping local DoFs to global DoFs
real(8) :: length ! Element physical length
real(8) :: cx, cy, cz ! Transformation direction cosines
real(8) :: P, My1, Mz1, My2, Mz2 ! Active internal element forces
real(8) :: Py, My_yield ! Cross-sectional plastic limits
real(8) :: k_factor ! Tangent stiffness softening modifier
real(8) :: ratio ! Yield surface interaction ratio
logical :: is_plastic
end type element_type
contains
! Initialize element parameters and build its unique LM mapping array
subroutine setup_element(elem, id, ni, nj, xi, yi, zi, xj, yj, zj, py, my_y)
type(element_type), intent(out) :: elem
integer, intent(in) :: id, ni, nj
real(8), intent(in) :: xi, yi, zi, xj, yj, zj, py, my_y
real(8) :: dx, dy, dz
integer :: d
elem%id = id
elem%node_i = ni
elem%node_j = nj
elem%Py = py
elem%My_yield = my_y
elem%P = 0.0d0; elem%My1 = 0.0d0; elem%Mz1 = 0.0d0; elem%My2 = 0.0d0; elem%Mz2 = 0.0d0
elem%k_factor = 1.0d0; elem%ratio = 0.0d0; elem%is_plastic = .false.
dx = xj - xi; dy = yj - yi; dz = zj - zi
elem%length = sqrt(dx**2 + dy**2 + dz**2)
elem%cx = dx / elem%length; elem%cy = dy / elem%length; elem%cz = dz / elem%length
! Populate Location Matrix (LM) destination vector fields
do d = 1, 6
elem%LM(d) = (ni - 1) * dof_per_node + d ! Node I Global DoFs
elem%LM(d + 6) = (nj - 1) * dof_per_node + d ! Node J Global DoFs
end do
end subroutine setup_element
! Formulate elemental Ke + Kg, and generate index coordinate tuples for global assembly
subroutine assemble_sparse_triplets(elem, trip_count, rows, cols, values)
type(element_type), intent(inout) :: elem
integer, intent(inout) :: trip_count
integer, intent(inout) :: rows(max_triplets), cols(max_triplets)
real(8), intent(inout) :: values(max_triplets)
real(8) :: Ke(12,12), Kg(12,12), Kt(12,12)
real(8) :: EA, EI, L, P
integer :: r, c, global_r, global_c
L = elem%length
P = elem%P
EA = 2.0d5
EI = 4.0d3
! Update localized cross-section material plasticity ratios
elem%ratio = (elem%P / elem%Py)**2 + (max(abs(elem%My1), abs(elem%My2)) / elem%My_yield)**2
if (elem%ratio >= 1.0d0) then
elem%is_plastic = .true.
elem%k_factor = 0.02d0 ! Post-yield tangent stiffness drop
else
elem%is_plastic = .false.
elem%k_factor = 1.0d0 ! Fully elastic
end if
Ke = 0.0d0; Kg = 0.0d0
! 1. Local Elastic Tangent Stiffness Matrix (Ke)
Ke(1,1) = (EA / L) * elem%k_factor; Ke(1,7) = -(EA / L) * elem%k_factor
Ke(7,1) = -(EA / L) * elem%k_factor; Ke(7,7) = (EA / L) * elem%k_factor
Ke(5,5) = (4.0d0 * EI / L) * elem%k_factor; Ke(5,11) = (2.0d0 * EI / L) * elem%k_factor
Ke(11,5) = (2.0d0 * EI / L) * elem%k_factor; Ke(11,11) = (4.0d0 * EI / L) * elem%k_factor
! 2. Local Geometric Stiffness matrix formulation (Kg) for axial P-Delta tracking
Kg(5,5) = (2.0d0 * P * L / 15.0d0); Kg(5,11) = -(P * L / 30.0d0)
Kg(11,5) = -(P * L / 30.0d0); Kg(11,11) = (2.0d0 * P * L / 15.0d0)
Kt = Ke + Kg
! 3. Map sparse coordinate rows/cols using Location Matrix (LM) pointers
do r = 1, 12
global_r = elem%LM(r)
do c = 1, 12
global_c = elem%LM(c)
if (global_r <= total_g_dof .and. global_c <= total_g_dof) then
trip_count = trip_count + 1
rows(trip_count) = global_r
cols(trip_count) = global_c
values(trip_count) = Kt(r, c)
end if
end do
end do
end subroutine assemble_sparse_triplets
! Modernized Intel MKL PARDISO Engine wrapper
subroutine solve_sparse_pardiso(n, nnz, ia, ja, a, rhs, displacement)
integer, intent(in) :: n, nnz
integer, intent(in) :: ia(n+1), ja(nnz)
real(8), intent(in) :: a(nnz)
real(8), intent(inout) :: rhs(n)
real(8), intent(out) :: displacement(n)
! PARDISO control variables
type(MKL_PARDISO_HANDLE) :: pt(64)
integer :: maxfct, mnum, mtype, phase, nrhs, error, msglvl
integer :: iparm(64)
integer :: idum(1)
real(8) :: ddum(1)
integer :: i
! Initialize PARDISO internal state pointers
do i = 1, 64
pt(i)%DUMMY = 0
iparm(i) = 0
end do
! Set PARDISO parameters for real structurally symmetric indefinite arrays
iparm(1) = 1 ! Explicit parameters configuration (no defaults)
iparm(2) = 2 ! Parallel fill-in reordering from METIS
iparm(3) = 0 ! Automatically query environment OpenMP threads
iparm(4) = 0 ! Direct solver execution (no iterative refinement)
iparm(6) = 0 ! Output solution directly into displacement array
iparm(10) = 13 ! Pivot perturbation tolerance threshold
iparm(11) = 1 ! Enable scaling/equilibration vectors
iparm(18) = -1 ! Output number of non-zero factorization profiles
iparm(21) = 1 ! Bunch-Kaufman pivoting for indefinite structures
mtype = -2 ! Real structurally symmetric matrix configuration
nrhs = 1 ! Number of right-hand side force vectors
maxfct = 1 ! Maximum number of numerical factorizations
mnum = 1 ! Matrix number in solver memory stack
msglvl = 0 ! Silent diagnostic messaging flag (set to 1 to debug)
! Step 1: Symbolic Analysis / Reordering Phase
phase = 11
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
if (error /= 0) then
print *, "[PARDISO ERROR] Symbolic analysis failure code: ", error
return
end if
! Step 2: Numerical Factorization Phase
phase = 22
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
if (error /= 0) then
print *, "[PARDISO ERROR] Matrix factorization failure code: ", error
return
end if
! Step 3: Back Substitution and Solution Vector Evaluation Phase
phase = 33
call pardiso(pt, maxfct, mnum, mtype, phase, n, a, ia, ja, idum, nrhs, iparm, msglvl, rhs, displacement, error)
if (error /= 0) then
print *, "[PARDISO ERROR] Back-substitution execution failure code: ", error
return
end if
! Step 4: Solver Memory Release Cleanup Phase
phase = -1
call pardiso(pt, maxfct, mnum, mtype, phase, n, ddum, ia, ja, idum, nrhs, iparm, msglvl, ddum, ddum, error)
end subroutine solve_sparse_pardiso
end module global_fea_mod
! Main structural solver execution routine
program main
use global_fea_mod
implicit none
type(element_type) :: framework(2)
integer :: csv_unit, step, i, j, r, c
integer, parameter :: total_steps = 50
! Sparse triplet tracking structures
integer :: triplet_count
integer :: trip_rows(max_triplets), trip_cols(max_triplets)
real(8) :: trip_vals(max_triplets)
! Final CSR storage matrices passed into PARDISO
integer :: csr_ia(total_g_dof + 1)
integer, allocatable :: csr_ja(:)
real(8), allocatable :: csr_a(:)
real(8) :: dense_temp(total_g_dof, total_g_dof)
real(8) :: F_global(total_g_dof)
real(8) :: U_global(total_g_dof)
! Displacement control tracking parameters
real(8) :: lambda, reference_P, reference_M, target_displacement_increment
integer :: active_nnz
! Node 1: (0,0,0) [Fixed base], Node 2: (4,0,0) [Joint], Node 3: (8,3,0) [Free end]
call setup_element(framework(1), 1, 1, 2, 0.0d0, 0.0d0, 0.0d0, 4.0d0, 0.0d0, 0.0d0, 600.0d0, 90.0d0)
call setup_element(framework(2), 2, 2, 3, 4.0d0, 0.0d0, 0.0d0, 8.0d0, 3.0d0, 0.0d0, 450.0d0, 70.0d0)
reference_P = 20.0d0
reference_M = 4.5d0
target_displacement_increment = 0.005d0 ! Prescribe fixed 5mm structural movement steps
lambda = 0.0d0
open(newunit=csv_unit, file='global_assembly_history.csv', status='replace', action='write')
write(csv_unit, '(A)') "Step,Lambda,ElemID,Length,P_Internal,My_Internal,Node2_U_Disp,Node2_V_Rot,YieldRatio,State"
print *, "====================================================================================="
print *, " 3D ADAPTIVE FEA FRAME ENGINE - INTEL MKL PARDISO SPARSE SOLVER "
print , "====================================================================================="
write(, '(A5, A8, A4, A12, A12, A12, A8)') "Step", "Lambda", "ID", "P-Internal", "Node2-Disp", "Yield-R", "Status"
print *, "-------------------------------------------------------------------------------------"
! Execution step loop
do step = 1, total_steps
triplet_count = 0
F_global = 0.0d0
dense_temp = 0.0d0
! 1. Assemble elements into sparse triplet coordinate vectors
do i = 1, 2
call assemble_sparse_triplets(framework(i), triplet_count, trip_rows, trip_cols, trip_vals)
! 2. Compress triplets onto a temporary matrix to eliminate duplicate index slots
do i = 1, triplet_count
dense_temp(trip_rows(i), trip_cols(i)) = dense_temp(trip_rows(i), trip_cols(i)) + trip_vals(i)
end do
! 3. ADAPTIVE DISPLACEMENT CONTROL BLOCK: Evaluate active diagonal tangent stiffness
if (dense_temp(7,7) > 1.0d-4) then
lambda = lambda + (dense_temp(7,7) * target_displacement_increment) / reference_P
else
lambda = lambda + 0.01d0
end if
! Map loads based on adaptive load factors
F_global(7) = reference_P * lambda
F_global(11) = reference_M * lambda
! 4. Convert structural dense components into standard Compressed Sparse Row (CSR) format
active_nnz = 0
do r = 1, total_g_dof
do c = 1, total_g_dof
if (abs(dense_temp(r,c)) > 1.0d-10) then
active_nnz = active_nnz + 1
end if
end do
end do
end do
if (allocated(csr_ja)) deallocate(csr_ja)
if (allocated(csr_a)) deallocate(csr_a)
allocate(csr_ja(active_nnz), csr_a(active_nnz))
active_nnz = 1
do r = 1, total_g_dof
csr_ia(r) = active_nnz
do c = 1, total_g_dof
if (abs(dense_temp(r,c)) > 1.0d-10) then
csr_a(active_nnz) = dense_temp(r,c)
csr_ja(active_nnz) = c
active_nnz = active_nnz + 1
end if
end do
end do
csr_ia(total_g_dof + 1) = active_nnz
! 5. Call the multi-threaded Intel PARDISO Direct Sparse Solver module
call solve_sparse_pardiso(total_g_dof, active_nnz - 1, csr_ia, csr_ja, csr_a, F_global, U_global)
! Post-convergence state processing
do i = 1, 2
framework(i)%P = reference_P * lambda
framework(i)%My1 = reference_M * lambda
write(csv_unit, '(I0,A,F6.2,A,I0,A,F6.2,A,F7.1,A,F7.1,A,E11.4,A,E11.4,A,F6.2,A,A)') &step, ',', lambda, ',', framework(i)%id, ',', framework(i)%length, ',', &framework(i)%P, ',', framework(i)%My1, ',', U_global(7), ',', U_global(11), ',', &framework(i)%ratio, ',', merge("PLASTIC", "ELASTIC", framework(i)%is_plastic)
write(*, '(I4, F8.2, I4, F12.1, E12.3, F12.2, 2X, A7)') step, lambda, framework(i)%id, framework(i)%P, &U_global(7), framework(i)%ratio, &merge("PLASTIC", "ELASTIC", framework(i)%is_plastic)
end do
end do
close(csv_unit)
print *, "-------------------------------------------------------------------------------------"
print *, "PARDISO execution complete! Sparse matrix tracking saved to 'global_assembly_history.csv'."
end program main
!****************************************************************************
!
! PROGRAM: Sulphur
!
! PURPOSE: Entry point for the console application.
!
! Technical Implementation Details
! Passivation ModelingMechanism:
! As solid {CaSO}_{4}) precipitates, it adheres to the solid {CaCO}_{3} particle core.
! Equation: passivation_factor = exp(-beta * state(3))
! Behavior: The reaction slows down exponentially relative to the product mass.
! It mimics physical shell blinding.4th-Order Runge-Kutta
! SolverPrecision: RK4 cuts the localized truncation error down to (O(dt^5)).
! Process: It samples the derivative curves four distinct times across each time step segment.
! Stability: It avoids numerical oscillations or artificial negative masses when a reaction terminates abruptly.
!****************************************************************************
program caco3_h2so4_rk4_passivation
implicit none
! Molar masses (g/mol)
real, parameter :: MM_CaCO3 = 100.09
real, parameter :: MM_H2SO4 = 98.08
real, parameter :: MM_CaSO4 = 136.14
real, parameter :: MM_CO2 = 44.01
! Physical and Thermodynamic Constants
real, parameter :: R = 8.314
real, parameter :: A_freq = 1.2e4
real, parameter :: Ea = 35000.0
real, parameter :: volume = 1.0
! Passivation Parameter
! Higher values simulate faster surface choking by insoluble CaSO4
real, parameter :: beta = 0.08
! Simulation Settings
real :: temperature
real :: t, dt, t_max
integer :: step, total_steps
integer, parameter :: file_unit = 15
! State Vector: [1]=CaCO3, [2]=H2SO4, [3]=CaSO4, [4]=CO2
real :: y(4)
real :: k1(4), k2(4), k3(4), k4(4), y_tmp(4), dydt(4)
real :: k_arrhenius
! 1. Initialize State and Conditions
temperature = 298.15 ! 25 °C
y(1) = 5000.0 ! Initial CaCO3 (g)
y(2) = 49.0 ! Initial H2SO4 (g)
y(3) = 0.0 ! Initial CaSO4 (g)
y(4) = 0.0 ! Initial CO2 (g)
! Time Configuration
t = 0.0
dt = 0.1 ! RK4 allows a larger, stable time step than Euler
t_max = 60.0
total_steps = int(t_max / dt)
! Calculate global Arrhenius constant
k_arrhenius = A_freq * exp(-Ea / (R * temperature))
! 2. Setup File Export
open(unit=file_unit, file='passivation_kinetics.csv', status='replace', action='write')
write(file_unit, '(A)') "Time_s,CaCO3_g,H2SO4_g,CaSO4_g,CO2_g"
print '(A, F6.2, A)', "--> Simulating Passivation with RK4 at: ", temperature - 273.15, " C"
print '(A10, A12, A12, A12, A12)', "Time(s)", "CaCO3(g)", "H2SO4(g)", "CaSO4(g)", "CO2(g)"
print '(A)', "------------------------------------------------------------"
! 3. RK4 Integration Loop
do step = 0, total_steps
! Export current state to CSV
write(file_unit, '(F8.2, A1, F10.3, A1, F10.3, A1, F10.3, A1, F10.3)') &
t, ',', y(1), ',', y(2), ',', y(3), ',', y(4)
! Print to console every 5 seconds
if (mod(step, 50) == 0) then
print '(F10.2, F12.2, F12.2, F12.2, F12.2)', t, y(1), y(2), y(3), y(4)
end if
! Early exit conditions
if (y(1) <= 0.001 .or. y(2) <= 0.001) then
print '(A, F5.2, A)', "--> Reaction stopped at ", t, " s (reactant depletion or passivation halt)."
exit
end if
! RK4 Step Evaluation
! Evaluation 1
call get_derivatives(y, dydt, k_arrhenius, beta, volume)
k1 = dydt * dt
! Evaluation 2
y_tmp = max(0.0, y + 0.5 * k1)
call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
k2 = dydt * dt
! Evaluation 3
y_tmp = max(0.0, y + 0.5 * k2)
call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
k3 = dydt * dt
! Evaluation 4
y_tmp = max(0.0, y + k3)
call get_derivatives(y_tmp, dydt, k_arrhenius, beta, volume)
k4 = dydt * dt
! Final Weighted Update
y = max(0.0, y + (k1 + 2.0*k2 + 2.0*k3 + k4) / 6.0)
t = t + dt
end do
close(file_unit)
print '(A)', "------------------------------------------------------------"
print '(A)', "--> Simulation complete. Saved to 'passivation_kinetics.csv'"
contains
! Subroutine to calculate derivatives for the state vector
subroutine get_derivatives(state, rates, k_arr, b_coeff, vol)
real, intent(in) :: state(4)
real, intent(out) :: rates(4)
real, intent(in) :: k_arr, b_coeff, vol
real :: c_h2so4, r_molar, passivation_factor
! Safeguard against negative mass states during intermediate RK4 steps
if (state(1) <= 0.0 .or. state(2) <= 0.0) then
rates = 0.0
return
end if
! Calculate acid concentration (mol/L)
c_h2so4 = (state(2) / MM_H2SO4) / vol
! Passivation Model: Exponential decay of active surface area due to solid CaSO4 build-up
passivation_factor = exp(-b_coeff * state(3))
! Total Molar Rate combining Shrinking Core and Passivation Blockage
r_molar = k_arr * (state(1)**(2.0/3.0)) * c_h2so4 * passivation_factor
! Convert to Mass Rates (g/s)
rates(1) = -r_molar * MM_CaCO3 ! d(CaCO3)/dt
rates(2) = -r_molar * MM_H2SO4 ! d(H2SO4)/dt
rates(3) = r_molar * MM_CaSO4 ! d(CaSO4)/dt
rates(4) = r_molar * MM_CO2 ! d(CO2)/dt
end subroutine get_derivatives
end program caco3_h2so4_rk4_passivationThis one models limestone with impurities, as long as you are not over about 150 lines it is ok, once past that it is a world of hurt.
Link Copied
9 Replies
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
AI makes stupid mistakes such as
integer k
REAL K(30)That one occurred in 2 of the three programs.
When you point out the error politely, it says, well it came in the supplied code or it should have known or some such excuse.
GOOGLE AI also suggests GFORTRAN.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
When converting python code into Fortran with ai, I need to pay attention that ifvthere are variables that has same name but with/without capitals.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
This sort of thing is not my area of expertise. If the program failed to compile, that would be another thing. I am not a fan of having AI write code for you - it tends to look good but be overly complex and have numerous flaws. Often it will hallucinate functions and libraries that don't exist.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Steve:
I agree with you, it proved to be a circular task, and yes numerous flaws.
But I was just interested to understand what it produced. It is fixable, but the cost would be as high as just developing existing code or writing it from scratch. Anyone who wastes corporate funds on this is not going to get a good return.
That is why I miss the two amigos so much.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
do i = 1, 2
call assemble_sparse_triplets(framework(i), triplet_count, trip_rows, trip_cols, trip_vals)
! 2. Compress triplets onto a temporary matrix to eliminate duplicate index slots
do i = 1, triplet_countcheck line 236-242.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Yes, it is a start to a 3D plastic program, but it has a long way to go, and I do not have the spare time or a job to test it on.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
module frame_types
implicit none
integer, parameter :: dp = selected_real_kind(15, 300)
type :: Element3D
integer :: node_i, node_j
real(kind=dp) :: E, G, A, J, Iy, Iz ! Material & Section properties
real(kind=dp) :: Mp_x, Mp_y, Mp_z ! Plastic moment capacities
real(kind=dp) :: Pp ! Plastic axial capacity
logical :: hinge_i = .false. ! Yield status at Node I
logical :: hinge_j = .false. ! Yield status at Node J
end type Element3D
end module frame_types
program plastic_frame_3d
use frame_types
implicit none
! integer, parameter :: dp = selected_real_kind(15, 300)
integer :: num_nodes, num_elements, ndof, i, step
integer :: max_steps = 20
real(kind=dp) :: load_factor, d_lf = 0.05_dp
type(Element3D), allocatable :: elements(:)
real(kind=dp), allocatable :: coords(:, ! Size (3, num_nodes)
real(kind=dp), allocatable :: K_global(:, ! Size (ndof, ndof)
real(kind=dp), allocatable :: F_inc(:), F_tot(:) ! Force vectors
real(kind=dp), allocatable :: U_tot(:), dU(:) ! Displacement vectors
integer, allocatable :: ipiv(:) ! Solver pivot array
integer :: info
! --- 1. Initialization ---
num_nodes = 2
num_elements = 1
ndof = num_nodes * 6
load_factor = 0.0_dp
allocate(elements(num_elements))
allocate(coords(3, num_nodes))
allocate(K_global(ndof, ndof))
allocate(F_inc(ndof), F_tot(ndof), U_tot(ndof), dU(ndof), ipiv(ndof))
! Example Geometry & Boundary conditions (A cantilever beam along X-axis)
coords(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
coords(:, 2) = [5.0_dp, 0.0_dp, 0.0_dp]
! Basic cross-section and material assumptions
elements(1) = Element3D(node_i=1, node_j=2, E=2.0e11_dp, G=7.7e10_dp, &
A=0.01_dp, J=2.0e-5_dp, Iy=5.0e-5_dp, Iz=5.0e-5_dp, &
Mp_x=1.5e5_dp, Mp_y=2.5e5_dp, Mp_z=2.5e5_dp, Pp=2.5e6_dp)
F_inc = 0.0_dp
F_tot = 0.0_dp
U_tot = 0.0_dp
F_inc(12) = -10000.0_dp ! Apply incremental vertical load at tip (Node 2, Z-dir)
print '(A)', "Starting 3D Plastic Hinge Incremental Analysis..."
! --- 2. Incremental Load Loop ---
do step = 1, max_steps
load_factor = load_factor + d_lf
! Assemble Tangent Stiffness Matrix
call assemble_global_stiffness(elements, coords, num_nodes, num_elements, ndof, K_global)
! Enforce boundary conditions (Simple fixity at Node 1: DOFs 1 to 6)
do i = 1, 6
K_global(i, = 0.0_dp
K_global(:, i) = 0.0_dp
K_global(i, i) = 1.0_dp
F_inc(i) = 0.0_dp
end do
! Solve for displacement increments dU using LAPACK routine DGESV
dU = F_inc
call dgesv(ndof, 1, K_global, ndof, ipiv, dU, ndof, info)
if (info /= 0) then
print '(A, I2)', "Error or singularity encountered in LAPACK solver. Info code: ", info
exit
end if
! Update state
U_tot = U_tot + dU
F_tot = F_tot + F_inc
! Check elements for new plastic hinge formation
call evaluate_plastic_hinges(elements, coords, dU)
print '(A, I2, A, F5.2)', " Step ", step, " completed. Load Factor = ", load_factor
end do
contains
subroutine assemble_global_stiffness(elems, nodes, n_n, n_e, nd, K_g)
type(Element3D), intent(in) :: elems(:)
real(kind=dp), intent(in) :: nodes(:, :)
integer, intent(in) :: n_n, n_e, nd
real(kind=dp), intent(out) :: K_g(:, :)
real(kind=dp) :: K_local(12, 12), L, E, G, A, J, Iy, Iz
integer :: ei, i, ji, dofs(12)
K_g = 0.0_dp
do ei = 1, n_e
! Local aliases for scannability
L = sqrt(sum((nodes(:, elems(ei).node_j) - nodes(:, elems(ei).node_i))**2))
E = elems(ei).E
G = elems(ei).G
A = elems(ei).A
J = elems(ei).J
Iy = elems(ei).Iy
Iz = elems(ei).Iz
K_local = 0.0_dp
! 1. Axial Stiffness (X-axis)
K_local(1,1) = E*A/L; K_local(1,7) = -E*A/L
K_local(7,7) = E*A/L; K_local(7,1) = -E*A/L
! 2. Torsional Stiffness (Rx-axis)
K_local(4,4) = G*J/L; K_local(4,10) = -G*J/L
K_local(10,10) = G*J/L; K_local(10,4) = -G*J/L
! 3. Bending in X-Y Plane (bending about Iz, translations along Y)
K_local(2,2) = 12.0_dp*E*Iz / L**3
K_local(2,6) = 6.0_dp*E*Iz / L**2
K_local(2,8) = -12.0_dp*E*Iz / L**3
K_local(2,12) = 6.0_dp*E*Iz / L**2
K_local(6,2) = 6.0_dp*E*Iz / L**2
K_local(6,6) = 4.0_dp*E*Iz / L
K_local(6,8) = -6.0_dp*E*Iz / L**2
K_local(6,12) = 2.0_dp*E*Iz / L
K_local(8,2) = -12.0_dp*E*Iz / L**3
K_local(8,6) = -6.0_dp*E*Iz / L**2
K_local(8,8) = 12.0_dp*E*Iz / L**3
K_local(8,12) = -6.0_dp*E*Iz / L**2
K_local(12,2) = 6.0_dp*E*Iz / L**2
K_local(12,6) = 2.0_dp*E*Iz / L
K_local(12,8) = -6.0_dp*E*Iz / L**2
K_local(12,12) = 4.0_dp*E*Iz / L
! 4. Bending in X-Z Plane (bending about Iy, translations along Z)
K_local(3,3) = 12.0_dp*E*Iy / L**3
K_local(3,5) = -6.0_dp*E*Iy / L**2
K_local(3,9) = -12.0_dp*E*Iy / L**3
K_local(3,11) = -6.0_dp*E*Iy / L**2
K_local(5,3) = -6.0_dp*E*Iy / L**2
K_local(5,5) = 4.0_dp*E*Iy / L
K_local(5,9) = 6.0_dp*E*Iy / L**2
K_local(5,11) = 2.0_dp*E*Iy / L
K_local(9,3) = -12.0_dp*E*Iy / L**3
K_local(9,5) = 6.0_dp*E*Iy / L**2
K_local(9,9) = 12.0_dp*E*Iy / L**3
K_local(9,11) = 6.0_dp*E*Iy / L**2
K_local(11,3) = -6.0_dp*E*Iy / L**2
K_local(11,5) = 2.0_dp*E*Iy / L
K_local(11,9) = 6.0_dp*E*Iy / L**2
K_local(11,11) = 4.0_dp*E*Iy / L
! Apply Plasticity Hinge Modifier Reductions
if (elems(ei).hinge_i) then
K_local(4:6, = K_local(4:6, * 1.0e-5_dp
K_local(:, 4:6) = K_local(:, 4:6) * 1.0e-5_dp
end if
if (elems(ei).hinge_j) then
K_local(10:12, = K_local(10:12, * 1.0e-5_dp
K_local(:, 10:12) = K_local(:, 10:12) * 1.0e-5_dp
end if
! Map element local DOFs to Global DOFs
dofs(1:6) = [( (elems(ei).node_i - 1) * 6 + i, i = 1, 6 )]
dofs(7:12) = [( (elems(ei).node_j - 1) * 6 + i, i = 1, 6 )]
do i = 1, 12
do ji = 1, 12
K_g(dofs(i), dofs(ji)) = K_g(dofs(i), dofs(ji)) + K_local(i, ji)
end do
end do
end do
end subroutine assemble_global_stiffness
subroutine evaluate_plastic_hinges(elems, nodes, dU_vec)
type(Element3D), intent(inout) :: elems(:)
real(kind=dp), intent(in) :: nodes(:, :)
real(kind=dp), intent(in) :: dU_vec(:)
! Inside real frameworks, element end-forces P, Tx, My, Mz are calculated here.
! Example placeholders for internal interaction calculations:
real(kind=dp) :: P = 0.0_dp, Tx = 0.0_dp, My = 1.2e5_dp, Mz = 1.0e5_dp
real(kind=dp) :: yield_val
integer :: e
do e = 1, size(elems)
! AISC/plastic design interaction criterion wrapper
yield_val = (P/elems(e).Pp) + (Tx/elems(e).Mp_x) + (My/elems(e).Mp_y) + (Mz/elems(e).Mp_z)
if (yield_val >= 1.0_dp .and. .not. elems(e).hinge_j) then
elems(e).hinge_j = .true.
print '(A, I2, A)', ">>> PLASTIC HINGE FORMED at Node J of Element ", e, "!"
end if
end do
end subroutine evaluate_plastic_hinges
end program plastic_frame_3dThis took 3 goes to sort out the AI errors, it runs, not sure if correct I would need to check against Strand7 and work through the code, then fix to Pardiso, add eigen solver and Monte Carlo the whole things, simple. LOL And fix the stiffness matrix for Ekhande's work.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Line 115 to the end of the stiffness matrix needs to be fixed, all in all it is interesting.
- Mark as New
- Bookmark
- Subscribe
- Mute
- Subscribe to RSS Feed
- Permalink
- Report Inappropriate Content
Reply
Topic Options
- Subscribe to RSS Feed
- Mark Topic as New
- Mark Topic as Read
- Float this Topic for Current User
- Bookmark
- Subscribe
- Printer Friendly Page