Arlington: The Rap

FEM Laplace Problem

In der Zwischenzeit haben wir in unserer FEM Vorlesung die Elementmatrizen für das 2D Laplaceproblem hergeleitet.

Die Bilder sind das Ergebnis meines eigenen FE-Solvers das ich dann mit Exel geplottet habe.

Für alle die sich jetzt Fragen was man sieht- diese kurze Antwort:

Es handelt sich quasi um eine Platte (zum Beispiel Blech) die ringsherum isoliert ist bis auf zwei Stellen am Rand wo einmal die Temperatur bei 100C und einmal bei 0C festgehalten wird.

An zwei Knoten im Feld wird einmal ein positiver und einmal ein Negativer Wärmestrom eingeleitet ( das sind die spitzen Dellen).

Bob Dylan lebt!

Beinahe unglaublich aber wahr. 

Bob Dylan tourt noch immer durch die Lande. Am Samstag hat er in Washington halt gemacht und ein grandioses Konzert hingelegt.

Ich habe mich eigentlich ein bisschen über die  Location gewundert, denn so groß war die Arena  gar nicht. 

Führung durch das Capitol

Diesen Dienstag Morgen haben wir eine Führung durch das Capitol genommen. 

Erster Crash-Test

Neulich konnten wir unseren ersten Crash Test miterleben.

Es scheint irgendwie ein teures Spielzeug für große Kinder zu sein 

Die Vorbereitung kann mehrere Tage in Anspruch nehmen  und das Event ist spätestens nach einer Sekunde vorbei.  Danach kommt dann die Datenauswertung. In diesem Fall darf ein Mitstudent die Dummy Kinematik auswerten.

Hausaufgaben

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

Massen Demo in Washington

Am Samstag wollte ich auf der bis zur National Mall joggen. Dort bin ich "ganz zufällig" auf eine Massendemonstration gegen die Tea Party Bewegung gestoßen.

Selbst Cat Stevens hat wohl dort ein Ständchen geträllert von dem ich allerdings nicht viel gesehen habe da ich wohl einer der letzten zehn tausend war die gekommen waren um diese Spaßveranstaltung mitzuerleben.

Wer genaueres über die Veranstaltung wissen will ist unter Spiegel Online richtig. 

Halloween

Gestern war hier Halloween und ich muss sagen das hat nicht sehr viel mit Halloween in Deutschland zu tun hat.
Hier finden nämlich auch Kostümparties weit über das Kindesalter hinaus statt.
Einer meiner Mitbewohner (der ist auch schon über zwanzig) hat sogar noch an Häusern geklingelt um Süßigkeiten zu bekommen.
Ich jedenfalls habe mich dann als das Gecko von Gaiko (das ist ein Autoversicherer) verkleidet.
Ein paar Videoclips vom Gecko finden sich (mal wieder auf YouTube) hier und hier.

Shenandoah National Park

Am vergangenen Freitag waren wir (Steffen, Alex und ich, (die deutschen Austauschstudenten), die Schwester von Alex und zwei französische Austauschstudenten) gemeinsam im Shenandoah National Park. Das ist übrigens ungefähr der Ort wo Salle sein FSJ gemacht hat.

Im Prinzip ist es wie das Ermstal nur dass alles ein bisschen größer ist und unverbauter.

Wir sind ein paar zig Meilen auf dem Skyline Drive gefahren. Hier wandert man nämlich nur im Notfall, es ist ja viel bequemer wenn man sich die Landschaft aus dem Auto heraus ansieht.

Wir haben auch Geräusche gehört die wir einem Bären zuschreiben würden, gesehen haben wir ihn aber leider nicht. 

Im übrigen handelt auch der Country Roads Song von genau dieser Gegend. Hier nochmals zum nachhören: "....Blue rich mountain, Shenendoah River ...."

Besuch von meinem Vater

 

Letzte Woche habe ich Besuch von meinem Vater erhalten. Nachdem er geschäftlich in Detroit war hat er das Wochenende bei mir in Washington verbracht.

Gemeinsam waren wir in der Kunstgalerie und in einem Burger Restaurant in dem hin und wieder auch Obama isst. 

Übrigens das Bild im Hintergrund stammt von Archimbaldo und entstand im 16. Jahrhundert.

Wer mehr über seine faszinierenden Bilder wissen will ist hier richtig.