Télécharger topoptim.procedur

Retour à la liste

Numérotation des lignes :

  1. * TOPOPTIM PROCEDUR KK2000 14/07/08 21:15:00 8088
  2. * TOPOPTIM PROCEDUR
  3.  
  4. ************************************************************************
  5. ** Topology optimization procedure
  6. ** Version 1.0
  7. ** Guenhael Le Quilliec (LMR - Polytech Tours)
  8. ** 2014/04/28
  9. ************************************************************************
  10.  
  11. DEBP TOPOPTIM tab0*'TABLE' ;
  12.  
  13. **********************************************************************
  14. * INPUT *
  15. **********************************************************************
  16.  
  17. * Mesh, model, boundary conditions and loading
  18. * ********************************************
  19.  
  20. * Mesh
  21. msh0 = tab0.'MAILLAGE' ;
  22. zact0 = EXIS tab0 'ZONE_ACTIVE' ;
  23. SI zact0 ;
  24. msh2 = tab0.'ZONE_ACTIVE' ;
  25. msh1 = DIFF msh0 msh2 ;
  26. SINO ;
  27. msh1 = msh0 ;
  28. FINS ;
  29.  
  30. * Thickness for shell and plate elements
  31. EL2D = EXIS tab0 'EPAISSEUR' ;
  32. SI EL2D ;
  33. ep0 = tab0.'EPAISSEUR' ;
  34. FINS ;
  35.  
  36. * Nu
  37. SI (NON (EXIS tab0 'NU')) ;
  38. tab0.'NU' = 0.3 ;
  39. FINS ;
  40. nu0 = tab0.'NU' ;
  41.  
  42. * Loading
  43. load0 = tab0.'CHARGEMENT' ;
  44. mltload0 = EGA (TYPE load0) 'TABLE' ;
  45.  
  46. * Boundary conditions
  47. ther0 = EXIS tab0 'BLOCAGES_THERMIQUES' ;
  48. SI ther0 ;
  49. bc0 = tab0.'BLOCAGES_THERMIQUES' ;
  50. SINO ;
  51. bc0 = tab0.'BLOCAGES_MECANIQUES' ;
  52. FINS ;
  53. mltbc0 = EGA (TYPE bc0) 'TABLE' ;
  54.  
  55. * Model
  56. SI (NON (EXIS tab0 'MODELE')) ;
  57. SI ther0 ;
  58. tab0.'MODELE' = MODE msh0 'THERMIQUE' 'CONDUCTION' ;
  59. SINO ;
  60. tab0.'MODELE' = MODE msh0 'MECANIQUE' 'ELASTIQUE' 'ISOTROPE' ;
  61. FINS ;
  62. FINS ;
  63. mod0 = tab0.'MODELE' ;
  64. SI zact0 ;
  65. mod1 = REDU mod0 msh1 ;
  66. mod2 = REDU mod0 msh2 ;
  67. SINO ;
  68. mod1 = mod0 ;
  69. FINS ;
  70.  
  71. * General optimization options
  72. * ****************************
  73.  
  74. * Volume fraction
  75. SI (NON (EXIS tab0 'FRACTION_VOLUME')) ;
  76. tab0.'FRACTION_VOLUME' = 0.4 ;
  77. FINS ;
  78. volfrac0 = FLOT tab0.'FRACTION_VOLUME' ;
  79.  
  80. * Maximum number of iterations
  81. SI (NON (EXIS tab0 'MAX_ITERATIONS')) ;
  82. tab0.'MAX_ITERATIONS' = 500 ;
  83. FINS ;
  84. maxit0 = tab0.'MAX_ITERATIONS' ;
  85.  
  86. * Filter rate
  87. SI (NON (EXIS tab0 'TAUX_FILTRAGE')) ;
  88. tab0.'TAUX_FILTRAGE' = 1 ;
  89. FINS ;
  90. frate0 = ENTI (ABS (tab0.'TAUX_FILTRAGE')) ;
  91.  
  92. * Convergence criterion
  93. SI (NON (EXIS tab0 'PRECISION')) ;
  94. tab0.'PRECISION' = 0.01 ;
  95. FINS ;
  96. crit0 = tab0.'PRECISION' ;
  97.  
  98. * Penalty factor (P)
  99. * ******************
  100.  
  101. * Start value of P
  102. SI (NON (EXIS tab0 'P_FACTEUR')) ;
  103. tab0.'P_FACTEUR' = 3.0 ;
  104. FINS ;
  105. p0 = FLOT tab0.'P_FACTEUR' ;
  106.  
  107. * Number of iterations to hold P constant from start
  108. SI (NON (EXIS tab0 'P_STABLE')) ;
  109. tab0.'P_STABLE' = 1e9 ;
  110. FINS ;
  111. pStab0 = ENTI tab0.'P_STABLE' ;
  112.  
  113. * Increment P by this amount
  114. SI (NON (EXIS tab0 'P_INCREMENT')) ;
  115. tab0.'P_INCREMENT' = 0.2 ;
  116. FINS ;
  117. pIncr0 = FLOT tab0.'P_INCREMENT' ;
  118.  
  119. * Increment P every 'P_NITER' iterations
  120. SI (NON (EXIS tab0 'P_NITER')) ;
  121. tab0.'P_NITER' = 1 ;
  122. FINS ;
  123. pNit0 = ENTI tab0.'P_NITER' ;
  124.  
  125. * Maximum value of P
  126. SI (NON (EXIS tab0 'P_MAX')) ;
  127. tab0.'P_MAX' = 3.0 ;
  128. FINS ;
  129. pMax0 = FLOT tab0.'P_MAX' ;
  130.  
  131. * Grey Scale Factor (Q)
  132. * *********************
  133.  
  134. * Start value of Q
  135. SI (NON (EXIS tab0 'GSF_FACTEUR')) ;
  136. tab0.'GSF_FACTEUR' = 1.0 ;
  137. FINS ;
  138. gsf0 = FLOT tab0.'GSF_FACTEUR' ;
  139.  
  140. * Number of iterations to hold Q constant from start
  141. SI (NON (EXIS tab0 'GSF_STABLE')) ;
  142. tab0.'GSF_STABLE' = 1e9 ;
  143. FINS ;
  144. gsfStab0 = ENTI tab0.'GSF_STABLE' ;
  145.  
  146. * Increment Q by this amount
  147. SI (NON (EXIS tab0 'GSF_INCREMENT')) ;
  148. tab0.'GSF_INCREMENT' = 0.05 ;
  149. FINS ;
  150. gsfIncr0 = FLOT tab0.'GSF_INCREMENT' ;
  151.  
  152. * Increment Q every 'GSF_NITER' iterations
  153. SI (NON (EXIS tab0 'GSF_NITER')) ;
  154. tab0.'GSF_NITER' = 1 ;
  155. FINS ;
  156. gsfNit0 = ENTI tab0.'GSF_NITER' ;
  157.  
  158. * Maximum value of Q
  159. SI (NON (EXIS tab0 'GSF_MAX')) ;
  160. tab0.'GSF_MAX' = 5.0 ;
  161. FINS ;
  162. gsfMax0 = FLOT tab0.'GSF_MAX' ;
  163.  
  164. * Other options
  165. * *************
  166.  
  167. * Draw criteria
  168. SI (NON (EXIS tab0 'TRAC')) ;
  169. tab0.'TRAC' = VRAI ;
  170. FINS ;
  171. trac0 = tab0.'TRAC' ;
  172.  
  173. **********************************************************************
  174. * PREPROCESSING *
  175. **********************************************************************
  176.  
  177. * Precision
  178. prec0 = 1e-9 ;
  179.  
  180. * Materials
  181. * Zone 1
  182. SI ther0 ;
  183. mat1 = MATE mod1 'K' 1.0 ;
  184. SINO ;
  185. SI El2D ;
  186. mat1 = MATE mod1 'YOUN' 1.0 'NU' nu0 'EPAI' ep0 ;
  187. SINO ;
  188. mat1 = MATE mod1 'YOUN' 1.0 'NU' nu0 ;
  189. FINS ;
  190. FINS ;
  191. * Zone 2
  192. SI zact0 ;
  193. SI ther0 ;
  194. mat2 = MATE mod2 'K' 1.0 ;
  195. SINO ;
  196. SI El2D ;
  197. mat2 = MATE mod2 'YOUN' 1.0 'NU' nu0 'EPAI' ep0 ;
  198. SINO ;
  199. mat2 = MATE mod2 'YOUN' 1.0 'NU' nu0 ;
  200. FINS ;
  201. FINS ;
  202. FINS ;
  203.  
  204. * Various fields
  205. un1 = MANU 'CHML' mod1 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  206. volEl1 = INTG 'ELEM' un1 mod1 mat1 ;
  207. SI zact0 ;
  208. un2 = MANU 'CHML' mod2 'SCAL' 1.0 'TYPE' 'SCALAIRE' 'GRAVITE' ;
  209. volEl2 = INTG 'ELEM' un2 mod2 mat2 ;
  210. FINS ;
  211.  
  212. * In thermal mode, 'GRAVITE' option is not respected
  213. * This problem is solved by INTG 'ELEM' divided by it self
  214. SI ther0 ;
  215. un1 = volEl1 / volEl1 ;
  216. SI zact0 ;
  217. un2 = volEl2 / volEl2 ;
  218. FINS ;
  219. FINS ;
  220. zero1 = un1 * 0.0 ;
  221. move0 = un1 * 0.2 ;
  222. void0 = un1 * 0.001 ;
  223. * Max and min stiffness
  224. Kmax0 = un1 * 1.0 ;
  225. Kmin0 = Kmax0 * prec0 ;
  226.  
  227. * Initial volumes
  228. SI zact0 ;
  229. vol1 = INTG un1 mod1 mat1 ;
  230. vol2 = INTG un2 mod2 mat2 ;
  231. vol0 = vol1 + vol2 ;
  232. SINO ;
  233. vol0 = INTG un1 mod1 mat1 ;
  234. vol1 = vol0 ;
  235. vol2 = 0.0 ;
  236. FINS ;
  237.  
  238. * Volume fraction for the optimized zone
  239. volfrac1 = ((volfrac0 * vol0) - vol2) / vol1 ;
  240.  
  241. * Initial topology
  242. SI (EXIS tab0 'TOPOLOGIE_CH') ;
  243. x0 = REDU tab0.'TOPOLOGIE_CH' mod1 ;
  244. * Add some rigidity everywhere
  245. x0 = BORN x0 'SCAL' 'MINIMUM' 0.01 ;
  246. * Filter
  247. REPE loop1 (MAXI (LECT frate0 2)) ;
  248. x0 = CHAN 'CHAM' (CHAN 'CHPO' mod1 x0 'MOYE')
  249. mod1 'GRAVITE' ;
  250. FIN loop1 ;
  251. * In thermal mode, 'GRAVITE' option is not respected
  252. * This problem is solved by INTG 'ELEM'
  253. SI ther0 ;
  254. x0 = (INTG 'ELEM' x0 mod1 mat1) / volEl1 ;
  255. FINS ;
  256. * Scale
  257. x0 = x0 * (vol1 * volfrac1 / (INTG x0 mod1 mat1)) ;
  258. SINO ;
  259. x0 = un1 * volfrac1 ;
  260. FINS ;
  261.  
  262. * Rigidity operator
  263. SI ther0 ;
  264. RIGI0 = MOT 'COND' ;
  265. SINO ;
  266. RIGI0 = MOT 'RIGI' ;
  267. FINS ;
  268.  
  269. * Multiple load condition cases with one boundary condition
  270. * Or multiple boundary condition cases with one loading
  271. SI (mltload0 ET (NON mltbc0)) ;
  272. tmp0 = TABL ;
  273. REPE loop1 (DIME load0) ;
  274. tmp0.&loop1 = bc0 ;
  275. FIN loop1 ;
  276. bc0 = tmp0 ;
  277. FINS ;
  278. SI (mltbc0 ET (NON mltload0)) ;
  279. tmp0 = TABL ;
  280. REPE loop1 (DIME bc0) ;
  281. tmp0.&loop1 = load0 ;
  282. FIN loop1 ;
  283. load0 = tmp0 ;
  284. FINS ;
  285. mltcase0 = mltload0 OU mltbc0 ;
  286.  
  287. * Number of cases
  288. SI mltcase0 ;
  289. nbcase0 = DIME load0 ;
  290. FINS ;
  291.  
  292. * Isovalues
  293. isoval0 = PROG 0.0 'PAS' (1.0 / 56.0) 1.0 ;
  294.  
  295. * Message
  296. MESS ' it. | obj. | max change | Penal. fac. |'
  297. ' GSF fac. | nb. el.' ;
  298.  
  299. **********************************************************************
  300. * OPTIMIZATION LOOP *
  301. **********************************************************************
  302.  
  303. * Optimization loop
  304. REPE loop0 maxit0 ;
  305.  
  306. * P coefficient
  307. SI (&loop0 > pStab0) ;
  308. SI (MULT (&loop0 - pStab0 - 1) pNit0) ;
  309. p0 = p0 + pIncr0 ;
  310. SI (p0 > pMax0) ;
  311. p0 = pMax0 ;
  312. FINS ;
  313. FINS ;
  314. FINS ;
  315.  
  316. * GSF coefficient
  317. SI (&loop0 > gsfStab0) ;
  318. SI (MULT (&loop0 - gsfStab0 - 1) gsfNit0) ;
  319. gsf0 = gsf0 + gsfIncr0 ;
  320. SI (gsf0 > gsfMax0) ;
  321. gsf0 = gsfMax0 ;
  322. FINS ;
  323. FINS ;
  324. FINS ;
  325.  
  326. * Stiffness
  327. SI ther0 ;
  328. K0 = Kmin0 +
  329. ((Kmax0 - Kmin0) * ((un1-void0) * (x0**p0) + void0)) ;
  330. SINO ;
  331. K0 = Kmin0 + ((Kmax0 - Kmin0) * (x0**p0)) ;
  332. FINS ;
  333.  
  334. * Material
  335. SI ther0 ;
  336. SI zact0 ;
  337. mat1 = MATE mod1 'K' K0 ;
  338. mat0 = mat1 ET mat2 ;
  339. SINO ;
  340. mat0 = MATE mod0 'K' K0 ;
  341. mat1 = mat0 ;
  342. FINS ;
  343. SINO ;
  344. SI zact0 ;
  345. SI El2D ;
  346. mat1 = MATE mod1 'YOUN' K0 'NU' nu0 'EPAI' ep0 ;
  347. SINO ;
  348. mat1 = MATE mod1 'YOUN' K0 'NU' nu0 ;
  349. FINS ;
  350. mat0 = mat1 ET mat2 ;
  351. SINO ;
  352. SI El2D ;
  353. mat0 = MATE mod0 'YOUN' K0 'NU' nu0 'EPAI' ep0 ;
  354. SINO ;
  355. mat0 = MATE mod0 'YOUN' K0 'NU' nu0 ;
  356. FINS ;
  357. mat1 = mat0 ;
  358. FINS ;
  359. FINS ;
  360.  
  361. * Active mesh, model and material
  362. SI zact0 ;
  363. mshA = (x0 + un2) ELEM 'SUPERIEUR' prec0 ;
  364. SINO ;
  365. mshA = x0 ELEM 'SUPERIEUR' prec0 ;
  366. FINS ;
  367. modA = REDU mod0 mshA ;
  368. matA = REDU mat0 mshA ;
  369.  
  370. * Number of active elements
  371. nbel0 = NBEL mshA ;
  372.  
  373. * Solve FE problem(s)
  374. SI mltcase0 ;
  375. resA = TABL ;
  376. REPE loop1 nbcase0 ;
  377. resA.&loop1 = RESO ((RIGI0 modA matA) ET bc0.&loop1)
  378. load0.&loop1 ;
  379. FIN loop1 ;
  380. SINO ;
  381. resA = RESO ((RIGI0 modA matA) ET bc0) load0 ;
  382. FINS ;
  383.  
  384. * Active mesh, model, material and resolution belonging to msh1
  385. mshA = mshA ELEM 'APPU' 'STRI' msh1 ;
  386. modA = REDU modA mshA ;
  387. matA = REDU matA mshA ;
  388. SI mltcase0 ;
  389. REPE loop1 nbcase0 ;
  390. resA.&loop1 = REDU resA.&loop1 mshA ;
  391. FIN loop1 ;
  392. SINO ;
  393. resA = REDU resA mshA ;
  394. FINS ;
  395.  
  396. * Objective and sensitivity
  397. SI mltcase0 ;
  398. obj0 = 0.0 ;
  399. dc0 = zero1 ;
  400. REPE loop1 nbcase0 ;
  401. SI ther0 ;
  402. gradA = GRAD modA resA.&loop1 ;
  403. cmp0 = EXTR gradA 'COMP' ;
  404. eneA = (EXCO (EXTR cmp0 (DIME cmp0)) gradA 'SCAL')**2 ;
  405. REPE loop2 ((DIME cmp0) - 1) ;
  406. eneA = eneA +
  407. ((EXCO (EXTR cmp0 &loop2) gradA 'SCAL')**2) ;
  408. FIN loop2 ;
  409. * In thermal mode, 'GRAVITE' option is not respected
  410. * This problem is solved by using INTG 'ELEM'
  411. eneA = (INTG 'ELEM' eneA modA matA) / (REDU volEl1 modA) ;
  412. obj0 = obj0 + (INTG (eneA * (REDU K0 modA)) modA matA) ;
  413. ene0 = zero1 + (REDU eneA mod1) ;
  414. dc0 = dc0 +
  415. (ene0 * ((void0 - un1) * p0 * (x0**(p0 - 1.0)))) ;
  416. SINO ;
  417. epsA = EPSI 'LINE' modA resA.&loop1 matA ;
  418. sigA = SIGM 'LINE' modA matA resA.&loop1 ;
  419. eneA = ENER modA epsA sigA ;
  420. * eneA = CHAN 'GRAVITE' eneA modA ;
  421. eneA = (INTG 'ELEM' eneA modA matA) / (REDU volEl1 modA) ;
  422. obj0 = obj0 + (INTG eneA modA matA) ;
  423. ene0 = zero1 + (REDU eneA mod1) ;
  424. dc0 = dc0 + (ene0 * ((-1.0 * p0 * (x0**(p0 - 1.0)))/K0)) ;
  425. FINS ;
  426. FIN loop1 ;
  427. SINO ;
  428. SI ther0 ;
  429. gradA = GRAD modA resA ;
  430. cmp0 = EXTR gradA 'COMP' ;
  431. eneA = (EXCO (EXTR cmp0 (DIME cmp0)) gradA 'SCAL')**2 ;
  432. REPE loop1 ((DIME cmp0) - 1) ;
  433. eneA = eneA +
  434. ((EXCO (EXTR cmp0 &loop1) gradA 'SCAL')**2) ;
  435. FIN loop1 ;
  436. * In thermal mode, 'GRAVITE' option is not respected
  437. * This problem is solved by using INTG 'ELEM'
  438. eneA = (INTG 'ELEM' eneA modA matA) / (REDU volEl1 modA) ;
  439. obj0 = INTG (eneA * (REDU K0 modA)) modA matA ;
  440. ene0 = zero1 + (REDU eneA mod1) ;
  441. dc0 = ene0 * ((void0 - un1) * p0 * (x0**(p0 - 1.0))) ;
  442. SINO ;
  443. epsA = EPSI 'LINE' modA resA matA ;
  444. sigA = SIGM 'LINE' modA matA resA ;
  445. eneA = ENER modA epsA sigA ;
  446. * eneA = CHAN 'GRAVITE' eneA modA ;
  447. eneA = (INTG 'ELEM' eneA modA matA) / (REDU volEl1 modA) ;
  448. obj0 = INTG eneA modA matA ;
  449. ene0 = zero1 + (REDU eneA mod1) ;
  450. dc0 = ene0 * ((-1.0 * p0 * (x0**(p0 - 1.0)))/K0) ;
  451. FINS ;
  452. FINS ;
  453.  
  454. * Sensitivity filtering
  455. SI (frate0 > 0) ;
  456. tmp0 = dc0 * x0 ;
  457. REPE loop1 frate0 ;
  458. tmp0 = CHAN 'CHAM' (CHAN 'CHPO' mod1 tmp0 'MOYE')
  459. mod1 'GRAVITE' ;
  460. FIN loop1 ;
  461. * In thermal mode, 'GRAVITE' option is not respected
  462. * This problem is solved by INTG 'ELEM'
  463. SI ther0 ;
  464. tmp0 = (INTG 'ELEM' tmp0 mod1 mat1) / volEl1 ;
  465. FINS ;
  466. dc0 = tmp0 / (BORN x0 'SCAL' 'MINIMUM' 0.001) ;
  467. FINS ;
  468.  
  469. * Optimality criteria update
  470. xOld0 = x0 ;
  471. l1 = 0.0 ;
  472. l2 = 1.0e5 ;
  473. x0 = zero1 ;
  474. REPE loop1 2000 ;
  475. SI ((((l2 - l1) / (l1 + l2)) > 1.0e-8) ET (l2 > 1.0e-40)) ;
  476. lmid0 = 0.5 * (l2 + l1) ;
  477. tmp1 = xOld0 * ((dc0 * (-1.0 / lmid0))**0.5) ;
  478. SI (gsf0 > 1) ;
  479. tmp1 = tmp1**gsf0 ;
  480. FINS ;
  481. tmp2 = xOld0 + move0 ;
  482. ma1 = tmp1 MASQ 'INFERIEUR' tmp2 ;
  483. tmp1 = (tmp1 * ma1) - (tmp2 * (ma1 - un1)) ;
  484. tmp1 = BORN tmp1 'SCAL' 'MAXIMUM' 1.0 ;
  485. tmp2 = xOld0 - move0 ;
  486. ma1 = tmp1 MASQ 'SUPERIEUR' tmp2 ;
  487. tmp1 = (tmp1 * ma1) - (tmp2 * (ma1 - un1)) ;
  488. x0 = BORN tmp1 'SCAL' 'MINIMUM' 0.0 ;
  489. SI ((INTG x0 mod1 mat1) > (vol1 * volfrac1)) ;
  490. l1 = lmid0 ;
  491. SINO ;
  492. l2 = lmid0 ;
  493. FINS ;
  494. SINO ;
  495. QUIT loop1 ;
  496. FINS ;
  497. FIN loop1 ;
  498.  
  499. * Compute the change
  500. change0 = MAXI (ABS (x0 - xOld0)) ;
  501.  
  502. * Write iteration history to screen
  503. MESS &loop0 obj0 change0 p0 gsf0 nbel0 ;
  504.  
  505. * Plot to screen
  506. SI trac0 ;
  507. SI zact0 ;
  508. tmp1 = x0 + un2 ;
  509. SINO ;
  510. tmp1 = x0 ;
  511. FINS ;
  512. tmp2 = REDU mod0 (tmp1 ELEM 'SUPERIEUR' prec0) ;
  513. tmp1 = REDU tmp1 tmp2 ;
  514. TRAC tmp1 tmp2 isoval0 'NCLK' ;
  515. FINS ;
  516.  
  517. * Convergence test
  518. SI (change0 < crit0) ;
  519. QUIT loop0 ;
  520. FINS ;
  521.  
  522. FIN loop0 ;
  523.  
  524. **********************************************************************
  525. * OUTPUT *
  526. **********************************************************************
  527.  
  528. * Topology field and final mesh
  529. SI zact0 ;
  530. tmp0 = x0 + un2 ;
  531. SINO ;
  532. tmp0 = x0 ;
  533. FINS ;
  534. tab0.'TOPOLOGIE_CH' = tmp0 ;
  535. tab0.'TOPOLOGIE_MAIL' = tmp0 ELEM 'SUPERIEUR' prec0 ;
  536.  
  537. FINP ;
  538.  
  539.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales