#This script performs struture optimization by material removal
#
#Author: Luigi Giaccari
#Date: 22/12/2010

#PARAMETERS
max_it=100#maximum number of iterations
DELETE_PERCENTAGE=0.02#at each iteration delete a given % of  elements (1=100%)
stiff=1000#stifness of the soft springs, high stifness make calculs more stable but less accurate

#GLOBALS
nonu=0#node numbers as global initialized to zero
elnu=0#element numbers as global initialized to zero
numtetra=0#number of tetraedrons that compose themodel
deleted=False#array flag for deleted elements
#LOADING MODULES
import sys
sys.path.append( '/opt/aster/13.6/lib/aster/Utilitai/' )
from Utilitai import partition
#from partition import *
from Utilitai import Table
import numpy as N
#import numeric as N
import math

#FUNCTIONS


def GetVmisesArray(VmisesValue,NodeLabel):#return an array with von mises stresses where the memory position points to the label-1
	#This suppose that nodes are labeled consecutively, for example N1,N2,N3 will have memeory position 0,1,2

        Vmises=[None]*nonu
        n=len(VmisesValue)#loop trough all compute von mises
        print n,' Nodes has von mises value'
       	for i in range(n):     
             value=VmisesValue[i]
             string=NodeLabel[i]
             mem_id = int(string[1:])-1 #making N1=0
             #print string
             #print mem_id 
             Vmises[mem_id]=value

        return Vmises


def median(arr):#get the median value of a list
        maximum=-10e100
        minimum=10e100
	n=len(arr)
        c=0
	for i in range(n):
             if arr[i]!=None:
                  if arr[i]>maximum:
                        maximum=arr[i]
                  if arr[i]<minimum:
                        minimum=arr[i]
        

	median=(maximum-minimum)/2
	print 'median value is: ',median
	return median



def MarkUnTouchElem(group_no,connex):#mark as untouchable a series of elements with unouchable nodes 
#group_no is a nodal group with memory address and "n" is total number of node in the model
	
        UnTouchNode=[False]*nonu
        UnTouchElem=[False]*elnu
	for i in group_no:
		UnTouchNode[i]=True

        for i in range(elnu):#loops trough elemnts
            numnod=len(connex[i])
            if (numnod!=4 and numnod!=10):#only apply to tetra #WARNING THIS LINE PREVENTS THE ALGORITHM TO WORK WITH EXA AND MIXED MESH
                        continue  
            for j in connex[i]:#loop trough elments nodes
                  if UnTouchNode[j]:
                       UnTouchElem[i]=True#elementis untouchable  
                       break#go to nxt element
	return UnTouchElem


def GetMeanVmises(Vmises,Connex,elnu):
#Computes mean von mises stress for each element
       
	MeanVmises=[None]*elnu
	for i in range(elnu):#loop trough all elements
                 if deleted[i]:
                       continue#jump deleted elements
                 numnod=len(Connex[i])#number of nodes of the elements 
                 #print 'numnod: ',numnod
              
                 if (numnod!=4 and numnod!=10):#only apply to tetra #WARNING THIS LINE PREVENTS THE ALGORITHM TO WORK WITH EXA AND MIXED MESH
                        continue 
                 mean=0#reset values before a new element
                 c=0
		 for j in Connex[i]:#loop trough all  nodes of the elemnet
                              mean=mean+Vmises[j]
                              c=c+1
                 if mean!=None:
			mean=mean/c#compute mean value inside element 
			MeanVmises[i]=mean
                        #print 'Element has:',len(Connex[i]),'nodes'
 
        #print MeanVmises 
       
	return MeanVmises

def BuildNewModel(MeanVmises,UnTouchElem):#Selects elements belonging to new model
 
       sortVmises=[]#list with sorted vmises value
       NewModel=[]#new model elemnts list    

       for i in range(elnu):
            if (deleted[i]==False and UnTouchElem[i]==False and MeanVmises[i] is not None):
                 sortVmises.append(MeanVmises[i])#fill the list with VonMises of remaining not untouchable tetras
       sortVmises.sort()#sort values in ascending order 
       #print sortVmises

       
        
       mem_id=int(elnu*DELETE_PERCENTAGE)#memory id value for treshold 
       
       treshold=sortVmises[mem_id]
       
       for i in range(elnu):
              
           if (deleted[i]==True):
               continue #skip deleted elements        
           if (MeanVmises[i]>treshold) or (UnTouchElem[i]):#Add Untouchable elements and high stressed ones
               NewModel.append('M'+str(i+1))#turn 0 to M1
           else:
               deleted[i]=True#element not added are deleted

       #compute some data
       maxstress=max(sortVmises)
       meanstress=sum(sortVmises)/len(sortVmises)
       minstress=min(sortVmises)
       uc=0#utilization coefficient
       for i in sortVmises:
             uc=uc+i/maxstress
       uc=uc/len(sortVmises)    
       #print NewModel 
       print 'ITERATION REPORT' 
       print 'The old model has: ',len(sortVmises),' tetra' 
       print 'tetra that are going to be deleted: ',mem_id+1  
       print 'treshold values is: ',treshold  
       print 'The new model has: ',len(NewModel),' tetra'
       print 'Maximum Stress is: ',maxstress
       print 'Minimum Stress is: ',minstress
       print 'Mean Stress is: ',meanstress
       print 'Utilization Coefficient is ',uc
       return NewModel
  

#START ASTER
DEBUT(PAR_LOT='NON',);


# initialise model parameters
msh=[None]*max_it
mod=[None]*max_it
nmh=[None]*max_it
sft=[None]*max_it
mat=[None]*max_it
res=[None]*max_it
bc=[None]*max_it
tab=[None]*max_it

#Material properties
Steel=DEFI_MATERIAU(ELAS=_F(E=2.1e11,
                            NU=0.27,
                            RHO=7800.0,),);

#Read MED MESH File

msh[0]=LIRE_MAILLAGE(UNITE=20,
                      FORMAT='MED',
                      INFO_MED=1,);

msh[0]=DEFI_GROUP(reuse =msh[0],
                MAILLAGE=msh[0],
                CREA_GROUP_NO=_F(GROUP_MA='model',NOM='TOUT',) ,);

#Since all other meshes are build on the first one we can pick connectivity data only on the first one
mesh1 = partition.MAIL_PY()#   Creating empty new MAIL_PY "class istantiation"
mesh1.FromAster(msh[0])#    Reading mesh Data From Aster using object referencing to the   class function FromAster
#    Get Number of Nodes and Number of Element
nonu=mesh1.dime_maillage[0]
elnu=mesh1.dime_maillage[2]
print 'the mesh has: ',nonu,' nodes and ',elnu,' elements'
#help(mesh1)


NodeLabel = list(mesh1.correspondance_noeuds)
ElemLabel = list(mesh1.correspondance_mailles)

ncoor = mesh1.cn#   Get Node coordinates referenced from node sequence position

#      Data are elements/nodes sequence positions
NodeGroup = mesh1.gno#nodes groups
ElemGroup = mesh1.gma#maillage groups
Connex   = mesh1.co#Get element connection as class CONNEC a two  numpy arrays


UnTouchElem=MarkUnTouchElem(NodeGroup['untouch'],Connex )#mark untouchable elemnts some points as untouched

numtetra=len((ElemGroup['model']))#number of strating tetraedrons
print 'Number of Tetraedrons:',numtetra
deleted=[False]*elnu

print os.getcwd()

for i in range(max_it):#start loop
        #Assigns model 
        nmh[i]=CREA_MAILLAGE(MAILLAGE=msh[0],
                             CREA_POI1=_F(NOM_GROUP_MA='spElmt',
                                          GROUP_NO='TOUT',),);
	mod[i]=AFFE_MODELE(MAILLAGE=nmh[i],
		           AFFE=(_F(GROUP_MA=('model','pres',),#assign the mod to all 3D elements not deleted plus pressure tria elements
		                   PHENOMENE='MECANIQUE',
		                   MODELISATION='3D',),
                                  _F(GROUP_MA='spElmt',
  		                  PHENOMENE='MECANIQUE',
 		                  MODELISATION='DIS_T',),),);

        #soft springs prevents orphan elements to create singularity 
        sft[i]=AFFE_CARA_ELEM(MODELE=mod[i],
                                DISCRET=_F(CARA='K_T_D_N',
                                            GROUP_MA='spElmt',
                                            VALE=(stiff,stiff,stiff,),),);

 

	#Assign Material properties to Elements
	mat[i]=AFFE_MATERIAU(MAILLAGE=nmh[i],
		          MODELE=mod[i],
		          AFFE=_F(GROUP_MA='model',
		                  MATER=Steel,),);

	#Boundary conditions
	bc[i]=AFFE_CHAR_MECA(MODELE=mod[i],
		            DDL_IMPO=(_F(GROUP_NO='fixed',
		                        DX=0.0,
		                        DY=0.0,
		                        DZ=0.0,),),
		            PRES_REP=_F(GROUP_MA='pres',
		                        PRES=1000.0,),);
        

	#Finite Element resu[i]

	res[i]=MECA_STATIQUE(MODELE=mod[i],
                               CARA_ELEM=sft[i],
		               CHAM_MATER=mat[i],
		               EXCIT=_F(CHARGE=bc[i],),);

	#Compute VonMises Stress / Strain

	res[i]=CALC_CHAMP(reuse =res[i],
		           RESULTAT=res[i],
                           EXCIT=_F(CHARGE=bc[i],),
		           CONTRAINTE=('SIGM_ELNO','EFGE_ELNO','SIGM_NOEU',),
                           CRITERES=('SIEQ_ELNO','SIEQ_NOEU'),);

        #create tab of von mises calues at nodes
        tab[i]=POST_RELEVE_T(ACTION=_F(OPERATION='EXTRACTION',
                              INTITULE='VonMises',
                              RESULTAT=res[i],
                              NOM_CHAM='SIEQ_NOEU',
                              GROUP_NO='TOUT',
                              NOM_CMP='VMIS',),);
 
	#Write Results to MED file

	IMPR_RESU(FORMAT='MED',
		  UNITE=80,
		  RESU=(_F(MAILLAGE=nmh[i],
		           RESULTAT=res[i],
		           NOM_CHAM='DEPL',
		           NOM_CMP=('DX','DY','DZ',),
                           GROUP_MA='model',),
		        _F(RESULTAT=res[i],
		           NOM_CHAM=('SIEQ_NOEU',),
		           NOM_CMP='VMIS',
                           GROUP_MA='model',),),);
	
        

       

           
		
         
	#extracting values from tables
	tab0 =  tab[i].EXTR_TABLE()
	#print tab
	pytab=tab0.values()
	#print pytab    

	Vmises=GetVmisesArray(pytab['VMIS'],pytab['NOEUD'])#getting an array with von mises sresses


	MeanVmises=GetMeanVmises(Vmises,Connex,elnu)#get MeanVmises for each element
	 
	NewModel=BuildNewModel(MeanVmises,UnTouchElem)#Select the elements belongin to new model
        #create the new mesh
        if i<max_it-1:      
	
                #create the new mesh copying the current one             
              
                #we also need to redefine the TOUT nodal group for table extraction from resu
                # and most of all the model element group for model definition
                msh[0]=DEFI_GROUP(reuse =msh[0],
                MAILLAGE=msh[0],
                DETR_GROUP_MA=_F(NOM='model',),
                DETR_GROUP_NO=_F(NOM='TOUT',),
                CREA_GROUP_MA=_F(NOM='model',MAILLE=NewModel,),
                CREA_GROUP_NO=_F(GROUP_MA='model',NOM='TOUT',) ,);
               

		#destroy all old concepts (be careful to keep the order in which they are created)
		DETRUIRE(CONCEPT = _F (NOM = nmh[i],) ,INFO=1, );
		DETRUIRE(CONCEPT = _F (NOM = mod[i],) ,INFO=1, );
                DETRUIRE(CONCEPT = _F (NOM = sft[i],) ,INFO=1, ); 
		DETRUIRE(CONCEPT = _F (NOM = mat[i],) ,INFO=1, );
		DETRUIRE(CONCEPT = _F (NOM = bc[i],) ,INFO=1, );
		DETRUIRE(CONCEPT = _F (NOM = res[i],) ,INFO=1, );
		DETRUIRE(CONCEPT = _F (NOM = tab[i],) ,INFO=1, ); 

	 
                del Vmises
                del MeanVmises
                del NewModel
                del pytab

     
	print 'Iteration: ',i,' ended'

DETRUIRE(CONCEPT = _F (NOM = msh[0],) ,INFO=1, );

FIN(FORMAT_HDF='OUI',);