!===============================================================================
! PROGRAM:          MMSM 1 Lab No 09
! PROJEKT:          REMODELING: Iterative adaptation of local material
!                   properties to local mechanical environment.
! TASK:             2D model of a plate with hole
!                   with initially constant material properties undergoes an
!                   iterative change of element material properties similar
!                   to bone remodeling or growth of wood.  The changing rate
!                   is modeled dependent on the local strain energy density.
! AUTHOR:           Simon, Wieland
! LAST CHANGES:     13.01.2012
!
! COMPUTER:
! OPERATING SYSTEM:
! ANSYS Version:    License: Teaching, 13.0
!
! CALL:             /input,Ue09-Remod.inp
!
!===============================================================================
 
!=== Beginning
finish              ! Stops all modules (preprocessor, solution, postprocessor)
/clear              ! Start a new analysis, delete the old database (file.db)
 
/PLOPTS, INFO, AUTO ! Good old contour labeling style
 
 
!===============================================================================
!  Parameters (lenght in mm, forces in N => stiffnesses and stresses in N/mm)
!===============================================================================
 
!--- Zur Programmsteuerung
nIter    =     10        ! Anzahl der Iterationen zum Materialumbau
 
!--- Zur Geometrie
Length   =     500.0        ! Length (x) of plate
Height   =     500.0        ! Height (y) of plate
Thickn   =      20.0        ! Thickness of plate
Diameter =     400.0        ! Diameter of the hole
x_hole   =     250.0        ! Verticale position of the holes center (from left)
y_hole   =     250.0        ! Horizontale position of the holes center (from bottom)
 
Netzfeinheit =     1        ! Globale Netzfeinheit 1(fein)..10(grob)
 
!--- Zur Last
Force    =  -100000.0       ! Load (Force, vertical)
 
!--- Zum initialen Werkstoff
Young_ini =  210000.0       ! Initial Young's modulus of the plates material (steel)
Poiss_ini =       0.3       ! Initial Poisson's ratio of the plates material (steel)
 
!--- Zur Remodelling-Funktion
SollWert  =  0.5            ! Sollwert fuer strain energy desnsity 
Epsilon   =  SollWert/10    ! Schaltbreite, (threshold)
FaktorAbb =  1000000            ! Steigung der Umbaufunktion Abbau
FaktorAuf =  FaktorAbb   ! Steigung der Umbaufunktion Aufbau
Limit_U  =  SollWert - Epsilon    ! Untere Schaltgrenze
Limit_O  =  SollWert + Epsilon    ! Obere Schaltgrenze
EMod_min =  Young_ini/100    ! Minimal zugelassene Werkstoffsteifigkeit
EMod_max =  Young_ini*2    ! Minimal zugelassene Werkstoffsteifigkeit






 
!===============================================================================
! A. Preprocessor (Setting up the model)
!===============================================================================
/prep7                       ! Switch to the preprocessor modul
 
 
!=== Build the geometry by means of bottom-up method
 
!--- Create plate
k,1,    0.0,    0.0, 0.0     ! Define keypoint via its coordinates x, y, z
k,2, Length,    0.0, 0.0     ! Define keypoint via its coordinates x, y, z
k,3, Length, Height, 0.0     ! Define keypoint via its coordinates x, y, z
k,4,    0.0, Height, 0.0     ! Define keypoint via its coordinates x, y, z
lstr,1,2                     ! Connect keypoints 1 and 2 to craete line 1
lstr,2,3                     ! Connect keypoints to craete a line
lstr,3,4                     ! Connect keypoints to craete a line
lstr,4,1                     ! Connect keypoints to craete a line
al,all                       ! Create area by all selected lines (1, 2, 3, 4)
 
!--- Create circular area
lsel,none                    ! Unselect all previous defined lines
k, 5, x_hole, y_hole, 0.0    ! Create keypoint at circle center
circle, 5, Diameter/2        ! Create circle lines
al,all                       ! Create area by all selected lines
allsel                       ! Select everything
 
!--- Loch ausschneiden (Boolsche Operation)
asba,1,2,,delete,delete      ! cut out hole (2) from plate (1) and get plate with hole (3)
/pnum,area,1                 ! Switch numbering of lines on
aplot                        ! Plot areas


!=== Meshing
!et,1,plane42                 ! Define the local element type 1 as a PLANE42 element
et,1,plane182                 ! Define the local element type 1 as a PLANE42 element
keyopt,1,3,3                  ! Sets key option 3 (of elem type 1) to 3 (plane stress with thickn)
smrtsize, Netzfeinheit        ! Globale Netzfeinheit 1(fein)..10(grob)
amesh, all                    ! Meshing areas with this number
eplot                         ! Plot elements
*GET, nElem, ELEM, 0, COUNT   ! Get number of elements


!=== Initial definition of material properties, one for each of the FE

!ToDo: Schleife ueber alle Elemente 

   mp,ex,1,Young_ini              ! E-Modul zuweisen
   mp,prxy,1,Poiss_ini            ! Querkontraktionszahl zuweisen

!ToDo  emodif,...                     ! Element el Mat.typ el zuweisen
 
!=== Define real constants
r,1,Thickn                   ! Define thickness of the plate into real set number 1

!=== Apply Load and Boundary Conditions
!--- Displacement boundary conditions
nsel,s,loc,x,0.0             ! Select all nodes at x = 0 (left side)
d,all,ux,0.0                 ! Set displacement ux to zero for selected nodes
d,all,uy,0.0                 ! Set displacement uy to zero for selected nodes
allsel                       ! Select all entities
!--- Applying the load as a single force
nn = node(Length,Height,0.0) ! Get node number at location x, y, z (bottom right corner)
f,nn,fy,Force                ! Apply load at node with number nn
 
!===============================================================================
!=== Arrays (fuer elementbez. Daten) anlegen
*dim,AR_IstWert,array,nElem     ! Ein Feld (Vektor) anlegen
*dim,AR_EMod,array,nElem        ! Ein Feld (Vektor) anlegen
*dim,AR_DeltaE,array,nElem      ! Ein Feld (Vektor) anlegen



 
!===============================================================================
! B. START REMODELING LOOP
!===============================================================================
iter = 1
!ToDo Schleife ueber Umbauschritte
   !============================================================================
   !====== Solution
   /solu                              ! Switch to the solution module
   antype,0                           ! Select the static analysis type
   nlgeom,on
   solve                              ! Solve current load step
   !============================================================================
   !====== Postprocessor
   /post1                             ! Switch to the postprocessor module
   
   !============================================================================
   !=== SED Auslesen
   etable,ET_SED,SEND,ELASTIC         		! Elementweises Erg. in der "E-Table" speichern
   *DO,el,1,nElem                     		! Schleife ueber alle Elemente
      *GET,AR_IstWert(el),ELEM,el,etab,ET_SED   ! ... und dann Element-Erg auslesen
   *ENDDO
 
   !============================================================================
   !=== Zwischenergebnisse ausgeben
   /title,'Strain-Energy-Density (SED) Schritt = %iter%'
   /dscale,1,1.0
   /contour,,,0,,2*Limit_O
   plesol, SEND, ELASTIC
   
   !============================================================================
   !=== Pause einlegen und evtl. abbrechen
   !*ask,aa,'= Continue! 0 = Stop',1
   !*IF,aa,eq,0,:ENDE
 
   !============================================================================
   !=== Umbaufunktion auswerten

   !ToDo ...
   !ToDo AR_DeltaE(el) = ...			! E-Modul-Aenderung berechnen
   !ToDo ...
 
   !============================================================================
   !=== Remodeling: Werkstoff aendern
   /prep7                			! Switch to Preprocor

   !ToDo ...
   !ToDo   *GET,mat_no,mat,el        		! Get material number of each elem
   !ToDo   *GET,AR_EMod(el),ex,mat_no        	! Get E modulus of that material number
   !ToDo ...
   !ToDo   AR_EMod(el) = ...                         ! E-Modul aendern
   !ToDo   mp,ex,el,AR_EMod(el)                 	! E-Modul zuweisen
   !ToDo ...


!===============================================================================
! END OF LOOP
!===============================================================================

 




!===============================================================================
! C. POSTPROCESSOR
!===============================================================================
/post1                    ! Switch to Postprocesser

!=============================================================================== 
!=== Plot des E-Moduls
!--- Element-Table anlegen; willk. mit epto,x; spaeter ueberschreiben
!etable,ET_emod,epto,x            ! zur Visualisierung der E-Module

!--- ElemTable mit E-Modeul fuellen
!*vput,AR_Emod(1),elem,1,etab,ET_emod,0    ! E-Modul

!--- Plot
!/contour,,,EMod_min,,EMod_max
!/title,'E-Modul-Verteilung nach Schritt = %nIter%'
!pletab, ET_emod, noavg

!===============================================================================  
!=== Pause einlegen und evtl. abbrechen
!*ask,aa,'= Continue! 0 = Stop',1
!*IF,aa,eq,0,:ENDE
 
!=============================================================================== 
!=== SED
/title,'Strain-Energy-Density (SED) Schritt = %nIter%'
/dscale,1,1.0
/contour,,,0,,2*Limit_O
plesol, SEND, ELASTIC


 

 
 
!===============================================================================
! END LABLE
!===============================================================================
:ENDE