Here is the serial code:
PROGRAM bmsolve_optimized
USE DISLIN
IMPLICIT NONE
!Compiler options: -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal -blas -fstack-arrays-lblas -fcoarray=single
!gfortran -ldislin -I /usr/local/dislin/g95 -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal-blas -fstack-arrays -lblas -fcoarray=single /usr/local/dislin/g95/dislin.f90 question_2_optimized.f90 -o question_2_optimized_gfort
!Possible Parrallelized version
!gfortran -ldislin -I /usr/local/dislin/g95 -march=skylake -mtune=skylake -O3 -mfpmath=sse -funroll-loops -fexternal-blas -fstack-arrays -lblas -fcoarray=single /usr/local/dislin/g95/dislin.f90 question_2_optimized.f90 -o question_2_optimized_gfort -fopenmp
!Data type to store details of the beam cross-section
TYPE data
REAL, ALLOCATABLE :: disp(:),moment(:)
REAL :: maximum_disp,maximum_values_disp(2)
REAL :: pb,pc
REAL :: In,A,d1,d2,d3,d4
REAL :: fa,fb,fc
REAL :: maximum_normal_stress,maximum_moment,maximum_values_moment(2)
INTEGER :: id
END TYPE data
!Parameters
REAL, PARAMETER :: q=8.,disp_all=5,pi=4.0*ATAN(1.0),l=6000,E=2E5
INTEGER,PARAMETER :: lsize=6001,dmax=500
!Variables
INTEGER :: i,j,k,r,t,maximum_disp_loc
INTEGER(KIND=8)::total
REAL, ALLOCATABLE :: div(:),divl(:)
REAL :: t1,t2
TYPE(data) :: s
TYPE(data) :: best_minimum_disp,best_minimum_bends,best_arr(5),best
!File name
CHARACTER(LEN=*),PARAMETER :: plt_file="question_2_best_plot.gp"
!Starting CPU clock and initializing calculation counter
CALL CPU_TIME(t1)
total=0
!Allocating arrays
ALLOCATE(div(dmax))
ALLOCATE(divl(lsize))
ALLOCATE(s%disp(lsize))
ALLOCATE(s%moment(lsize))
!Output format
100 FORMAT (X,A,25X,A)
110 FORMAT (X,A,G23.8,3X,A)
120 FORMAT (X,A,G36.8,3X,A)
130 FORMAT (X,A,G25.8,2X,A)
140 FORMAT (X,A,G27.8,X,A)
150 FORMAT (X,A,G26.8,2X,A)
160 FORMAT (X,A,G30.8,X,A)
170 FORMAT (X,A,G26.8,2X,A)
180 FORMAT (X,A,G18.8,1X,A)
190 FORMAT (X,A,G20.8,X,A)
191 FORMAT (X,A,G20.8,X,A)
192 FORMAT (X,A,G33.8,X,A)
DO i=1,lsize
divl(i)=REAL(i)
END DO
DO i=1,dmax
div(i)=REAL(i)
END DO
s%pb=2000.
s%pc=4000.
CALL find_react(s%fa, s%fb, s%fc, s%pb, s%pc)
PRINT *, " "
PRINT 192, "Initial reaction at A: ",s%fa,"N"
PRINT 192, "Initial reaction at B: ",s%fb,"N"
PRINT 192, "Initial reaction at C: ",s%fc,"N"
PRINT *, " "
s%id=1
s%d1=div(80)
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
CALL find_disp_moment(s%pb, s%pc, s%fb, s%fc, s%In, s%disp, s%moment)
CALL find_maximum_disp_moment(s%disp, s%moment, s%maximum_values_disp, s%maximum_values_moment, s%maximum_disp,&
s%maximum_moment)
IF (MAXLOC(s%maximum_values_disp,1)==1) THEN
maximum_disp_loc=MAXLOC(s%disp,1)
ELSE IF (MAXLOC(s%maximum_values_disp,1)==2) THEN
maximum_disp_loc=MINLOC(s%disp,1)
END IF
PRINT 192, "Location of maximum displacement: ",(maximum_disp_loc-1),"mm"
PRINT *, " "
!Solid circular cross-section
s%id=1
best_arr(1)%maximum_disp=0.
s%disp=0.
s%d2=0.
s%d3=0.
s%d4=0.
DO i=20,dmax
s%d1=div(i) !Radius
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
s%maximum_normal_stress=(s%maximum_moment*s%d1)/s%In
IF (s%maximum_disp>best_arr(1)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
best_arr(1)=s
END IF
total=total+1
END DO
CALL find_disp_moment(best_arr(1)%pb, best_arr(1)%pc, best_arr(1)%fb, best_arr(1)%fc, best_arr(1)%In, best_arr(1)%disp,&
best_arr(1)%moment)
PRINT 100, "SHAPE: ","Circle"
PRINT 110, "Moment of Inertia: ",best_arr(1)%In, "mm^4"
PRINT 120, "Area: ",best_arr(1)%A,"mm^2"
PRINT 130, "Maximum deflection: ",best_arr(1)%maximum_disp,"mm"
PRINT 140, "Global maxima: ",MAXVAL(best_arr(1)%disp,1),"mm"
PRINT 150, "Location: ",MAXLOC(best_arr(1)%disp,1)-1,"mm"
PRINT 160, "Global minima: ", MINVAL(best_arr(1)%disp,1),"mm"
PRINT 170, "Location: ",MINLOC(best_arr(1)%disp,1)-1,"mm"
PRINT 180, "Maximum normal stress: ",best_arr(1)%maximum_normal_stress,"Pa"
PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(1)%moment),"N-mm"
PRINT 191, "Global minima moment: ",MINVAL(best_arr(1)%moment),"N-mm"
PRINT 192, "Radius: ",best_arr(1)%d1,"mm"
PRINT *, " "
PRINT *, "Total calculations so far: ",total
PRINT *, " "
s%id=2
best_arr(2)%maximum_disp=0.
s%disp=0.
s%d3=0.
s%d4=0.
DO i=20,dmax-1
DO j=i+1,dmax
s%d1=div(i)
s%d2=div(j)
IF ((s%d2-s%d1)<10.) THEN
CYCLE
END IF
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
s%maximum_normal_stress=(s%maximum_moment*s%d2)/s%In
IF (s%maximum_disp>best_arr(2)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
best_arr(2)=s
END IF
total=total+1
END DO
END DO
CALL find_disp_moment(best_arr(2)%pb, best_arr(2)%pc, best_arr(2)%fb, best_arr(2)%fc, best_arr(2)%In, best_arr(2)%disp,&
best_arr(2)%moment)
PRINT 100, "SHAPE: ","Hollow Circle"
PRINT 110, "Moment of Inertia: ",best_arr(2)%In, "mm^4"
PRINT 120, "Area: ",best_arr(2)%A,"mm^2"
PRINT 130, "Maximum deflection: ",best_arr(2)%maximum_disp,"mm"
PRINT 140, "Global maxima: ",MAXVAL(best_arr(2)%disp,1),"mm"
PRINT 150, "Location: ",MAXLOC(best_arr(2)%disp,1)-1,"mm"
PRINT 160, "Global minima: ", MINVAL(best_arr(2)%disp,1),"mm"
PRINT 170, "Location: ",MINLOC(best_arr(2)%disp,1)-1,"mm"
PRINT 180, "Maximum normal stress: ",best_arr(2)%maximum_normal_stress,"Pa"
PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(2)%moment),"N-mm"
PRINT 191, "Global minima moment: ",MINVAL(best_arr(2)%moment),"N-mm"
PRINT 192, "Inner Radius: ",best_arr(2)%d1,"mm"
PRINT 192, "Outer radius: ",best_arr(2)%d2,"mm"
PRINT *, " "
PRINT *, "Total calculations so far: ",total
PRINT *, " "
s%id=3
best_arr(3)%maximum_disp=0.
s%disp=0.
s%d3=0.
s%d4=0.
DO i=20,dmax-1
DO j=i+1,dmax
s%d1=div(i)
s%d2=div(j)
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
s%maximum_normal_stress=(s%maximum_moment*(s%d2/2.))/s%In
IF (s%maximum_disp>best_arr(3)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
best_arr(3)=s
END IF
total=total+1
END DO
END DO
CALL find_disp_moment(best_arr(3)%pb, best_arr(3)%pc, best_arr(3)%fb, best_arr(3)%fc,best_arr(3)%In, best_arr(3)%disp,&
best_arr(3)%moment)
PRINT 100, "SHAPE: ","Solid Rectangle"
PRINT 110, "Moment of Inertia: ",best_arr(3)%In, "mm^4"
PRINT 120, "Area: ",best_arr(3)%A,"mm^2"
PRINT 130, "Maximum deflection: ",best_arr(3)%maximum_disp,"mm"
PRINT 140, "Global maxima: ",MAXVAL(best_arr(3)%disp,1),"mm"
PRINT 150, "Location: ",MAXLOC(best_arr(3)%disp,1)-1,"mm"
PRINT 160, "Global minima: ", MINVAL(best_arr(3)%disp,1),"mm"
PRINT 170, "Location: ",MINLOC(best_arr(3)%disp,1)-1,"mm"
PRINT 180, "Maximum normal stress: ",best_arr(3)%maximum_normal_stress,"Pa"
PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(3)%moment),"N-mm"
PRINT 191, "Global minima moment: ",MINVAL(best_arr(3)%moment),"N-mm"
PRINT 192, "Breadth: ",best_arr(3)%d1,"mm"
PRINT 192, "Height: ",best_arr(3)%d2,"mm"
PRINT *, " "
PRINT *, "Total calculations so far: ",total
PRINT *, " "
s%id=4
best_arr(4)%maximum_disp=0.
s%disp=0.
DO i=100,dmax-1
DO j=100,dmax-1
DO k=i+1,dmax
DO r=j+1,dmax
s%d1=div(i)
s%d2=div(j)
s%d3=div(k)
s%d4=div(r)
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
s%maximum_normal_stress=(s%maximum_moment*(s%d4/2.))/s%In
IF (s%maximum_disp>best_arr(4)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
best_arr(4)=s
END IF
total=total+1
END DO
END DO
END DO
END DO
CALL find_disp_moment(best_arr(4)%pb, best_arr(4)%pc, best_arr(4)%fb, best_arr(4)%fc,best_arr(4)%In, best_arr(4)%disp,&
best_arr(4)%moment)
PRINT 100, "SHAPE: ","Hollow Rectangle"
PRINT 110, "Moment of Inertia: ",best_arr(4)%In, "mm^4"
PRINT 120, "Area: ",best_arr(4)%A,"mm^2"
PRINT 130, "Maximum deflection: ",best_arr(4)%maximum_disp,"mm"
PRINT 140, "Global maxima: ",MAXVAL(best_arr(4)%disp,1),"mm"
PRINT 150, "Location: ",MAXLOC(best_arr(4)%disp,1)-1,"mm"
PRINT 160, "Global minima: ", MINVAL(best_arr(4)%disp,1),"mm"
PRINT 170, "Location: ",MINLOC(best_arr(4)%disp,1)-1,"mm"
PRINT 180, "Maximum normal stress: ",best_arr(4)%maximum_normal_stress,"Pa"
PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(4)%moment),"N-mm"
PRINT 191, "Global minima moment: ",MINVAL(best_arr(4)%moment),"N-mm"
PRINT 192, "Inner Breadth: ",best_arr(4)%d1,"mm"
PRINT 192, "Inner Height: ",best_arr(4)%d2,"mm"
PRINT 192, "Outer Breadth: ",best_arr(4)%d3,"mm"
PRINT 192, "Outer Height: ",best_arr(4)%d4,"mm"
PRINT *, " "
PRINT *, "Total calculations so far: ",total
PRINT *, " "
s%id=5
best_arr(5)%maximum_disp=0.
s%disp=0.
DO i=10,dmax
DO j=10,dmax
DO k=10,dmax
DO r=10,dmax
s%d1=div(i)
s%d2=div(j)
s%d3=div(k)
s%d4=div(r)
IF (s%d3>s%d1 .OR. (s%d2 .NE. s%d3)) THEN
CYCLE
END IF
CALL find_area_I(s%id, s%d1, s%d2, s%d3, s%d4, s%A, s%In)
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
s%maximum_normal_stress=(s%maximum_moment*((s%d4+(2*s%d2))/2.))/s%In
IF (s%maximum_disp>best_arr(5)%maximum_disp .AND. s%maximum_disp<disp_all) THEN
best_arr(5)=s
END IF
total=total+1
END DO
END DO
END DO
END DO
CALL find_disp_moment(best_arr(5)%pb, best_arr(5)%pc, best_arr(5)%fb, best_arr(5)%fc,best_arr(5)%In, best_arr(5)%disp,&
best_arr(5)%moment)
PRINT 100, "SHAPE: ","I-beam"
PRINT 110, "Moment of Inertia: ",best_arr(5)%In, "mm^4"
PRINT 120, "Area: ",best_arr(5)%A,"mm^2"
PRINT 130, "Maximum deflection: ",best_arr(5)%maximum_disp,"mm"
PRINT 140, "Global maxima: ",MAXVAL(best_arr(5)%disp,1),"mm"
PRINT 150, "Location: ",MAXLOC(best_arr(5)%disp,1)-1,"mm"
PRINT 160, "Global minima: ", MINVAL(best_arr(5)%disp,1),"mm"
PRINT 170, "Location: ",MINLOC(best_arr(5)%disp,1)-1,"mm"
PRINT 180, "Maximum normal stress: ",best_arr(5)%maximum_normal_stress,"Pa"
PRINT 190, "Global maxima moment: ",MAXVAL(best_arr(5)%moment),"N-mm"
PRINT 191, "Global minima moment: ",MINVAL(best_arr(5)%moment),"N-mm"
PRINT 192, "Flange Breadth: ",best_arr(5)%d1,"mm"
PRINT 192, "Flange Height: ",best_arr(5)%d2,"mm"
PRINT 192, "Web Breadth: ",best_arr(5)%d3,"mm"
PRINT 192, "Web Height: ",best_arr(5)%d4,"mm"
PRINT *, " "
PRINT *, "Total calculations so far: ",total
PRINT *, " "
best%maximum_disp=0.
DO i=1,5
IF (best_arr(i)%maximum_disp>best%maximum_disp) THEN
best=best_arr(i)
ELSE IF (best_arr(i)%maximum_disp == best%maximum_disp .AND. best_arr(i)%maximum_normal_stress<best%maximum_normal_stress) THEN
best=best_arr(i)
END IF
END DO
IF (best%id==1) THEN
PRINT *, "The best shape is the solid circle shape"
ELSE IF (best%id==2) THEN
PRINT *, "The best shape is the hollow circle shape"
ELSE IF (best%id==3) THEN
PRINT *, "The best shape is the solid rectangle shape"
ELSE IF (best%id==4) THEN
PRINT *, "The best shape is the hollow rectangle shape"
ELSE IF (best%id==5) THEN
PRINT *, "The best shape is the I-beam shape"
END IF
PRINT *, " "
!Assigns the best shape to the dummy variable
s=best
best_minimum_disp=best
best_minimum_bends=best
!Finding the lowest maximum displacement by finding best place for supports
DO i=2,lsize-1
DO j=i+1,lsize
s%pb=div(i)
s%pc=divl(j)
!Finding the reaction forces
CALL find_react(s%fa, s%fb, s%fc, s%pb, s%pc)
IF (s%fb .NE. s%fb .OR. s%fc .NE. s%fc) THEN
CYCLE
END IF
CALL find_disp_moment(s%pb, s%pc, s%fb, s%fc, s%In, s%disp, s%moment)
CALL find_maximum_disp_moment(s%disp, s%moment, s%maximum_values_disp, s%maximum_values_moment,&
s%maximum_disp, s%maximum_moment)
IF (MAXLOC(s%maximum_values_disp,1)==1) THEN
maximum_disp_loc=MAXLOC(s%disp,1)
ELSE IF (MAXLOC(s%maximum_values_disp,1)==2) THEN
maximum_disp_loc=MINLOC(s%disp,1)
END IF
s%maximum_disp=ABS(find_disp_point(s%pb, s%pc, s%fb, s%fc, s%In, maximum_disp_loc))
!Finding the maximum normal stress
IF (best%id==1) THEN
s%maximum_normal_stress=(s%maximum_moment*s%d1)/s%In
ELSE IF (best%id==2) THEN
s%maximum_normal_stress=(s%maximum_moment*s%d2)/s%In
ELSE IF (best%id==3) THEN
s%maximum_normal_stress=(s%maximum_moment*(s%d2/2.))/s%In
ELSE IF (best%id==4) THEN
s%maximum_normal_stress=(s%maximum_moment*(s%d4/2.))/s%In
ELSE IF (best%id==5) THEN
s%maximum_normal_stress=(s%maximum_moment*((s%d4+(2.*s%d2))/2.))/s%In
END IF
IF (s%maximum_disp<best_minimum_disp%maximum_disp) THEN
best_minimum_disp=s
END IF
IF (s%maximum_normal_stress<best_minimum_bends%maximum_normal_stress) THEN
best_minimum_bends=s
END IF
total=total+1
END DO
END DO
PRINT *, " "
PRINT *, "The best positions with the least displacment are: "
PRINT 192, "Support B at: ",best_minimum_disp%pb,"mm"
PRINT 192, "With reaction force: ",best_minimum_disp%fb,"N"
PRINT 192, "Support C at: ",best_minimum_disp%pc,"mm"
PRINT 192, "With reaction force: ",best_minimum_disp%fc,"N"
PRINT 192, "Maximum displacement: ",best_minimum_disp%maximum_disp,"mm"
PRINT 192, "Maximum normal stress: ",best_minimum_disp%maximum_normal_stress,"Pa"
PRINT *, " "
PRINT *, " "
PRINT *, "The best positions with the least maximum stress are: "
PRINT 192, "Support B at: ",best_minimum_bends%pb,"mm"
PRINT 192, "With reaction force: ",best_minimum_bends%fb,"N"
PRINT 192, "Support C at: ",best_minimum_bends%pc,"mm"
PRINT 192, "With reaction force: ",best_minimum_bends%fc,"N"
PRINT 192, "Maximum displacement: ",best_minimum_bends%maximum_disp,"mm"
PRINT 192, "Maximum normal stress: ",best_minimum_bends%maximum_normal_stress,"Pa"
PRINT *, " "
PRINT *, " "
PRINT *, "Total number of calculations are: ",total
PRINT *, " "
CALL CPU_TIME(t2)
PRINT 192, "Elapsed time: ",t2-t1,"Secs"
PRINT *, " "
PRINT *, "Printing graphs...."
CALL plot_with_dislin(divl, best%disp, lsize, 1,best%pb,best%pc)
CALL plot_with_dislin(divl, best_minimum_disp%disp, lsize, 2,best_minimum_disp%pb,best_minimum_disp%pc)
CALL plot_with_dislin(divl, best_minimum_bends%disp, lsize, 3,best_minimum_bends%pb,best_minimum_bends%pc)
PRINT *, "Done!"
DEALLOCATE(div); DEALLOCATE(divl); DEALLOCATE(s%disp); DEALLOCATE(s%moment)
CONTAINS
SUBROUTINE find_react(fa,fb,fc,pb,pc)
IMPLICIT NONE
!Data from main program
REAL,INTENT(IN) :: pb,pc
REAL,INTENT(OUT) :: fa,fb,fc
!Local variables
REAL :: detA
REAL :: A(2,2),A_in(2,2),X(2,1),B(2,1)
A(1,1)=((pb**2)*((3*pb)-pb))/6.
A(1,2)=((pb**2)*((3*pc)-pb))/6.
A(2,1)=((pb**2)*((3*pc)-pb))/6.
A(2,2)=((pc**2)*((3*pc)-pc))/6.
B(1,1)=(q*(pb**2)*((6*(l**2))-(4*l*pb)+(pb**2)))/24.
B(2,1)=(q*(pc**2)*((6*(l**2))-(4*l*pc)+(pc**2)))/24.
detA=1/((A(1,1)*A(2,2))-(A(1,2)*A(2,1)))
A_in(1,1)=+detA*A(2,2)
A_in(2,1)=-detA*A(2,1)
A_in(1,2)=-detA*A(1,2)
A_in(2,2)=+detA*A(1,1)
X(1,1)=(B(1,1)*A_in(1,1))+(A_in(1,2)*B(2,1))
X(2,1)=(B(1,1)*A_in(2,1))+(A_in(2,2)*B(2,1))
fb=X(1,1)
fc=X(2,1)
fa=(q*l)-fb-fc
END SUBROUTINE find_react
SUBROUTINE find_area_I(id,d1,d2,d3,d4,area,In)
IMPLICIT NONE
INTEGER, INTENT(IN) :: id
REAL, INTENT(IN) :: d1,d2,d3,d4
REAL, INTENT(OUT) :: area,In
!Perform checks and calculate area and I
IF (id == 1) THEN
area=pi*(d1**2)
In=(pi/4.)*(d1**4)
ELSE IF (id == 2) THEN
area=pi*((d2**2)-(d1**2))
In=(pi/4.)*((d2**4)-(d1**4))
ELSE IF (id == 3) THEN
area=d1*d2
In=(1./12.)*d1*(d2**3)
ELSE IF (id == 4) THEN
area=(d3*d4) - (d1*d2)
In=((1./12.)*d3*(d4**3))-((1./12.)*d1*(d2**3))
ELSE IF (id==5) THEN
area=(2*d1*d2)+(d3*d4)
In=(2*(((1/12.)*d1*(d2**3))+((d1*d2)*((d4/2)+d2)**2))) + ((1/12.)*d3*(d4**3))
ELSE
PRINT *, "INVALID SHAPE !!"
END IF
END SUBROUTINE find_area_I
SUBROUTINE find_disp_moment(pb,pc,fb,fc,In,disp,moment)
IMPLICIT NONE
REAL, INTENT(IN) :: pb,pc,fb,fc,In
REAL, INTENT(OUT) :: disp(:),moment(:)
disp=0.
moment=0.
!Load the displacment due to support B
disp(2:INT(pb))=disp(2:INT(pb))+(((fb*(divl(2:INT(pb))**2))*((3*pb)-divl(2:INT(pb))))/(6.*E*In))
disp(INT(pb)+2:lsize)=disp(INT(pb)+2:lsize)+(((fb*(pb**2))*((3*divl(INT(pb)+2:lsize))-pb))/(6.*E*In))
!Moment due to the support B
moment(1:INT(pb))=(moment(1:INT(pb)))+(fb*(pb-divl(1:INT(pb))))
!Load the displacement due to support C
disp(2:INT(pc))=disp(2:INT(pc))+(((fc*(divl(2:INT(pc))**2))*((3*pc)-divl(2:INT(pc))))/(6.*E*In))
disp(INT(pc)+2:lsize)=disp(INT(pc)+2:lsize)+(((fc*(pc**2))*((3*divl(INT(pc)+2:lsize))-pc))/(6.*E*In))
!Moment due to the support C
moment(1:INT(pc))=(moment(1:INT(pc)))+(fc*(pc-divl))
!Load the displacment due to the distributed load
disp=disp-((q*(divl**2))*((divl**2)+(6.*(l**2))-(4.*l*divl)))/(24.*E*In)
!Moment due to the distributed load
moment=moment-(q*((l-divl)**2))/2.
!Zeroing out the disp at the supports
disp(1)=0
disp(INT(pb)+1)=0
disp(INT(pc)+1)=0
moment(1)=0.
END SUBROUTINE find_disp_moment
SUBROUTINE find_maximum_disp_moment(disp,moment,maximum_values_disp,maximum_values_moment,maximum_disp,maximum_moment)
IMPLICIT NONE
REAL, INTENT(IN) :: disp(:),moment(:)
REAL, INTENT(OUT) :: maximum_values_disp(:),maximum_values_moment(:),maximum_disp,maximum_moment
!Finding max displacement
maximum_values_disp(1)=ABS(MAXVAL(disp,1))
maximum_values_disp(2)=ABS(MINVAL(disp,1))
maximum_disp=MAXVAL(maximum_values_disp,1)
!Finding the max moment
maximum_values_moment(1)=ABS(MAXVAL(moment,1))
maximum_values_moment(2)=ABS(MINVAL(moment,1))
maximum_moment=MAXVAL(maximum_values_moment,1)
END SUBROUTINE find_maximum_disp_moment
!Pure function to calculate the displacement at a single point
PURE REAL FUNCTION find_disp_point(pb,pc,fb,fc,In,point) RESULT(disp_point)
IMPLICIT NONE
REAL,INTENT(IN) :: pb,pc,fb,fc,In
INTEGER,INTENT(IN) :: point
disp_point=0.
IF (point<=INT(pb)) THEN
disp_point=disp_point+((fb*(divl(point)**2))*((3*pb)-divl(point)))/(6*E*In)
ELSE IF (point>INT(pb)) THEN
disp_point=disp_point+((fb*(pb**2))*((3*divl(point))-pb))/(6*E*In)
END IF
IF (point<=INT(pc)) THEN
disp_point=disp_point+((fc*(divl(point)**2))*((3*pc)-divl(point)))/(6*E*In)
ELSE IF (point> INT(pc)) THEN
disp_point=disp_point+((fc*(pc**2))*((3*divl(point))-pc))/(6*E*In)
END IF
disp_point=disp_point-((q*(divl(point)**2))*((divl(point)**2)+(6*(l**2))-(4*l*divl(point))))/(24.*E*In)
END FUNCTION find_disp_point
SUBROUTINE plot_with_dislin(x,y,n,check,p_b,p_c)
IMPLICIT NONE
INTEGER,PARAMETER :: wp=selected_real_kind(15)
REAL(KIND(wp)),ALLOCATABLE :: ox(:),oy(:),bx(:),cx(:),py(:)
INTEGER,INTENT(IN) :: n, check
REAL,INTENT(IN) :: p_b,p_c
REAL(KIND(wp)),INTENT(IN) :: x(:),y(:)
ALLOCATE(ox(n));ALLOCATE(oy(n));ALLOCATE(bx(3));ALLOCATE(cx(3));ALLOCATE(py(3))
ox=0.;oy=0.
bx=p_b;cx=p_c
py=[MINVAL(y,1),0.,MAXVAL(y,1)]
! Initialisation
CALL metafl("PDF") ! Set the file format
CALL setpag("DA3L") ! Set the page type
CALL scrmod('REVERS') ! Set the background
CALL disini() ! Initialize
! Set font
CALL triplx
! Labels
CALL name("Length (mm)","X") ! Label for the x axis
CALL name("Deflection (mm)","Y") ! Label for the y axis
IF (check .EQ. 1) THEN
CALL titlin("Deflection curve for best shape",1)
ELSE IF (check .EQ. 2) THEN
CALL titlin("Deflection curve for positions with least displacement",1)
ELSE IF (check .EQ. 3) THEN
CALL titlin("Deflection curve for positions with least stress",1)
END IF
! Axis set up
CALL labels('EXP','XYZ')
CALL graf(x(1)-1,REAL(n)+1,x(1)-1,1000.,MINVAL(y,1),MAXVAL(y,1)+0.01,MINVAL(y,1),((ABS(MAXVAL(y,1))+0.01)+ABS(MINVAL(y,1)))/10.)
!Create the title
CALL title
! Plot the graph
CALL incmrk(0)
CALL lintyp(2)
CALL curve(x,y,n)
CALL incmrk(0)
CALL lintyp(3)
CALL curve(x,oy,n)
CALL incmrk(0)
CALL lintyp(3)
CALL curve(bx,py,3)
CALL incmrk(0)
CALL lintyp(3)
CALL curve(cx,py,3)
!End dislin
CALL disfin
DEALLOCATE(ox);DEALLOCATE(oy);DEALLOCATE(bx);DEALLOCATE(cx);DEALLOCATE(py)
END SUBROUTINE plot_with_dislin
END PROGRAM bmsolve_optimized