Atom topic feed | site map | contact | login | Protection des données personnelles | Powered by FluxBB | réalisation artaban
You are not logged in.
Pages: 1
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)
Offline
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)
Offline
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/
Offline
Excellent work, generous of you to share that.
Offline
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 ?
Offline
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' )
Offline
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.
Offline
managed to get to iteration #34 before the "Matrice non factorisable" error.
Offline
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
Offline
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)
Offline
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
Offline
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)
fred
Offline
Cool! Thanks for the detailed info Frederic...
Offline
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?
Offline
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
Offline
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),
);
Offline
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)
Offline
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
Offline
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
Offline
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)
Offline
Hello,
I am glad this works, I fiddled around all day. The ever increasing mesh is a big problem, though. Also POST_RELEVE_T is VERY slow, the larger the mesh gets. The only way around is keep the number of iterations low and use a coarse mesh.
However, instead of doing the CREA_MAILLAGE
msh(i+1)=CREA_MAILLAGE(MAILLAGE=msh(i),
CREA_MAILLE=_F(NOM=str(i),
TOUT='OUI',
PREF_MAILLE=str(i),
PREF_NUME=i+1,),);
we could copy the mesh with
msh(i+1)=COPIER ( CONCEPT=msh(i) );
This should keep the mesh size within reasonable limits. I had to use () instead of [], please do not copy the above text.
Mario.
Offline
however if I understand it correctly this is not SIMP, or some algorithm where the stiffness of the elements just changes and there is no actual removal of the elements. I think that would be better.
Offline
Hello,
for sure, there are better methods out there, than just deleting elements with a certain VMIS. While I must say, I do not fully understand what is going on in this script, the number of elements is reduced.
Otherwise, there would not be a (local) maximum of VMIS in the attached image (arc in center of image, where elements were removed above it).
Correct me if I am wrong,
Mario.
Offline
Pages: 1