#1 2010-12-22 11:38:00

Bracchesimo

Member

From: Italy

Registered: 2010-08-17

Posts: 100

Website

Structure Optimization with Code_Aster

I share this script that I used to perform a simple structure optimization by material removal.

Feedback and collaboration is appreciated.

Luigi

A dedicated post can be found here:
http://www.advancedmcode.org/structure- … aster.html

#This script performs struture optimization by material removal
#
#Author: Luigi Giaccari
#Date: 22/12/2010
 
#PARAMETERS
max_it=1#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/SALOME-MECA-2010.1-x86_64/aster/STA10.1/bibpyt/Utilitai' )
from Utilitai.partition import *
from partition import *
from Table import*
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
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,);
 
#Since all other meshes are build on the first one we can pick connectivity data only on the first one
mesh1 = 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 
 mod[i]=AFFE_MODELE(MAILLAGE=msh[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_NO='TOUT',
 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_NO='TOUT',
 VALE=(stiff,stiff,stiff,),),);
 
 
 
 #Assign Material properties to Elements
 mat[i]=AFFE_MATERIAU(MAILLAGE=msh[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_ELEM(reuse =res[i],
 RESULTAT=res[i],
 OPTION=('SIGM_ELNO_DEPL','EQUI_ELNO_SIGM',),);
 
 #create tab of von mises calues at nodes
 tab[i]=POST_RELEVE_T(ACTION=_F(OPERATION='EXTRACTION',
 INTITULE='VonMises',
 RESULTAT=res[i],
 NOM_CHAM='EQUI_ELNO_SIGM',
 GROUP_NO='TOUT',
 NOM_CMP='VMIS',),);
 
 #Write Results to MED file
 
 IMPR_RESU(MODELE=mod[i],
 FORMAT='MED',
 UNITE=80,
 RESU=(_F(MAILLAGE=msh[i],
 RESULTAT=res[i],
 NOM_CHAM='DEPL',
 NOM_CMP=('DX','DY','DZ',),),
 _F(RESULTAT=res[i],
 NOM_CHAM=('EQUI_ELNO_SIGM',),
 NOM_CMP='VMIS',),),);
 
 
 
 
 
 
 
 
 #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             
 msh[i+1]=CREA_MAILLAGE(MAILLAGE=msh[i],);
 
 #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[i+1]=DEFI_GROUP(reuse =msh[i+1],
 MAILLAGE=msh[i+1],
 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 = msh[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'
 
 
 
 
 
 
 
FIN(FORMAT_HDF='OUI',);

Last edited by Bracchesimo (2010-12-22 14:11:37)

---

Attachments:
OptiAster.zip, Size: 196.65 KiB, Downloads: 1,194

razvan.curtean

Member

Registered: 2008-11-20

Posts: 119

Re: Structure Optimization with Code_Aster

Hi

I try to run your example and I get this error:

   

 Version séquentielle de Code_aster 
==========================================
==========================================
#This script performs struture optimization by material removal
#
#Author: Luigi Giaccari
#Date: 22/12/2010

#PARAMETERS
max_it=1#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/SALOME-MECA-2010.1-x86_64/aster/STA10.1/bibpyt/Utilitai' )
from Utilitai.partition import *
from partition import *
from Table import*
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
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,);

#Since all other meshes are build on the first one we can pick connectivity data only on the first one
mesh1 = 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 
	mod[i]=AFFE_MODELE(MAILLAGE=msh[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_NO='TOUT',
		                  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_NO='TOUT',
                                            VALE=(stiff,stiff,stiff,),),);

 

	#Assign Material properties to Elements
	mat[i]=AFFE_MATERIAU(MAILLAGE=msh[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_ELEM(reuse =res[i],
		           RESULTAT=res[i],
		           OPTION=('SIGM_ELNO_DEPL','EQUI_ELNO_SIGM',),);

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

	IMPR_RESU(MODELE=mod[i],
		  FORMAT='MED',
		  UNITE=80,
		  RESU=(_F(MAILLAGE=msh[i],
		           RESULTAT=res[i],
		           NOM_CHAM='DEPL',
		           NOM_CMP=('DX','DY','DZ',),),
		        _F(RESULTAT=res[i],
		           NOM_CHAM=('EQUI_ELNO_SIGM',),
		           NOM_CMP='VMIS',),),);
	
        

       

           
		
         
	#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             
		msh[i+1]=CREA_MAILLAGE(MAILLAGE=msh[i],);
              
                #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[i+1]=DEFI_GROUP(reuse =msh[i+1],
                MAILLAGE=msh[i+1],
                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 = msh[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'







FIN(FORMAT_HDF='OUI',);

==========================================
==========================================
JDC.py : ERREUR A L'INTERPRETATION DANS ACCAS - INTERRUPTION
>> JDC.py : DEBUT RAPPORT
CR phase d'initialisation
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! erreur non prevue et non traitee prevenir la maintenance fort.1 !
   ! Traceback (most recent call last):                              !
   !    File "./Python/Noyau/N_JDC.py", line 205, in exec_compile    !
   !     exec self.proc_compile in self.g_context                    !
   !    File "fort.1", line 20, in <module>                          !
   !     from partition import *                                     !
   !  ImportError: No module named partition                         !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fin CR phase d'initialisation

>> JDC.py : FIN RAPPORT
JDC.py : ERREUR A LA VERIFICATION SYNTAXIQUE - INTERRUPTION
>> JDC.py : DEBUT RAPPORT
DEBUT CR validation : fort.1
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! - Il faut au moins un mot-clé parmi : ('DEBUT', 'POURSUITE')                   !
   ! - Il faut au moins un mot-clé parmi : ('FIN',)                                 !
   ! - Il faut qu'au moins un objet de la liste : ('DEBUT', 'POURSUITE') soit suivi !
   ! d'au moins un objet de la liste : ('FIN',)                                     !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! erreur non prevue et non traitee prevenir la maintenance fort.1 !
   ! Traceback (most recent call last):                              !
   !    File "./Python/Noyau/N_JDC.py", line 205, in exec_compile    !
   !     exec self.proc_compile in self.g_context                    !
   !    File "fort.1", line 20, in <module>                          !
   !     from partition import *                                     !
   !  ImportError: No module named partition                         !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
FIN CR validation :fort.1

>> JDC.py : FIN RAPPORT
EXECUTION_CODE_ASTER_EXIT_13373=1
.

Can you help me?

Thanks
Razvan

Last edited by razvan.curtean (2011-01-11 15:04:25)

RichardS

Member

From: Munich, Germany

Registered: 2010-09-28

Posts: 563

Website

Re: Structure Optimization with Code_Aster

Hello Luigi,
nice work!
I tested your script on SALOME-MECA 2010.2 and it runs without any problems.

Regards,
Richard

---

Richard Szoeke-Schuller
Product Management
www.simscale.com
We are hiring! https://www.simscale.com/jobs/

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

Excellent work, generous of you to share that.

delmas

Administrator

From: EDF R&D

Registered: 2007-12-12

Posts: 837

Re: Structure Optimization with Code_Aster

Waouw! Congratulations, Brachesimo! This is a very nice illustration on what can be done with Python and Code_Aster.

---

Code_Aster release : unstable on (Ubuntu Precise Pangolin 12.04 64 bits) - GNU + Intel

Code_Aster. What else ?

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

I try to run your example and I get this error:
Razvan

Hi Razvan,

You just have to adapt a path, and set the following to the equivalent on your system:

import sys
sys.path.append( '/opt/SALOME-MECA-2010.1-x86_64/aster/STA10.1/bibpyt/Utilitai' )

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

Come to think of it, it might be interesting to add mesh refinement to the equation. I hope the find some time the next few weeks to mess with that... I'm surprised by how quick different iterations are converging actually. I've explored some topology optimization software and those that I've messed with ( more academic code than production software admittedly ) was about an order of magnitude slower.

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

managed to get to iteration #34 before the "Matrice non factorisable" error.

---

Attachments:
iteration_34.png, Size: 45.34 KiB, Downloads: 859

Bracchesimo

Member

From: Italy

Registered: 2010-08-17

Posts: 100

Website

Re: Structure Optimization with Code_Aster

Hi Guys!

I'll be glad if you post comments on: http://www.advancedmcode.org/structure- … aster.html cause that's the place I check almost everyday.

Replying to your question:

managed to get to iteration #34 before the "Matrice non factorisable" error.

Increase the stifness of the soft spring!

Come to think of it, it might be interesting to add mesh refinement to the equation. I hope the find some time the next few weeks to mess with that... I'm surprised by how quick different iterations are converging actually. I've explored some topology optimization software and those that I've messed with ( more academic code than production software admittedly ) was about an order of magnitude slower.

Consider However this a very brute method to perform opitization and the code is way too simple. Mesh refinement will be good, but it wont be easy to integrate it into the code.
I think what this code is in urgent need is the element "revival": currently elements are only deleted but sometimes is good to pick them back.

I'll be busy for a while cause I am involved in several others projects with higher priority so at least for now the project is parked.

Thanks for your feedback I really enjoy to exchange opinions.

Luigi

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

Hi,

I've tried to run the script on a denser mesh. Its interesting to see how the resolution affects the results of the optimization program.
I recreated the element in Salome such to be able to recompute the mesh density while keeping all the groups (probably the scale is quite different).
I've used the following parameters:

max_it=50 #maximum number of iterations
DELETE_PERCENTAGE=0.015 #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

I got to 15 iterations ( the result is the image shown ) before getting the "matrice non factorisable" error.
Hopefully I'll manage to get a bit further with a stiffness of 10k.
Luigi, would you perhaps be able to explain the stiffness parameter in greater detail? For instance, is there a relation to the dimension of the model?

By the way, would there be a way to continue the simulation from the last iteration?
If I would chance the stiffness parameter at that point, that would probably discard the validity of the results, but would give some insight into what parameter I should be using.
Any ideas on that perhaps?

Cheers,

-jelle

Last edited by jelle (2011-01-16 11:59:46)

Bracchesimo

Member

From: Italy

Registered: 2010-08-17

Posts: 100

Website

Re: Structure Optimization with Code_Aster

Luigi, would you perhaps be able to explain the stiffness parameter in greater detail? For instance, is there a relation to the dimension of the model?

Soft springs prevent orphan elements (elments that remains alone) to cause singularity. Each nodes of the model is linked with a soft spring in order to keep the system statically solvable. This is a very brute way, but I think the simpler in code_aster.

By the way, would there be a way to continue the simulation from the last iteration?

Sure but is not implemented yet. It would be enough to take the last result in .med file and extract results.

If I would chance the stiffness parameter at that point, that would probably discard the validity of the results, but would give some insight into what parameter I should be using.
Any ideas on that perhaps?

You can put the soft spings stiffness into a list so that at iteration i stifness will be k[i].
But I think it would be better to increase them only after a matrix singularity show up. Unfortunatly I know nothing about error handling in code_aster. Something like:

Try calculus1
catch calculus2

Anyway consider the if your model is very stiff, deplacements are low and soft spring has a poor influence on the model.

Luigi

Frederic.renou

Member

From: paris

Registered: 2009-02-26

Posts: 154

Website

Re: Structure Optimization with Code_Aster

SIMP optimisation

I post this here as the same kind of algo is used. There is two main differences: I use a young modulus field in order to follow the evolution step. And I use the elastic energy to govern the material removal.

The main difficulties which come from CA limitation are the following:
Field manipulation under CA is painful
Need to use STAT_NON_LINE for ENEL computation
Can not define density or plate thickness as a field or function
Can not calculate an intergale of a NOEU_NEUT_R field ( should define a SIEF_elga field as result and use POST_ELEM)

More information here:
http://code.google.com/p/ogaca/wiki/SIMPoptimisation

I attach the file here for a 2d example

Last edited by Frederic.renou (2011-03-21 12:23:28)

---

Attachments:
astersimp.zip, Size: 426.57 KiB, Downloads: 901

---

fred

jelle

Member

Registered: 2010-04-20

Posts: 35

Re: Structure Optimization with Code_Aster

Cool! Thanks for the detailed info Frederic...

wspellf

Member

From: Belo Horizonte (Brazil)

Registered: 2013-06-14

Posts: 30

Re: Structure Optimization with Code_Aster

Thanks a lot for this work!

But in Salome-Meca 2013.2 didn't work because the new Code Aster version (1.13.1). Does anyone know how to proceed?

Pisolino

Member

From: Italy

Registered: 2009-04-06

Posts: 64

Re: Structure Optimization with Code_Aster

Hi,
starting from the commfile written by Bracchesimo i tried to update the algorithm to the salome-meca 2015.2 version. I already come through some errors but now i'm stuck with an error (that i believe it is quite stupid but i can't solve it). Until now i had to switch from Numeric to numarray, then updating CALC_ELEM in CREA_CHAMP. The actual error is driven by CREA_MAILLAGE, for what i understood something can't update during the for loop, indeed the error appears at the second loop.

CR phase d'initialisation
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   ! <S> Exception utilisateur levee mais pas interceptee. !
   ! Les bases sont fermees.                               !
   ! Type de l'exception : error                           !
   !                                                       !
   !  Erreur données : maille déjà existante :  A2         !
   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
fin CR phase d'initialisation

this is the code part:
#create the new mesh copying the current one             
        msh[i+1]=CREA_MAILLAGE(MAILLAGE=msh[ i ],
                    CREA_MAILLE=_F(NOM='msh',
                    TOUT='OUI',
                                     PREF_MAILLE='A',
                    PREF_NUME=i+1,),);

Any suggestion??
Thanks

Andrea

---

Attachments:
files.tar.gz, Size: 44.71 KiB, Downloads: 597

Pisolino

Member

From: Italy

Registered: 2009-04-06

Posts: 64

Re: Structure Optimization with Code_Aster

Solved ! it was necessary to use an incremental name (counter int converted to string) for NOM and PREFI_MAILLE :

#create the new mesh copying the current one             
        msh[i+1]=CREA_MAILLAGE(MAILLAGE=msh[ i ],
                    CREA_MAILLE=_F(NOM=str(i),
                    TOUT='OUI',
                                        PREF_MAILLE=str(i),
                    PREF_NUME=i+1),
                            );

bpaillard

Member

Registered: 2015-10-16

Posts: 28

Re: Structure Optimization with Code_Aster

Hi all,

I tried to run Luigi's script modified by Andrea, but there is an issue I reckon.

Incrementing NOM and PREF_MAILLE is indeed a way to solve the crash, since when using NOM='msh' and PREF_MAILLE = 'A', code_aster complains that the GROUP_MA 'msh' exists already, when going from step 1 to step 2.

However incrementing it this way adds cells to the mesh at every step, since a new group with an incremented name, and the whole mesh in it, is created at every step.

I've tried removing the newly created group by changing :
DETR_GROUP_MA=_F(NOM='model'),
to
DETR_GROUP_MA=_F(NOM=('model',str(i))),

In that case the group seems to be removed, but still the new group has an amount of cells multiplied by 2 at every step...

Is there a way to use CREA_MAILLAGE without any option, just to duplicate the mesh ? The documentation makes me believe that it isn't, but I'd be glad to have somebody else's opinion...

Thanks ! And thanks for that great script.

Last edited by bpaillard (2017-08-04 14:47:31)

z-1616

Member

Registered: 2020-07-20

Posts: 10

Re: Structure Optimization with Code_Aster

Hello Everyone!

I'm really interested in this code but with the new versions of code_aster, I'm facing a lot of problems running it.
Did anyone know how to make this code runs on the new version of code_aster.

Thank you in advance,
Best Regards,
Zaki

lihan

Member

Registered: 2023-03-17

Posts: 3

Re: Structure Optimization with Code_Aster

Hello all,

After changing applied nodes for springs and a few small modifications, this seems working. (updated for code aster version 14.8.0)

Yet the mesh is still duplicated at each step, waiting for expert's opinion.

Lihan

---

Attachments:
OptiAster_2022.7z, Size: 162.34 KiB, Downloads: 5

jacob

Member

From: Dolní Benešov

Registered: 2022-03-07

Posts: 147

Re: Structure Optimization with Code_Aster

Hello, changing the beginning to

#LOADING MODULES
#import sys
#sys.path.append(r'C:\Users\AppData\Local\code_aster\v2021\14.8\lib\aster\Utilitai' )

#from code_aster.MacroCommands.Utils import partition
from code_aster.MacroCommands.Utils.partition import *

#from Utilitai.partition import *
#from partition import *
#from Table import*
import numpy as N
# import Numeric as N
import math

worked for me with 15.4

however only to 5th iteration, then it takes really long to solve.

Last edited by jacob (2023-03-18 11:08:30)