Weil nun der Eindruck entstehen könnte ich hätte zu viel Zeit übrig und müsste nichts lernen, veröffentliche ich nun hier meine aktuelle Hausaufgabe in FEM.
Program FEM ! Createt by Martin Werz
implicit none
Integer, parameter :: rk=selected_real_kind(14)
Integer:: Numnp ! Number of Nodal Points
integer:: Numel ! Number of Elements
integer:: NumBC ! Number of Boundary conditions
real :: Tens ! Tension
integer, allocatable, Dimension (:,:) :: NP ! Nodal Points of Element(2,Numel)
real, allocatable, Dimension(:) :: XCO ! X-Coor. of Nodalpoints (Numnp)
real(kind=rk), allocatable, Dimension(:,:) :: GlobStiff ! (Numel,Numel)
real(kind=rk), allocatable, Dimension(:,:) :: GlobStiffMod ! (Numel,Numel), Modified EQ.
real(kind=rk), allocatable, Dimension(:,:) :: InvStiff ! (Numel,Numel)
logical, allocatable, Dimension(:) :: BC ! Boundry Condition Force/ Displ.
real, allocatable, Dimension(:) :: Val_BC ! Value of Boundary condition
real, allocatable, Dimension(:) :: Val_BCMod! Boundary for the Modified EQ.
real, allocatable, Dimension(:) :: Displ , Force ! Solution Vectors.
character * 32 :: keyword !reading the line, searching for keyword.
integer :: len_to_key ! number of lines until the keyword apears.
integer :: io_stat
real:: act_XCO ! XCo which is working on
character *1 :: Kind_BC
real :: loc_Val_BC
integer::act_Node ! Node wich is working on.
integer :: act_Element ! Element wich is working on.
integer :: left_Node
integer :: right_Node
real :: Xlen
integer :: ii,jj,ll,kk,mm ! just for counting.
character *32 :: Input
character *32 :: FileIn, Fileout
character*1 :: contr
3001 continue
Print*,'Please type in the file name of the FEM input file'
read(*,'(A32)'), Input
FileIn=trim(input)
open(unit=1000, file=FileIn, status='old', action='read', iostat=io_stat) ! open the file
If(io_stat.ne.0) then
Print*, 'File not found'
goto 3001
endif
3002 continue
Print*,'Please type in the file name of the FEM output file'
read(*,'(A32)'), Input
Fileout=trim(input)
open(unit=2000, file=Fileout, status='new', iostat=io_stat) ! open the file
If(io_stat.ne.0) then
Print*, 'File exists already'
Print*, 'Please type in a name wich does not exist'
goto 3002
endif
Print*, 'When you type in 0 the system matrices are not shown'
Print*, 'Otherwise they are shown. Take 0 for large problems'
read*,contr
do while (index(keyword,'*TENSION').eq.0)
read(1000,'(A32)',iostat=io_stat) keyword
!Print*,keyword,'*'
enddo
read(1000,*) tens
print *,'tension=', tens !Printin tension
rewind(1000) ! reading procedure for the Nodes
Print*,'reading Nodes'
len_to_key=0
do while (index(keyword,'*NODE').eq.0)
read(1000,'(A32)',iostat=io_stat) keyword
len_to_key=len_to_key+1
!Print*,trim(keyword), len_to_key
enddo
numnp=0
do while(io_stat==0)
read(1000,*,iostat=io_stat) act_Node, act_XCO
numnp=numnp+1
!Print*,numnp,act_Node
enddo
numnp=numnp-1
print *, 'Number of Nodal Points =', numnp
allocate (xco(numnp))
rewind(1000)
do ii=1, len_to_key
read (1000,*)
enddo
do ii=1, numnp
read (1000,*) act_node, act_xco
xco(act_node)=act_xco
enddo
Print *,'The nodal Points with theier coordinates'
If(contr=='0') goto 5001
do ii=1, numnp
Print*, ii, xco(ii) !printing the list of read in Coordinates
enddo
5001 continue
rewind(1000) ! reading procedure for the Elements
Print*,'Reading Elements'
len_to_key=0
do while (index(keyword,'*ELEMENT').eq.0)
read(1000,'(A32)') keyword
len_to_key=len_to_key+1
enddo
numel=0
io_stat=0
do while(io_stat==0)
read(1000,*,iostat=io_stat) act_Element, left_node, right_node
!Print*,act_Element, left_node, right_node
numel=numel+1
enddo
numel=numel-1
Print*,'number of Elements=', numel
allocate (NP(2,numel))
rewind(1000) ! Going back in File
do ii=1, len_to_key
read (1000,*)
enddo
do ii=1, numel
read(1000,*,iostat=io_stat) act_element, left_node, right_node
!Print*,act_element, left_node, right_node, 'actual reading'
NP(1, act_element)=left_node
Np(2, act_element)=right_node
enddo
if(contr=='0') goto 5002
do ii=1, numel
Print *,'Element',ii,'has Nodal Points', NP(1,ii), 'and', NP(2,ii)
enddo
5002 continue
rewind(1000) ! reading procedure for the Boundary conditions
len_to_key=0
do while (index(keyword,'*BC').eq.0)
read(1000,'(A32)') keyword
len_to_key=len_to_key+1
enddo
numBC=0
do while(io_stat==0)
read(1000,*,iostat=io_stat) act_node, kind_BC, Loc_val_BC
!Print*, act_node, kind_BC, Loc_val_BC
numBC=numBC+1
enddo
numBC=numBC-1
Print *,'number of Boundary Conditions =', numBC
allocate (val_BC(numnp))
allocate (val_BCMod(numnp))
allocate(BC(numnp))
val_BC=0
BC=.false.
rewind(1000)
do ii=1, len_to_key
read (1000,*)
enddo
do ii=1, numBC ! Reading the Boundary conditions
read(1000,*) act_node, kind_BC, Loc_val_BC
!Print*,act_node, kind_BC, Loc_val_BC
val_BC(act_node)= Loc_val_BC
if (kind_BC=='D') then
Bc(act_node)=.true.
endif
enddo
if (contr=='0') goto 5010
do ii=1, numnp
if (BC(ii)) then
kind_bc='D'
else
Kind_bc='F'
endif
Print *,'Nodal point',ii,'has the BC ', kind_bc,' with the value', val_BC(ii)
enddo
5010 continue
!!!!!!!!!!!!!Here is the reading Part finsish.!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
allocate (GlobStiff(numnp,numnp))
allocate (GlobStiffMod(numnp,numnp))
allocate (invStiff(numnp,numnp))
allocate (Displ(numnp))
allocate (Force(numnp))
GlobStiff=0._rk
!!!!!!!!!
!Print*, BC
!!!!! creation of the local Stiffness matrix
Do act_Element=1, numel
xlen=xco(NP(2,act_Element))-xco(NP(1,act_Element)) !lenght of actual Element.
GlobStiff(NP(1,act_Element),NP(1,act_Element))=GlobStiff(NP(1,act_Element),NP(1,act_Element))+(Tens/xlen)
GlobStiff(NP(2,act_Element),NP(2,act_Element))=GlobStiff(NP(2,act_Element),NP(2,act_Element))+(Tens/xlen)
GlobStiff(NP(2,act_Element),NP(1,act_Element))=GlobStiff(NP(2,act_Element),NP(1,act_Element))-(Tens/xlen)
GlobStiff(NP(1,act_Element),NP(2,act_Element))=GlobStiff(NP(1,act_Element),NP(2,act_Element))-(Tens/xlen)
enddo
GlobStiffMod=GlobStiff ! Because GlobStiff is needed for missing Forces.
val_BCMod=Val_BC
Do ii=1,numnp
if (BC(ii)) then
if (abs(GlobStiff(ii,ii)).le.0.0001) then
Print*,'**************************************************************ERROR IN SYSTEM MATRIX'
Print*,'******************************************Look if all Nodes are conected by Elements'
Write(2000,*)'********************************************************ERROR IN SYSTEM MATRIX'
Write(2000,*)'************************************Look if all Nodes are conected by Elements'
endif
GlobStiffMod(ii,:)=0
GlobStiffMod(:,ii)=0
GlobStiffMod(ii,ii)=1
Do jj=1,numnp
Val_BCMod(jj)=Val_BCMod(jj)-val_BC(ii)*GlobStiff(jj,ii) ! First part of mod. the vektor.
enddo
endif
enddo ! GlobStiffMod is finished here.
Do ii=1, numnp ! second part of modifiing the vektor.
if(BC(ii))then!!!!!!!!!!!!!!!!!!!!!!!!!!! Not sure if this is correct
val_BCMod(ii)=val_BC(ii)
endif
enddo
!!!!!!!!!Now we need the inverse of the GlobStiffMod Matrix
InvStiff=0._rk
Do ii=1, numnp
InvStiff(ii,ii)=1
Enddo
!Print*, 'inverse Stiffness at beginning'
!Do ii=1, numnp ! zum Testen
!Print*, invStiff(ii,:)
!enddo
if (contr=='0') goto 5003
Print*, 'System Matrix:'
Do ii=1, numnp
Print*, GlobStiffMod(ii,:)
enddo
5003 Continue
!!!!!!!!!!!!!Here Starts the inversion of the System Matrix
Do ii=1, numnp
!Solving the problem with zero on the main diagonal
6000 continue
if (abs(globStiffmod(ii,ii)).le.0.0001) then
Print*,'zero on main diagonal'
Do kk=ii, numnp
if (abs(globStiffmod(kk,ii)).gt.0.01)then
invStiff(ii,:)=invStiff(ii,:)+invStiff(kk,:)
GlobStiffMod(ii,:)=GlobStiffMod(ii,:)+GlobStiffMod(kk,:)
goto 6000
endif
enddo
Print*,'kk=',kk
if(kk==numnp) then
Print*, 'untere schleife'!!!!!!!!Test
Do ll=1, ii-1
if (abs(globStiffmod(ll,ii)).gt.0.01)then
invStiff(ii,:)=invStiff(ii,:)+invStiff(ll,:)
GlobStiffMod(ii,:)=GlobStiffMod(ii,:)+GlobStiffMod(ll,:)
Do jj=1, ii-1
invStiff(ii,:)=invStiff(ii,:)-invStiff(jj,:)*GlobStiffMod(ii,jj)
GlobStiffMod(ii,:)=GlobStiffMod(ii,:)-GlobStiffMod(jj,:)*GlobStiffMod(ii,jj)
enddo
if (abs(globStiffmod(ii,ii)).le.0.001) then
Print*,'Matrix is not invertable'
Write(2000,*) 'Matrix is not invertable'
contr='1'
goto 7000
endif
goto 6000
endif
enddo
endif
endif
invStiff(ii,:)=invStiff(ii,:)/GlobStiffMod(ii,ii)
GlobStiffMod(ii,:)=GlobStiffMod(ii,:)/GlobStiffMod(ii,ii)
Do jj=ii+1, numnp
invStiff(jj,:)=invStiff(jj,:)-invStiff(ii,:)*GlobStiffMod(jj,ii)
GlobStiffMod(jj,:)=GlobStiffMod(jj,:)-GlobStiffMod(ii,:)*GlobStiffMod(jj,ii)
enddo
enddo
! Here is the Matrix in a Triangular form.
Do jj=numnp, 1, -1
!Print*,'jj=',jj
do ii=1, jj-1
!Print*,'ii=',ii
invStiff(ii,:)=invStiff(ii,:)-invStiff(jj,:)*GlobStiffMod(ii,jj)
GlobStiffMod(ii,:)=GlobStiffMod(ii,:)-GlobStiffMod(jj,:)*GlobStiffMod(ii,jj)
!Print*,'Inverse'
!Do ll=1, numnp ! zum Testen
!Print*, invStiff(ll,:)
!enddo
!Print*,'GlobStiffMod'
!Do ll=1, numnp
!Print*, GlobStiffMod(ll,:)
!enddo
enddo
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!End of Inversion Process
!Print*, 'Global Stiff Mod'
!Do ii=1, numnp
!Print*, GlobStiffMod(ii,:)
!enddo
if(contr=='0') goto 5004
Print*,'The Inverse System Matrix is'
Do ll=1, numnp ! zum Testen
Print*, invStiff(ll,:)
enddo
5004 continue
Displ=0
Do ii=1, numnp !Matrixproduct between modif.inverseStiffnes and modif. val_BC
do jj=1, numnp
Displ(ii)=Displ(ii)+invStiff(ii,jj)*Val_BCMod(jj)
enddo
enddo
if(contr=='0') goto 5005
Print*, 'Displacement Vector='
Do ii=1, numnp
Print*,'Node',ii,':',Displ(ii)
enddo
5005 continue
Force=0
Do ii=1, numnp !Matrixproduct between Stiffnes Matrix and Displacements
do jj=1, numnp !to get the missing Forces
Force(ii)=Force(ii)+GlobStiff(ii,jj)*Displ(jj)
enddo
enddo
if(contr=='0') goto 5006
Print*, 'Force Vector='
Do ii=1, numnp
Print*,'Node',ii,':',Force(ii)
enddo
5006 continue
!!!!!!!!!!!!!!!!!!!!!!!!!!This part is to create an output file
7000 continue
Write(2000,*) 'This is the sulution to the Inputfile', Filein
Write(2000,*) 'The Stiffnes Matrix is' !Write Stiffness Matrix
Write(2000,*)
if(contr=='0') goto 5007
Do ii=1, numnp
Do jj=1, numnp
Write(2000,3010,Advance='no'), GlobStiff(ii,jj)
3010 Format(F10.5,' ; ')
enddo
Write(2000,*)
enddo
Write(2000,*) 'The modified System Matrix is' !Write the System Matrix
Write(2000,*)
Do ii=1, numnp
Do jj=1, numnp
Write(2000,3011,Advance='no'), GlobStiffMod(ii,jj)
3011 Format(F10.5,' ; ')
enddo
Write(2000,*)
enddo
5007 continue
!!!!!!!!!!!!!!!!!!!!!!!!!!Writeing the Solution
Write(2000,*) '**********************************************************'
Write(2000,*) 'These are complete Information ablut the system:'
Write(2000,*) 'Node ; X-Ccoordinate ; Force ; Deflection'
Do ii=1, numnp
Write(2000,3012) ii, xco(ii),Force(ii),Displ(ii)
3012 Format(I8,' ; ' ,F15.5,' ; ', F15.5,' ; ', F15.5)
enddo
Print*,'Program run successfuly'
end program FEM
