From 8bd7efb48198f4af9dc51eeb9d7b78eb49404376 Mon Sep 17 00:00:00 2001 From: Su Tian Date: Wed, 12 Aug 2020 13:48:36 -0400 Subject: [PATCH] 200812 init --- __init__.py | 0 __init__.pyc | Bin 0 -> 157 bytes __pycache__/__init__.cpython-37.pyc | Bin 0 -> 157 bytes __pycache__/analysis.cpython-37.pyc | Bin 0 -> 3034 bytes __pycache__/cross_section.cpython-37.pyc | Bin 0 -> 3081 bytes __pycache__/presg.cpython-37.pyc | Bin 0 -> 6285 bytes __pycache__/sg.cpython-37.pyc | Bin 0 -> 3888 bytes __pycache__/utils.cpython-37.pyc | Bin 0 -> 6451 bytes analysis.py | 124 +++++ analysis.pyc | Bin 0 -> 3625 bytes cross_section.py | 131 +++++ io/__init__.py | 0 io/__init__.pyc | Bin 0 -> 160 bytes io/__pycache__/__init__.cpython-37.pyc | Bin 0 -> 160 bytes io/__pycache__/__init__.cpython-38.pyc | Bin 0 -> 164 bytes io/__pycache__/iosc.cpython-37.pyc | Bin 0 -> 11875 bytes io/__pycache__/iovabs.cpython-37.pyc | Bin 0 -> 9501 bytes io/__pycache__/utils.cpython-37.pyc | Bin 0 -> 2418 bytes io/iosc.py | 639 +++++++++++++++++++++++ io/iosc.pyc | Bin 0 -> 14796 bytes io/iovabs.py | 482 +++++++++++++++++ io/iovabs.pyc | Bin 0 -> 11394 bytes io/utils.py | 88 ++++ io/utils.pyc | Bin 0 -> 2946 bytes ms/__init__.py | 0 ms/__pycache__/__init__.cpython-37.pyc | Bin 0 -> 160 bytes ms/__pycache__/analysis.cpython-37.pyc | Bin 0 -> 2710 bytes ms/__pycache__/beam.cpython-37.pyc | Bin 0 -> 2603 bytes ms/__pycache__/iogebt.cpython-37.pyc | Bin 0 -> 5907 bytes ms/__pycache__/prebeam.cpython-37.pyc | Bin 0 -> 7564 bytes ms/analysis.py | 105 ++++ ms/beam.py | 98 ++++ ms/iogebt.py | 419 +++++++++++++++ ms/prebeam.py | 475 +++++++++++++++++ presg.py | 329 ++++++++++++ presg.pyc | Bin 0 -> 8074 bytes process.py | 81 +++ sg.py | 241 +++++++++ sg.pyc | Bin 0 -> 4843 bytes utils.py | 281 ++++++++++ utils.pyc | Bin 0 -> 8526 bytes 41 files changed, 3493 insertions(+) create mode 100644 __init__.py create mode 100644 __init__.pyc create mode 100644 __pycache__/__init__.cpython-37.pyc create mode 100644 __pycache__/analysis.cpython-37.pyc create mode 100644 __pycache__/cross_section.cpython-37.pyc create mode 100644 __pycache__/presg.cpython-37.pyc create mode 100644 __pycache__/sg.cpython-37.pyc create mode 100644 __pycache__/utils.cpython-37.pyc create mode 100644 analysis.py create mode 100644 analysis.pyc create mode 100644 cross_section.py create mode 100644 io/__init__.py create mode 100644 io/__init__.pyc create mode 100644 io/__pycache__/__init__.cpython-37.pyc create mode 100644 io/__pycache__/__init__.cpython-38.pyc create mode 100644 io/__pycache__/iosc.cpython-37.pyc create mode 100644 io/__pycache__/iovabs.cpython-37.pyc create mode 100644 io/__pycache__/utils.cpython-37.pyc create mode 100644 io/iosc.py create mode 100644 io/iosc.pyc create mode 100644 io/iovabs.py create mode 100644 io/iovabs.pyc create mode 100644 io/utils.py create mode 100644 io/utils.pyc create mode 100644 ms/__init__.py create mode 100644 ms/__pycache__/__init__.cpython-37.pyc create mode 100644 ms/__pycache__/analysis.cpython-37.pyc create mode 100644 ms/__pycache__/beam.cpython-37.pyc create mode 100644 ms/__pycache__/iogebt.cpython-37.pyc create mode 100644 ms/__pycache__/prebeam.cpython-37.pyc create mode 100644 ms/analysis.py create mode 100644 ms/beam.py create mode 100644 ms/iogebt.py create mode 100644 ms/prebeam.py create mode 100644 presg.py create mode 100644 presg.pyc create mode 100644 process.py create mode 100644 sg.py create mode 100644 sg.pyc create mode 100644 utils.py create mode 100644 utils.pyc diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/__init__.pyc b/__init__.pyc new file mode 100644 index 0000000000000000000000000000000000000000..02f92be9ffec9fe71871da698646012ca1bc72da GIT binary patch literal 157 zcmZSn%*&OhAs3g-00oRd+5w1*S%5?e14FO|NW@PANHCxg#cn_`XRDad;?$zzn9{t= zlEj$e#GIU%Vgn=H;*!#o%9vs!-J;au)Wo9XjF_V0_{FKzg`kf;VCv^Y- literal 0 HcmV?d00001 diff --git a/__pycache__/analysis.cpython-37.pyc b/__pycache__/analysis.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..bcf268d782099c3388217b56b44818885c71ce92 GIT binary patch literal 3034 zcmdT`%WqsY7q@*MGjo&4lcoWc53dNbNI*qMC_)I@k|IfFS+w z{pR7XU4(uym05D&<2}go2?&N5PEpvDIK;M2LSjoNbZzN{zAal>JM7@5Z;&m7i_`CJ zwnPv{*>bot)vQi6YvDR*I?Q426CA$40_HIvzOJwqYd=BZ2Ezxax9}UR1nu9AT~k78 zL-rxxge<=Wq0lL=QH^VI=Jd(Lop_4WuBq2v?L77{QVzq9F@`>j+$UmC6Xrnqk#bMH z+GXxlbcQE>z6p8{^he4+Z8iGZJ0lan_79A$4{c?>$yznO?~8AsSKG`3Jz90^0I@bK zw-fyrBNPM_(x;J#GOoCgwz@q}wkn5*kFu2RQmKSh7?ZsuHtP;rr!Ta9bSo90Zmso)Vy@*$M6N|oH%+0l5UU1K;VJNI>GZm2QS zc(C8=Xs47z?UjNb+}4h$a`T#ors+?;7w z=GDE%tG%}No4cdiv$F?{-M((gSmbcV+Ji0$Oe!;$WSYoKx2CKA35iWW%^!O=cl-Aw zT+V(4#}@Tvl%{?8#`SHfDt6MB*SEodMAsz4~V8zr5-e5r>v}E;7>ol#C33&S7J*I=cmz zK=)6X2auQo=>BE&V{itJ0Eo|d>{gQTTmr4mP9ck80T3_BBsLsY#i|lOml+&R+MN5Y zWqLYA@h<(;a0Z=nTgIh~hP>gsU3F3d!!~f=rAI}P+KI!Zi>y7~yjZ7E$wg8CnHom0 z93HM*_VxIs%bIhu@wI2#hM62qW77E5XIkeI(x>}v{#mdXKhrV`5aX@OHnX5G3pCS! z&!W#vl;qzluv-_ z17RLTp_3$uFmDizZ`nb`%TW0P(;m?z&NYd0jUPhAmx>t%HYFEgzqhXas4RKTv;*rj zUN6%q=Gv=3E3|uD08u;gM2amly!TKfifg}#@!A`tl^j9)QK?;szo~X1R>FLR)Gfvz zDlkdqEhSgIuhFN~{@i#fWtfVRr`T!?ei3 zQ^57wi8H3Xz5RQ4?r6+yf;#;gP4O~~?H5NLxC0!)AL*CsBp(E#L{jRH#~a9kYV z6%yb@vWPd@KIrXOu^DET;2qV(2IF;xfGhF-QwuOb?0hM@4d-_wQt@aEdTEXy5n-0q z!*|X9s(1@3V#kQLLFhJvr+@;8YQmHbWOY=E0$vzWTsO@NHqz}t0H4NjITS5x*EE2Z z82`v&rn`+QDRyiGCrP1Q^CjSUjr8(Utv2%e%?m)7;Acak-Md}v<^35i>-Iaq5LL?G XH6D>0Ae!)M$f)DN>2oV$7Mq0x%T+;s2WoOfoH z?peca&QCZIkNFeqyyd^iL!Rap31#1SMk3CyYWDW_EJx{4|GKLBW2&mFtNT{3*RfEF z>#uxiowuy-i8(1Qh?}Tp2#}T(V=HP}5eeq)$mZ6GTyDMGw=H4ittb%I9ZTBMd1y&T zwZCtO9l@`?hs=`0!7?dl*)?-k0J+8-x*RwT8T6r-`Ky?r&xirP_y?2FnpYCQahrPr{) zPru*!c<@J~v>8l`tV#wZ8IK3%z3m&On#!Yr*}kEbQHf3u2HKAC?;9b%5V# zQdd)KX^(((=?hI^$WCg2%&A*(Y$_ec=fD`6VF_2vdMAc)AMcE&5J_n{D%4VC3;%?j zJBEditq0=Rk>Z}MzmxVOaV#iu{$gLT9#bIL;N2?}22SOI_O4jplH5DCrH?uLEcF*s zpSjk8YZX^9s)bR1WqELBRC_tWL{(cPyV`#`fO=Xb@^@d{xxM=thF-)f z<|BY50(|F1N35e=7aI%gu#EahvVr~%BT+PZ_CgbSXI9|THb5M6IvmjsBOPL%fMK;u3R8kZDm5^aPbsD0r`08k zK15v{Q(tiYnNF7ZXBTIj)a|?ab(;SM!Bqlw#jcO;T#J$)i)Uj0&~nwH;LRHm5RdP;++)zQJ%2(xF7c#$9%awvUqc3GGSr z#C}Hd7GuEJmKYUx7}prPj6LZ=zRvV>Ja>cf9N&=7Grquhp78?XMaD~vml>}xzR36z zd@%XDN^qVq(hSylq=H9p11J@b*vV-2YIe^!8 zIpAi_rWx@G9l?oQ7Q~OAYC1Ov>q}^uj(l7y*oirvx;aIA|CgHLJ?in%tyDf4;|Z`} z_Y>Ui$ShHZdvQ#rlc+N+;`_$8{{3-c%KNpQr?ovwBfCO1 zY3=5S(}*yQTqRWjOifOh5>`xCX_(s7fFUbs||PrZ!!bT8)7^1hj?IQ_y1p}xrH4u ze*oacV>fS}uIM_8*Uv)29MR8`L)x`5M*#8pHf9577gkfnO7ila%3Aa`k!2dk|2uur zBlFj!C|!Di=h~Vi{a)mx276P&mFCFQ9s*->q&Lx5U-%{;Zz@Ef+PvL_#9h#}$D(>| uRc5QI_M4o<_X1}}{=-UY=B&z3so<05D*6pw$Yz~1!WBo~^UlVN4gbGopUkHK literal 0 HcmV?d00001 diff --git a/__pycache__/presg.cpython-37.pyc b/__pycache__/presg.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..2bacc7d2ee362e241ebf2888eafae5aa5d1f5738 GIT binary patch literal 6285 zcmcIoO>7%UcJ4ot%^y(|^=sLd+a6mIXD!*5J;A_n29xpFv%B$n79$V5>3GAmSWT+s zW|OWi*_2R?1nW&s308;29)dYgL9o~aL6BpPLC!hk5bU9mLx2Fga1J@-=u^H|-Q?KL z0LdjJx_(~Od-dw|tM|UHPiAH+8h-OXS!w_Js;2!rbtXR(jgRrfe~N->Ob<0Tt94gb zeZw_XZMg+i7yBhc)3v_smi26ArC)Vxd3&Z`SEI9TBd5&el%_kMQx^J*?h@JtGk>hP z%gka0)GMsQY*zlU?w(=#uGXr)f=1fIj!qSRGk8jP?&69615~Dy1ZI6@CORvyA}f7m zB!~C3Qt9XOuIb9TZsQjK<1$G*QG@iB(Va1Doc*)Tg7> zJDM7+o$84-rIx{`KQ-c?y+tve8a$Q8o3t}pQW*V0qr_$nQIqydeOXJ&Vn)<+dDv5?_jne%8!z>w zA}g{aZD}*@kBzPYx&QfC8+`;>2BaB~c0I92o9OG1NBz~L3XZC@y9P8N-YW5(=+7QJ zrt;2XqE2i!ms_QI3Lkk|KGFHFl$Fp=YO)4f&QC}FQ8Xr$H>^9y=3)124>W(~-92Ef zIW~iJO!&gWLr!k+N)-lLEz>Or@j&pwFzC3|0T=tzju%ZwT1Kk-sSzEv%BkVsz_a5T z{!Q2H-PoamT01vgD>_7^P|iTP0+h$FY;hh8j# zPFn5^qF8v5h+D?1BIa=NR#l`XAG}D*nAYt^ejKNEEVv)_#J*eZ@Id%H7(}gRYR0|3 zJNrG42wu4BcLeCE+4DtO07o$7HmzOuyWRE>+m6#hHyn5(wYovXpq@Vx*y%7W9Q%BL zQSL=OA7TYhtuTm1YW6*h#=|f`JBS2tKwxVoHAjBmEz|be-O%f~Gx@ICp||IUsosal zhHlv(W`sC3h1Vm&-Bb_UBDRMOn*5AKD2-;YEMi^O+OXojI1x@{VUUAOgyaZX<|D*9!;saHT>H!FHwujy<0g5JR2 z(d$M7BMp63zijYL@OS8C)9(X3WNQ^@{u*x+u8rdk*EWPcHl)D}sh^sODJaNT%!KO} zC}gyYImMb%N;#!4rPw)zfl(WIYAWipR^wT=h!vM%!)3N~Mk7wxZ$&K2<(Ei}Pc5wTPw=f}hLE0*E{!gW6}hAm z%1UxZ&W+EIz2yp9{tAAqYbVCDUn1Ze&;9_3=otF0Kwpy|KwC=%>CE@)xmN!}LpATr z46mYM$4g=B>7aYp~f45ytNba*`GYtpr?aVyQVbxh}p?S%rV;+ zqk49hq2B{`0aWs2C--EF5ApVL%Jc*{Iql-aig>G+)T~QWzF*`q1C%M>MeJh^bKjF| z?7g#Eav9!o1vXrlm#dm|pzDU*=n=HJB(L1jM*pC!`uF5}%DXm&Dc+Nt^3qxO5#+SQ z9=W-!$rf1;QSbfHF+8ZHJSa-OC9lYf%Ybq#T5?ssFTX{zV7qLUtFY9%vidsa zT~!ucm+PRK2OIFbb@Di@eqOzc&3wekDDc_IjH}6&lfttV!~vyCvT>eMB*U`S7w{OVnQ~C~cyP)R(R9;gOZ_1mhz9nzu*0}~dT+iD# z(ca11x0HR9_E`OvT0P$l`ZuW`)ux<~0^9r@tU*BO%0mScGT>r+kry7t0l+{XX*EO+ zN7Brad_@2s#96y6{QgiiIpG^1E2Qdn=a4^7OF`5L51F6Z{T!JbTic_4=oUk-A4FbS z7W+ZxDFGW`CVX6Id<#Wdo=nHup?7pROpP#1i+vwBjVY4PU+cFO%4(T!FcyGXLO?bE zC5x~?Y7B=xNK`+53zb{y4?2VX-hfxodsRRIEoe6(RCF7?aIoiv?Yh_jwZ)pF`o!>r(%(cH5Bd~*&2XPnOn^TNCo5@dOY^i65t;n z3`ANYkO%0bkV=LG&CrkdDrBTZZ#eXU0nI}Zx^>0*fNXvj`*qid^Df!zlcNcQ`kV)V zP97mlw@KrZflrPe0_*Z^(#+1y9J^)0n`*mmVbJ$`URqG51jZa3MvSkJ5_T9w@zCqQ zRdP$jZY^&GEWW{V?;}~M_0&I#TOXu_5pTymfYLq-`q{P&^o^jC>%f3H$S<)7 z`a(WN^eL?X9?Z2sspu8V0H(HTs-MlNqkYbxk(yb{W&vSWio^#_11p(z)OF}W53Sp@ zGMP3>2Iqo7YcCkKVdv}B)VzkV1k_#wf~OtS@Q~hhAbXOP{YhKkeaxbN1$(tss|JZp zqiOKpB@>n77eQAbhzuWfD1Ok7K8B{K|BqlcY4Mr>?th7REA#_n477i$KQWR5&~-t@ z@{+K|1)?VS)E*aslYy5HO3Vay|0T1M3SCP`t!E}xR+7fkBHxqMu}#+!h~Z6ERy3h# z!0O0_G+LGTz^5Lojmr{~B z*>U4kW3$M0GuS~xf(z2-cD3l|M@=yYOy6YqR7Z}l19h!inwWoD!AE>O=LdS7&dE7{ zpiMb}S;1L(i}Tx4&PByZbEcfK&Mf3Yo&UaAkY!oH9%iM9{zxp2mw-P@r869Mp0e1#n$8+#`nT2TwVc4%!9Q?lTaYQbBJjbfLHT|h>+-oCg zIk$;b5s2k}5c&>%+mqDY&m1I&pcgrapSdJw6W^H*#d1%y6x$@Bz9Z?;zVA$j#J(q- z7eN?0d%n~0IASex28Y7&x=_y%m`iFSF(CY6>8U<65=2AH>_W5cDShABa`p%P0l0!= z#Wfvb7-Ih0gV(}A#|!hB*SOyqJomrebQc_X%5l105F$B|QPhRRRBYG(hFN~}+z$ss z|La-7H_SpZ3Gl994mOqRsJ&({CLV-P=N6@wXdz=Y;HXF_Tmi(KCm%8 zy&KtEIle-#bIU>C%RMECdhdA2_V)I%xgFznql4R#xE*H~ED9*O6f#Ntv9T4TmTDEY zsNyc7!m+s(cV6pTuXQKYU*Jlq7z>6Ms%oEe1!9iv-@Et47oUIeL+2sVSr6&#_`?b1 z%gEbI%UKk6VO<8&yTL?|jPKf)!_J?|rcg^qFux2k499K1g*#T_2OzHSM5e$R`2)qv)aYl8y*AmYD6T7O0b84G#1kJNhq z;g=5{aKd%mqvAVM&^6Q>=@EuOdbAPX<^DJ@U)2N zv<;}|``h<2u=)t?w4g|7^$nVOd%$LNK?ZluwwR&Uj1~;G$>_jP1J|M!AJA&EZ;X;> zr^R5P>Jr`&HQZ`8<4^=4K)b5=Zcm@_kHE<5M?4eu;3`vyw! z@#gRU5*WrmXmD5@bZ(*QOCZ7!W@d!u14CHCerX6>S}*L-#@G=pj9WZ*gt=w(+_%_a z^b?|o#YX2As{RRx#7V*kEn$YX@P#d$mu84_g!|G6TcRyIw65sjtfl&$hjArU8fRNF zsnW8TJPvLLd7P+HC*>$j0$r(bQjL`i@>0mm_XCvc_kyG>bd^@)Dt#uIckNy<%5c~x zdRJx{vp&5S?8rFh{#&|ZVu4Vh)(#~eJ@UXk{%Ep;EkyZ+yN479X_b8)DPOJUesrE;* z?m)D`ppaVEep2S6ERBlp7~0kEk{gl38`=n^x_XDjlXp=;t^mkp=9~Yw_C9A=c53a0xn%)ohKeNR@9%g>U)# zHtwQQGGkeUon$vwag@oT`HoY!bS~1oc8WNM2y!rpMyed?+Jk9jF|2lR7e~u%pjL6@ zpAB^$&3|s)y|Mk3mP&7ri?oWjb)03}`s%gIx*Ci9ZGG*slA12HyRCGx8zt}#y`Agf zDBae>E2I6|i=q@;qv$=Hr>}t+rfa&6Wi~(8wk*p8{e_Ht5chSE&r5J1X zYtx^6aDOnMsCXs=(wJ0LDs@l}f@)U=TMw{c>v(Vk&@=f_S=_t^FDSw`aX!K6(h1J!^yXhL-j8S)sU3db(E(xKDce1{!W-rG?Wx zI_>#|sU({}f(NpY2<3DlMer!i zON3+bqzvZ!kU?4_l@^uSpiLAZ?Q{% z2x7qCFuBz|GCFfs%q|?kGgr*hE-*Wz-!ob2UYlZ&!aj;>KZ=mR#~JmzQS@{iXR{f2 zn<$egQtSAtqF`0uCqijpfskh^3N{i(QzG+R$Mmmz%gjx{m={`h3sqkQIfh4!Z$Kl) zHK36z+Kfh?=%8&2pHXRfZd20y zD3B^h$cd;>OZVbD@i%lD$eCIMIa3A+csF8X9r+@%1@QqdpnxDymuZ`+Cun^K*lG`H z&WAirZz-kSU8c*tc7^FX%voi6W=Lz!GJTinIi}}%c8%#e(+i@7l^ZvE_JJEyxBE$=+Ano z&~^8y|EUNUM%Rvwt_xrw#9$DtJCxY@d*KpufyofCZk4$l#^JePR_;I?k!qlN!ourt zb;ea?71NtZhi7J-tc-<@=m`;a7g;A3;J9j*71$!`$1qEiHI9Ogb70qZMzo06xBC<$8EWeYr8J<2xI}NTrF-*_pED8){X_$X&|LkM=uo)kUydGEr4+d zAjg0M>l_9Wy8)7geXZ9f(!a6SdI12QfsQ-lxeL6*x%OV4?H0h1%&xJbNDh)`Np^v^+E20^-~ zoP07@Qxg9>s=ffyzzMw}2t3Z=*U-uNxsoW=UYSA2{Gm=+(?(^dpnn(7KcX!8+I&N} zy|!%~SO@mh-m|7AUIpE$wPif*0>SAOW%Zlqwar(pTX>@8njL=8`D-(+&HTjW>$7-J zi^C_*o&7#Q?4S5L^VcT5vfJDD^TJ-EEXn$?6K9h+uIO5dEdy>UIIe5iNzpoj5j%0f$YGY85w#Mz zq-IwPf9>ZI(M2m zZZM;nFORX5v(OwbPoT_&`L$7=6t1v^^EE3M#p@y`@+h4WW1@hzDe;Dw5R;&s7UQBQ zPSwpVCZ@z`P^QJSm;uj>m=$M`&x$!QkNk{yQ!I!@aLkFb;v6XR;uY~K@`d*eabCRk zwIN=Ui%-q+S-~C{r3=ro5@VyrD8rBAGx4qC3pbHO#xo|^A?vXib6zyYjhG4ZhJl*3 zNM&1>(Sow2%TC6X>vQQoyx#`CJjyis*ptnIi{DqB%9s95&wbZjmQ>p9xtH8dWw)SP z*V>-)nr>&!z1F?zc4RGTsk9+FXiK6u{4jD`%Jr1;4wf7KuGH-%{+!mLXi+Frrf8uh zuCwP^nND?>uBcPS)2C~mjIL!8p~Mp&v=jS@RILyd)+sxk$ZrVwt)cQ@I9gv{+x#e$ zD%{);{K(r3y+&g*yt29+M*HGmGhAI(GL)XGZEvcuwq2>UgbX*EVZH5d?n6?z(mqIR zEpiO#eMvaz7zZ^aX-#WhV???tm$g%Z&m#ZC6cx>#Lv2h4B zTr_&7uwql#=Zpp7wuE!YW~l@&;saOioH1`ic5FOl3r5fC^UoVSJGNrYV`J7_`Q7h}d8YBb@U#!OS<7%3XHvV%BQ;!fJQZ0YrH)OEe z>iD&xS@p*bBjpEm-2fHVDo;t7OEs+0iuAXR55K8fnlGyvx61sBrnZ5M;@*0$wcSxR zxE1MEVv1Iju!o5$52RMWQm)KZJ8{x^By2rlQP~Wd3V;|i6JEQbEh0@wmlL;!^XhN- zLw!*1;hgT1vFp=02d)qhY#8GOmgm##3Uod#O|W-R&T}=3)*2Zk{Ydj^idysUp&(ey z=K)*L9#fb-t9w0yK0jr5jK`NYdC!I}+QQnk)OF}E1E)=ThO+sL#jt@Nptg7Rkr7+a zGH2t7`Q??hBi2*nkgbv7h`F%)F&XSPrGn-q&Q?n`y$HKuZ;_HDV?s9krVnjP%&lfr zwnC%{3lruiEI7U=6V&Hn>w!MRcQ8s=L}FmA7Mo^yyXuXodenljM+xXN-Ru%JIN_e_+AhuX03%J^aHHY_JVptCaiXR z(x-I8!`6BWjB@qPZOkGh_%m36&$2nL-a=Vx+yJCVp-2n#QT-VOnzjHu30^k^O%u>F zGWC?I7t(L7n{^WlC0H{7J6~xToTQKp{2B60rsg!kfy3RgbL z)aV8(dMC*5^X}z{+yHyZVK(PsHZY{fjLq133KP;hX~%jW_vTj!S8m_xwV??RR{BCDl@R3u6wG(gd9ZMn3w zl+%0$PNLN&+l9LrlA3iM5~^ghuYy$GcP)j0eN7S8uepGP_Q8riY%Moc0z1XfQ`q6e zeB!s$^B@%{If;|5C9(DBS_8C3sPCb@%xnH}WjI61@P4|leI&hhrv~m&2eA_cHVed@ zVrm(*!C{l}4GudE@^HLlNKo}fpeBMjozCJoHHE@E9sZFg1`YLNTtpe+04Uu2iKS{Ys% zEsNZM#+51xq{jm|1kh0=3mTTCb^B%nq|gpwwP!-p`!yFjPFR63g(BM+vF;(ZN7oSL z?p=Rm>@8*5ZhDIxBMzbuQZC^u%$q_GASh7LO_klP%IFIz14JZna<1okfpDvt+7k4x zeZJlgg}sWu>2B8GXgA$zW#{t|<@WkPFX=Mbisx1V3Zc8?Rjyuky~-*j*M{Sktjma) z>TxmEzL+YBsgjuLsF9*7hPz#_*>1?SRPN22H?g93 z-A{Gn%F2oknLfR`a@oD2zsHwHTZo#rp45`IJ|1*QOEgIJd>X$d};? zS>(D`UGlJT>Fy7QVIv@L5$c&nQ?`QTOYI3tcEY@bRm&!%s8_&}u$|*}pUX!1IhxZ5 z;gtgHobClhc7dsDpbxCwM1?GVotRgVq_z#02RlbR%@K5yUr zE3;vel^QZvuc4e+;hu_0X2RQ7@vSC2Ts=P7yZ|zLK~4USCO|VW*b-9|s|+GHGTXtO zrHu+O#@NlPIk<3`E5&THdPXsfHLdUR@KVHn0oM$}y%Axy#HCI7GGaIc3S_h}S){-* zZNmK;R&CDEoxAOm`lzJPp#9{Fc7T#G+>=0)n z%TBo8#FprTJH8P-48KLEETn@!1~^k#cbTc1pbtX28B_+F@lQC?xOC7H zIU?7~b#F%qgY{-!hsV$N#^N!3wfjbl^>Qc|5H0m#H4d(5DqXYgA+gNbOHJ$QcNylFM7r9Cs8ggYKixH zr*N6nQg$8i{W_EAh;hiD>Mq4onB!FMbj&cv(=@kl*k0$65fn8=_sppG;R*ip$cXaM zSX7`W31hMD3dX{YLdqMxX&6?v|1|e!yo3GrIvz&;0G9A$BnHd|Xve3xh1h0J-;VOOrHu)}E(O0FuOoI5 z6yls_6V^2H1@P0$387?>=zM_}O`_z@qOCrHbUuRe z!Mqr*5GI8N70_{v+_Ne8{irQ)A$)sXVW|>+-jxm6ltJWv;YZtU-R~d@h`hRv06Jb{ zU#^UJW-J}dYc8(C)G6z%UC)|+_Ubcqn~fD2Ap}^tn-1`Bs5~MPqal37keP>mS*GU5 zUJ%luFbX}r(frojk16}^-Wk`>3j}VI>V528!ex{c9_&}MDUu>|BV0F~bQOt-krP`B zD4P}FEJb)}qTdRp_VsRl^1;1K!oSg;hJlE2>0YcTj;-Mv{1B5FxA2AJM5RFwv}I0L zIXs2}5`~4_5;jiENs9Mdjf0@o^u5OFod@>^T0s7eKCK0FD8Ouu3g8td4+uut6b2;F z-C12myAJP@oEPfAKXGa;xwRF-hghL(Y^5v};ybp$dP-(ITnHKfRn^sZoxnpphJS# znt{V4blbZKY?}qPnTUkR<8YMlM7QJ6eGkhN8686_=J> zveoCv64p&@9j~UCp%LQc9f3z}R4UW7BNOOFW6?0tWU_w}xI0^#g4+wd{2-vy4?w!i(r^e0 za4zAkFpY|AefiTHP=<;J!V3`{iR8?Q#+umLev|Pb16G7{NnXDMuY`5PRdu|>)SLRP zO}~2SXQ1?DMc?V^hNtgBI{MS$l@2R3tZC~=oK|%Qw-UXUla`(2Z=}0*6H7P33rOhQ ffc}6`Ge2E;wNS{96&9vv@?-h8@-}Ke$xr_eGvc+2 literal 0 HcmV?d00001 diff --git a/analysis.py b/analysis.py new file mode 100644 index 0000000..620a15b --- /dev/null +++ b/analysis.py @@ -0,0 +1,124 @@ +import os +import sys +import traceback +import datetime as dt +import subprocess as sbp +import msgpi.presg as psg +import msgpi.sg as sgm +import msgpi.io.iosc as miosc +import msgpi.io.iovabs as miovabs + + +def solve(sg_xml, analysis, solver, scrnout=True): + """ + + Parameters + ---------- + sg_xml : str + File name of SG design parameters (XML format) + analysis : str + Analysis to be carried out + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + Format of the generated input file ('vabs' or 'swiftcomp') + """ + + # Preprocess + sg_in, smdim = psg.preSG(sg_xml, analysis, solver) + + # Solve + run(sg_in, analysis, solver, smdim, scrnout) + + # Parse results + print(' - reading results...') + if analysis == 'h': + if solver == 'vabs': + sm = miovabs.readVABSOutHomo(sg_in + '.k') + elif solver == 'swiftcomp': + sm = miosc.readSCOutHomo(sg_in + '.k', smdim) + return sm + elif analysis == 'd': + pass + elif 'f' in analysis: + results = miosc.readSCOutFailure(sg_in, analysis) + return results + + +def run(input_name, analysis, solver, smdim, scrnout=True): + """ Run codes + + Parameters + ---------- + solver : str + Excect command string of the code ('vabs' or 'swiftcomp') + input_name : str + Name of the input file + analysis : str + Analysis to be carried out + h - homogenization + d or l - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + smdim : int + Dimension of the macroscopic structural model + + :param scrnout: Print solver messages + :type scrnout: bool + """ + try: + analysis_long = { + 'h': 'homogenization', + 'ha': 'homogenization aperiodic', + 'd': 'recover', + 'l': 'dehomogenization', + 'la': 'dehomogenization aperiodic', + 'lg': 'dehomogenization gmsh format', + 'lag': 'dehomogenization aperiodic gmsh format', + 'f': 'initial failure strength', + 'fe': 'initial failure envelope', + 'fi': 'initial failure indices strength ratios' + } + msg = ' - running {cn} for {an}...\n' + + if solver == 'vabs': + solver = solver + 'iii' + cmd = [solver, input_name] + + if solver == 'swiftcomp': + cmd.append(str(smdim) + 'D') + if 'd' in analysis: + analysis = analysis.replace('d', 'l') + cmd.append(analysis.upper()) + + cmd = ' '.join(cmd) + + # try: + if scrnout: + sys.stdout.write(msg.format( + cn=solver, an=analysis_long[analysis] + )) + sys.stdout.flush() + sbp.call(cmd) + else: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + # except: + # sys.stdout.write('failed\n') + # sys.stdout.flush() + # return + # else: + # sys.stdout.write('done\n') + # sys.stdout.flush() + except: + e = traceback.format_exc() + print(e) + + +def runBatch(): + + return diff --git a/analysis.pyc b/analysis.pyc new file mode 100644 index 0000000000000000000000000000000000000000..18137af8c4d02f89c8ec6c366831aaea45f7a29c GIT binary patch literal 3625 zcmdT`OK%)S5U!cE*UrX{6O-74R~jVnB4Gm&ijaj6;NTz>44O?KvQ|Xn>E7Mh%wsj( zlVD<>%!LDIF8mWD4qQ0%TlfL+Rn5#A8-tL51G78bJ>B){>gunm^~+)>|M|_&BdQ-C z{X2N=w`er}9wkJ()3HapUNi60ZmSv3(C%zAZqs;CV#cy^*04CnE-@owR@KcaC& zV~-Aq4$%5^)S|vmhck3IOU|c*nR@Uoo27BPp@7sQXaf=Sx2eZs+d}vqWiODPW9<*xRQw=(_!N7X z2je;F>3MAN9Pi~PkEqj8cs>b>Fg4B;RzBCKk9^wU{V&o)Z7S;unP8UfI58>%m&ymK zcS~tw<6)*IQ%S0O@1xslkQZs_UKX}63zK~t+f$O>ZKhqW_Kb?cqKJ)FdFiIgMrvJ+ z@-!csEdDZdah^?=G?uhEx%g_5M`2P^zFL?le`t;^AAlpyVi$*r8ia9D7KV3bvY{JI z#SY9FR3>|9l6+#0sp2!J;!MYpu`0~;RIV!6>dLz|Phid|Cv1zK*q3f(6!vMr0WZYa zq;zV)&hD-~4EOAs$_uq-AH@R~<>_Rta;<$52g=8?t_l)8&vyToO&}5YdOeFuhwcQwhBu!u(e!bYR} zQI4?gr*=4r`;9#tll@8~n{=KPKEpbmc)GRHn)Md&TfuwPn?c*cyX-~mLH0KLO;7~z zL2-xC29R(AP=LFCLj&R3=;jDl^oBl1x`PhDwt#Mt^aXSPwTtL3kzPgzG`mbU06Y2$ zBLbt!ia@~w8L!gnp)WWAR3aKy5CEblh3IF7gE0;oC?wxf}-~S>qk%N#ltuHuVw))<(=)b(Da-HuBY@USEof=8`_8;GQ z0h%W8P9erWg-yUn0=W2}VJBa6I0CnR3c@$z)MOSpp*e*#j0%8HK8Yg+;4*S$0jy7P zI7xNxo06tf4JDiElL7)&4a}*rHXND?tT)~M1nHr|oK3Zt=ZPqG7Z6h?vrd2+5+FV~ zufmBb;v7jq@M#p*_#{~=V3c(Xil+}kusE?#ow0OEViF1jofe>msU6h`STrX?`_f-Z zJHf$*CE-vu^k?ZZ;GA+MMUyH;@xa(M5T2IOB+4fF)=9EriS0G8)G4jZGQMSxqwJ73 zq#lRaArf0B;LIjzDURcS(Orn zMJDX3k7|;HkrC7{v0O+h9p{J{c9Pw<0sB?f9~H4PBB-i-WO9&{c2v>sO@tZcNg@MO z6P1Cv6m1@ay<(f5BVy7SOIDIbQKr{3XBTwUH(- z`-|R1f7w%Bs}->Xe)xd5@LG;E4WD#*X6JYV$D@BU!OuepDxLR27mZ3*&#^^*cNJ~$ z1~=JMiX26t4PIwDE)%4R*gzq5TvFL#z@>qRsg|)OHB5D8 zc9=>9RjtMOhSXF{5a$v-{2g#;NE!A!GAYeChN>zp;Qa{4ouV|s6@1BQ1#e+=>_Hj* g+irXdxLYRXZ9X@984Z_6{9D3r*;{TctSxN)2Cm_{FKzg`kf;g`kf;I zQ<|4qk{DB*n3EGzY+$5YTvD1+8B=VeTa;Rynpl*a5mQv0oDrX#pORV}lUtl#kQtMi g9}^#+3Dg!JuUAlci^C>2KczG$)edCEXCP((01=2N-v9sr literal 0 HcmV?d00001 diff --git a/io/__pycache__/iosc.cpython-37.pyc b/io/__pycache__/iosc.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..984207085acddd2d605e5e6ed4b39dd8dd453f22 GIT binary patch literal 11875 zcmcIqYit`=cAh6Wq)3W-STD==#BmbKjaA1<(^PO0uOHbsO_VfF)+ymcp?QZAC2~mJ z8QBsubWvyiD1zPk+QKbKh#i9s`qFoe45fsIs{|XeVAO#8($d6!w z_B&^WLsFJoZ;DcyJNI$!xvz82{m!{_DWA_O`0c76zccxcqWmi*x<4I}i+H@>M<5iT z)|FE9R!bVE^@gD-s?sn^DK#odH>^?yaZTvoR!Ui62ovuikrHXVbHWlCyz?R}hVU+k zoXF!nEDB;6@1iJ*5xhsls2Ib0RE&!wc#nw*F^TuMI4Y*a%(vCj5usjJ=4SWNGi9-& z(F?zmcrth{tIK;6?jrqDVmp?A`F zEIb)JS>z4j$pxy=*EIQtP&+x7QU+4*3Zs+5llN1BBFt@VC-tGyDeRb=U;62wu&rKE zcFebbigt%P#h|!r?J4vYsWr{Ji+374ez<-&K&!v3bW)(mFhwfJ1^G#3P2JI?9kslp z%1cCt8ejNXf5^`T#cgd{Wq#VU5|~SQ_YF{>U!Xn7G=+O;hS91H>Y^~Fm1PaJitcgL zQn}Xm4%R~Xi2HIhT5$bd)C0s?_w565i~4XFw_*^FP3Th!xB1Wa;WQfe)La^idW*h{ zBc}3LkI<~ev!uzp_u(-S@#u_196H1k^*hAT!h3xEF8HK_5uWK%cl1IRncu>k()LkszptO2H2_PLEydm4Zk8^!UE&670-^|rk|(5kqv;&!8?n_Lac45gwL~t-ol0`CUz6PFU^18rrk8c}?!-NP$OPlu=3n)- zIm2Y zv50HuIQI8Ldo*(=gQMViEQ(Dls}GYLhywD4xSgT(yYABzXYTKI#~ST0|I~qc!~OMw zV{xDU8udx?Q}QB@izTLWg5%M6?cl@#PDS)$1l;UH@<%TyVyuIn?Wk@Wtx(@m@~OBL zQ@%LR%M-x~lF>u>jD!CP-t{~W$>^zVlDZSL0*BNiS?kTt=ZKS*Fy)U2+I=W^=sxXE zqW(koX?Kdt4sUmQpk0Y}t$0>R?{#&W9oYA2m!!wOPrD=q_I<`X+nsOd_ha3agnk$j zW6-*H=_ClUEG=1*%JsXwRfcYQWc}_9<%S|MzoxvdyoLG1*>tRZ*M8k83tN=^vR#$U zhV8F7_Vqips(-oJXxTNlwc*>l-Ua((4zU`J9j@_)@zarQ3 zw^~lOn&-=!yUhId>HL-|f_)C?kPKBBs@YK4H;`DYP(J;rh;pvu@um{r-j=5xlpZ1I&V4n>gwf%n{RrK^lomrHNSk*E7$8cy>sW! zdj5vky6K%iD;>`%%gV}4=~Y&iDox>dHyhq^t9G;2#J^XWZ*6^*y-=@teg#wXoS^G;I%7g|@N0=6l>(%_NFjHx|o?mu-FI1gS zb2s)YM5Q@T;5m=y$xvT@;yeWuJO5;8xEn}Q$no<}&OI7xO)oTBWq&2KyjH#DJDYxJ z-70$ytx%|Y%Z)I9-Ip5`e?vOobX+Gink~l-Gmg8_aAetcLX(!Jg%{Ni4eENRH_EM0 zuep9`HZU3W(5%**Wk1Zh z!dZ6QC8zE*9B|D0vg~?Ea>8x8Se9y2HkN#+(V|WfBhy0JFu3Zyem@2n0 z1QA+so)_kNrRPPMX~|}b`G>}uv*ls1XnKiq!y%3~WsQ1bHJVt>TbRP(S1=P21H0~2 z{8|(9(}**ucD2fty^=wT7n;q6vs|7#iAkX?fu1&~r)i{U<~f3OaIsXh2pW-8nj5ZP z%JOT8hE&Q&@wn%u;V4-zZ#j|#rBsMAy8S5S>owO|N|IT$Aj?ZDkiF30P_Ne9l2)(L z61q+)g9YZ6yi$&$aYJ4yPpK{+uT-F9(ySNiPEC}|`j+%tp;`{rTcKJh8FjxVAQdP| zL0S0aCrf67yJ=B%S~occ@s;Bw!|HNx^T{0LSlvA`jFdM++dExUb80~|)c!w9O{o^= z6x0dCimIkrco(#+YH0=3$Rfq3aZJ;wlxU~$XJ8CSqy1U>GrE+ArD*I~ zjmBtzqox}AEZ{t5C$SugNFDM@|nq8DGSBJbTr+<`SwGIFRUta9e0>-C z9TrhK>LU>6rWl45qy_q}2@IM7=D$Lt0-nwIsoiv-0ijuwgwrU*Lzr++4brPXn?M#Z zvcVkmIJ^%Igov{-hb-5P!h1Q4M%W0m)#LCfIGBg@VeHWMfqXEx-P)SOXhiCe_V2At zZBy-ZU+s&iokr~;KSKxzGXP{VsQ`)R_LSZHUSuzH$4|(L+RIOa0#cHSsKf+jIfdFO zF&Wha{_oZSD%#QAA4IDSl>2=#-LbZw0hb4xci1%DdB35!#|_0VtYxv5|DwV|R~3I) z%xuF3`E|wB-cp$FgRKVfjaQ?OuP&cT_#W-!>!ZIJj5)Vi#SFmShkWEjJjd8Ki{07a zcRrdC%*SjzlSSCT`N0|L6oLZJ67}R5p5s8Qkm4PE{aYl*$yjyyH%XsP9OzT_|G7_g z(x=Y{`eewPNuM6-_37u54jMN}2U!BV!(j;&SvsD2VCfhFkK%>?*)rspNzeKvgLpoe zWDrl-DSazo%9j$J{Zi2a&m2}@c=HxUKFm@vEKWxI7ET3>?<|EJdYZU4Fr%f{w&vR29tf9 z;wk`|oGDG5cJ1SDz#O?+tr8Tw4a3B*RjaP!dG;%1Al1$GiNSKRS?aa2TXA9nN%lp1 z{Eb?p>8&{KTC>NZeWWk9SFwHMIv|v^F9Q`i(ta5*o@Pyw>Cqm$=AM16UUnV(H6Tc! zvzp_z$6lyAn;Q~C>SanY2-?%!J9ilv42aTQwl8;aetV=hs4I~buqZ9m2{T@Ij4%u8 zWzVa@jtGGkyCd6Yvb1!7@f;BjVO%g*xVK-a&+h*hesTny_5l!-_Nm3D4bU3TU&7l4 z=)?TkRbWvu-+p9(s@K*a8JksdlzN$lVMz1PK8Bh#_bexDg4NNa#Ee6$8qW*Q@)6L+l90g}6 zI7`7}6r4s7TEFUSU6ryaL$ieyf{jye-f`qyG1RKAyb6+1h9I;@TdI_Y32~<>g?G&r z&Zc~sinCO_xY4-fNcjq})VxMXcN##`HKIyHwlvJ?o}`5qXDoS@(Bv2zFW({q>vp-m z0jnUUk4(dZFghQ>=vO!V7qPwn4lR41BTzEHkwB}XK(i6v9SnQ9z1W17HL^gR5xPBX z48Bjr>Evw^4n7W)YiJr{-^YQQ3-`olP)3NlsM<;S!MuZ+T{w|D)0Eig=xdbBvD;tcHHQ?meJLcv~z{w`#MV%b0^Ti{h33!XTTc>V_W4`aAL$UAjWNVS?AV%hI4We+%Ph5`_ASi=C<(l zxb3LxFWpcA^2p>C+xmR9#yUMztL^D)k?+8Mpt%CWISs9-)wSr153V5|V-3dKkkqkxUxVzOJ~5!WH2-g6*Qcngv7 z7jIe^fMnoeYiR|I-|7s?YzRG>r(|5tci);^2ZsvH75&J;L|Xqc9`8H630AQUcrh0P zTr=U(2&w}v5MYN%wTp`;m{!Ifj4J?eb2k+b6$BE-lBNRCZ9Ye}AXGo=n7+kaVBl<% zpd0dcyJ?1F55RAL16@&lmT(}m>F4E;S~DQ1Qwa|;o*sZlLQj`14?w^EGB(E3n@^!$ zny$bw8;?S17=VO_&tNK=yY&(ruJBmSr572AL~=^YwNLPdPke?75Zhm(j9)>}@4|XQ zgytOyUtX9ZuNV&PGEt0h-V2=b0(s*+`84J0b@=m4kT<-!w7f=HMzwis72dyUje=^Z zyEP#fQBq2|l_n0egQA*a-q$ZLlC$k&64P^_Rv>oc5WrE0A^BD))_sphX%45fh)Q-) zr$RfDegs((!Gy7=!17X7X&P_Bd!a%oRUbCwZ(~XW9rvB&i%?b(2WB=%^KL%so5+D4 zCecH^@z}SK3XH%M=A;7ceDMafexKAc(ner>rx1jwBrw8>Bx*vV~EKWr|{_H@HTj=p$ zqOBqH_|M}W5B)Da&K>A+p0>Mpd!Ww+Rsn5zOwu^_250h_q}Ta16PxK1^m~};Kk4r| z?i~8&8oK?1!M-EsWGw8(MOnat p)$VljE=SZu+i^7tzzEKt|`|@q1LmhZ@e}YR^ zU~EKNgIt6q({~{Db2;e|o#w(n-F9TtW8E#EMP>vDY41eaW2|p`Tmy_>hQD7fQ1KmN z(dQz5Qp{PrPvjOAw<(~_5gOzpjg)$TbSd&1^JosJ_ThNHT#bC%-j8U@Ttpq(G6bax za0CJd0hC24Lp=p;Pk9A69sBc7bN)E;X>-Zw@WAhmZ3O+uocPF13VaGSC?EwRZ&UD_2*T6~ZYC;g zbE$*V^d6B!iX`%(A5uFLvC(j{A%vpXi$G0x#KxkeN$g4EoFSj5jKK*we59dXAzwJ)!L8=O zfp#U}l~)1oafR^jG&*XEM)$ZXKR_%Q8{PjTV=I9YC+zFJmE2L^z7&toB(;e12`Y=6 zpBP(;#wPTtPAiFAtSumFiX4trjq+j&I#=6>(`lHF=&-2)^bo>Pv&^*a-*@?e`ty6J z;UFUg#v#BW0AkW6^5G~!2-hLP!ab$KoDT;pN=Le-?mKaNqpVzB++6fiVvLpfCBWo! zfTb4FgfvceGLU3)ZCa=^4nF90fF4Cs4OKG=U-AKwd9)=&{EH3qM(ILaos%>5W!}RZys|9+bIayt!cI@@19xpCojT zqPC(=rLwrY7>D|?7&=<&ar_a0q;vZuDIQgR2+Bk#35*aV zA$sTIH@=QSYMAqi+v;Ab+eXxyKN6rt@{w*EQQL9;Xq3OL_veBWtUrA1qDtQ%yH8DO z9W8;5KDH_MD{9!)1k#sy%X&ZmTcCNC!gOl|_r4WRo~NuQDR_zk2D@|^Gj3(r%ljwf zD@dTJAeuDLH3Fh1{P zR8er+j%yRJsAjcUghdtZ5WT39MMf8L6pLIaaFCz`SB^SeGSX>IbQID<9OUnT-@v9y zG(G8s?xy-sz!gC!HEk?hKM8z@$d^+P8RWwnLq1>8_42WJxX9zJ!~(s?i<@xC401WJ zz~-XdLHm*3x=ArsS+S|iii&mA_;S5@3pe_rh7SPB6~3eQN|QWm{l-S*S__LYxorB; zvECDba3n5+aU#4W!cAeWKkI9)+Lc#m)#wBzY(*EWKn_1ZDw%UaNXeZ47>Q)g{}flF zw=u7ncScMrIAw7;%4aLi&`71=45Ant8C6zq!tAaLRsbiGNI~J@R00*w(|8s1-3*=1 zQBBa`gVJuT^(k_J&~;>7Q(4Up)B@NatZ`uNtZ41v3??&Bk}H}#um{!>#VhzSX9@Nn zfbt_m@C^$1IN}ASlL#M%#C{V7!oB5+@ZR^S#vdSvZAW?7mSR80cZu@8o9l%)X_>}B zgJps}Xs8+%s-VrNwmxCXkB}d#D~WMvB25n;K6L((n~sdaZOy~ykk{d>=-|4i14A&- z;Ep@z>rj9y8H%KO-Q&1`MZR|tSE*FK=wtu%$^-fcH{JdV`vy&*ef^t#A2oUKAIaXm z$H0w9*YUGpzX9o`3yf~gJYP`Q_?@*MTX}53CJ$`yf~CJyJZ7VT$Lb)u;tbxy^pjOl z-gSj_TjBD>w)x(|GosVB-dlL$>G`LxcH|%O^odZd&zWJi=GE9gTyY{__$idK8sdF1 zu*;sHqMm|0uzZ65QGV>3%FU@s5{T8NO>m_koq1~b3@i|D9_?KX`Q!h>$H`&Qr1!HdCRhn{5PP# B(9r+@ literal 0 HcmV?d00001 diff --git a/io/__pycache__/iovabs.cpython-37.pyc b/io/__pycache__/iovabs.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..2abba6321e7b57c44779dd6d2d6251ea6bd1c88b GIT binary patch literal 9501 zcmbtaO>i8?b)LVSot^!`{~rmG!yiG|ghE70QDlapB$JXXI)W6D4iy9pYPdbSTwr&0 z*)vN4>)DtIXh)TDEW~BIoJ+D$9@{>(~8WzwZ9}4Zobv+X{Zqes_Lj{7psqS0eg99hq101m8g-6rt9XwdAd? zsZ48Y8o%{5gWsvTsVS;bPp?^O!kMYtYgyzqp+8pEa$-)Tg!x!q%Zs$IP*M;XVdFg} zvLc7~xX6nF-VjN?5iCd4G(Q(|7sh*`9p7V0gfxNv|mDXS%wUie+WlfiQxPw)#< zQTEK8pM=V$+SMXuUtQNaYDZ&QkJNp&qlZQ}73o41+P>Z~gx*QvG4Z5@*0Jzp@Yr~= zD9Z_>lMl^E5vg5m&;00Kr?6*kzaOTf!mfG)ecio@9>zN3(RkO|SLiLwEiI_wo!-*q z-h)~Mxxb}!(%{H22WZs%G_=EPn2W}DwOy6l=?h9{0_`-sCm+Q^!QNc@<3Z_JM*b=pN#gyc^+?PemLUJg8zf3aGe<8I>~v=3&yhmU!C8+g0akY=k_tiyBFb? zi=9)^iRcvR4V&uHj81|=e?4lkGykr#r#~?7DA7D-gm~7&d6Ke=^)aV(=F$EnTW~g- zhptq&&^@uQY^A~butM4|a1Fl~lz}CP^Ca1a*0;!2bfn1fo9T)Ujz!`KPjyd63;XaY zqRG$b;qo(2l+*l;S~;<+AzzK`XcD$YZp3z{M=X5Eb|ZS7hD5A)%$v>gmW^h%RIHvH zTf~svzPxTmp(*q*UYGv{x8rk?Ovu-4)uF3I`@=*X)T?5O22bl zcYMF(mjmni5_#I^$f~r%_Vqm*oqbBb3uu4#DgDxHo=tYvC-gnpx3`~zwG&+*VSaOB z68n2?WKPdV=a1NO5p8k0L|OKVEcxAo+8yPNBJAH&-c#O3&RQ`9kfreA}fR18G$pHa} zSe3DA$I9S@CLRONEj+;@66_kv-RX`R0ZVM@(g>AE-Btfmdq6s=vHCm`z3PWWE!I|^ zkJajr6q!T2Pv5+LW&LY`Cxi8tUk%;$z^&ERgXb?l7lbXbvmRW2PI`go%F_LH8I!1YISo1F7xs2z9Sg%~VObO*K zzYrUK3pAyiM@d?Eeozf}ic7K92x6n@hWBGDXx6Hsw;jgTJvZ=ZzG6M7)Z_fEP_{~8 zOL}j5z84#frsu~Q&u`T|>4sjMqKRtaMeSpQ)QI)E+mr;C;#3{Rsl}G-yS1I58pP@5 zeY{ISoGRBEZWyP0vR-U^HLvdZVO0dN={C`ah^<~(5a))qmr-7Cz^3;w^b|&14`RJW zLu$%KbAuRTW6RqKFh=Azs5H)UwjrzN0Y+=adI_!dpz6nCZ$k-Lb!)f0Qdn(3PQ6z^ z63S&3fNUt;l6TK7ZN+xEA;GiZZI|Ly6HXcyFQ9PX@o^e?vLR<+Q6JVPb810l%`Vs6 z%0{UrA9`y!&NuGWst;PQW2w=?=xly%kX7P*lJ6-Wk0oTy-SH$lLm?^X4>-=(s=l`| zAdUO*vsf=fORZL=sqwwFEO(ly5En3u8v}XJlIUbYa+K}$9K?B|`hsCOuPWe>ULEsZ z+mT^2HnAgw?j_7TY7*jwSlx=%hv*lx3P+$J%SkweoFdOqE3rY2$ZIZ1@>FgQSH}eC z;6+;UX+s@Xb80~|RD+(8zk)ieTBy}D3vWv+P>OmRIeYMz&{R`UO@3_VBhG16#Mo){i}YvkXux?KMRu??(mbb5!!D2oi6-cRCTU~1U7|gQezda< zE-*26&`{~Bg03oBMd|2OrDKG+$_Q;&`{Ecg4nxugmm(|7bZxPS(dmHjr{)xl&D>XLi~#&e>FL9z6&-fWpuPDkw71dT`W4z| z1tfH{0L>;YDd!-~W@-NyJ+y)S&j6Ab0?-Tg$aV95TfidDoI)2A0O|f*flHjds_eXo z77Jgv#ks2lKVCAFu&`xg#Q%^0k;75}4@Y}Jn?5iRNN%aa?m=+x7xyxgnj~1ZbdK(pm8-y@jRrJ4tUopi8aV_Jnf_ zXQceRUV(5Kv=Sh6t~VM|RDA#h=hjYu^KVTHO6^Hj^j*|{4S3PH)pY$&)!#YzZ~WwG zvc`$)HJG9bpy%8QtL3uq!3b}=Pf_ANjF2hLP={Sp9@^ngx0(SByN;e#OCjPxieQ$nZo2o@UaW1)9i;a9MZV^jxwNEFj%D3ntt z0XGP09R72T&`GVW%>o)20130|A^^&hIlvR`3Tl?YW2qMbWR}ok8TnWn5t1~p8|u#jz*1<#T&4TS4S>$Tf)E9GY8-#^L*z@e!O@RuNJ&=k?~vL4m=pw9 z`XM5}>ZT^tP#4-?Kv5i=i2eeMK;_T?7mi~B5+a58?af0ai0H~ykp}p>L9yUIC!}8P zK`;}oQu1+>;;8RqL^2J;IITGaQ09+%R3@SBC}Fy5!Jq~k)an`3UOk|IJg&Ye!GN65 zb!IZgRU37YZ{LwOI6>Z9;^z(X!7L?c;+!x z1i{BO`CZhN=rE!mO&IM-JYn0&Y|9ak6fke9V-7qrm3ZWURyK_&MTm#}G}Sd>bi$}W zLA@NPeGM3yd|HQ311A$N(E`|T#G{SX^guV;0aV@pfzuY)1J5lpI>r zgA2qB)NEfHG2!#0i#kF~H`3?uMQ2W-1rC2DQz1qVGd*DrG6KhS5@$2yi&Q}CeQ1u} zk{|}=5~rl=V{>{$1+?Od_R&p>TulBUl31DRz5j*bf-B9DwB%^Sh!xYD##hLSO+G~HgQ zgmFN2HJMphG&z%a-k*s5QeY_{&*9i_p+n|^r!FW-)R*TBM_;z?mSkfCU_k^$P14%u zX)UUW3-I1a!6-?a;`eA8$;}YDT^gBXh`51(_;DfJqLFwv0V27HeZ6 zSQfPU2K)4jm=#3z6+&FQS};zoBJk}UU)2qSukZ~GU(>aZOhR7RgHw^lxGGIwYuU97 ze33WsMc(LG*y}CAWb%ChDBz!Cnh;ryXgro*;4<2++n}jc1|=D8^=;s`Os^NhZg!MO zc3pv>5|EiyF*-KNv%=cagP$RX??>2+hZIf+qGXJ4noVIf@Z^~%gBJhERG;U+m?wpA zmx*3BX+y0b13F-Le;G9eXpA2FeeoL-l5v5Xq$S+92i$pHge;Y9HHMHE=`-GYPr!D3 z5=E)*q_BAdnu>6?6VBZ=Fp6RMbW*;nA1a;cb4S{4YP2TOdZYPeINRqrB605ca*Pk! zKyOVh;kz=8#yH$7oj<)GAN}u+>}osp^(!xnPTTtWl}j(Jy!3h}RzDbmkm1Uyz~A`S z%Y5<*Z$!m~*skL90ey-odC6t-0{D{4=J2F>iOLwy@KSn?$csqgd=;Sre0A6eKcEjL z6Md4t47l{MB~CYz4+DH!C3hf(+I*@#B0+Q|VVKK6jqwbhA-r~q*<^hZQj_1K9xHz^C52Njfawl_#UzHE;!f<h~uyiEaYQYxQXfs}5%L#q zu9^*tFQKa^FmS?YM2;z+=MHyY(_@P-@G7st%I1sw4zI(?;Y<7jJ_9S4&+@s~AtU7sMuKox)}J- z6-KibaY^&F3cP;Pzkx0kbD=^hIJ!0fmFS`HgFY*+mUQ*-hSCy!FTXSv=TR7o7JJ@f zewvDygL8rfeTYM#g=D9v?fCSR`90j=>+$)_kD{p!#lU8q=#s5#o!umfjQm&o<^9$h zI>~AMz_W4}^@zMucs%WzcU+74&}$5{OTGpw7~;9VYFe4kOh_$N93czi(ppe2N7G zph7ETeumcR3Uy$Vb7cOt@BXg{f7F)h8vNu2Qh@Y6lLMz4{3D0i=j&H=^|v zL}z3|M_UC=^&4aMqT#wF42fUeSr2(F->-#oZJ}%FxZz9bkEaej8~KAY60KspySod;FBXb@ z!A{mTS(80j+hA*3>>=C6#Ur-E9d=a*jEs#EhzJ2STM~nMg literal 0 HcmV?d00001 diff --git a/io/iosc.py b/io/iosc.py new file mode 100644 index 0000000..9996dc7 --- /dev/null +++ b/io/iosc.py @@ -0,0 +1,639 @@ +import os +import numpy as np +import msgpi.io.utils as utl +import msgpi.sg as sgm + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readSCIn(fn_sg, smdim): + """ Read data from the SwiftComp input file + + :param fn_sg: File name of the SwiftComp input file + :type fn_sg: string + """ + + fn_base, fn_extn = os.path.splitext(fn_sg) + name = os.path.basename(fn_base) + sg = sgm.StructureGene(name, 3, smdim) + + count = 0 + count_node = 1 + count_element = 1 + count_layertype = 1 + count_material = 1 + line_material = 0 + extra_head = 0 + if smdim == 1: + extra_head = 3 + elif smdim == 2: + extra_head = 2 + head = 2 # at least 3 lines in the head (flag) part + with open(fn_sg, 'r') as fin: + for li, line in enumerate(fin): + line = line.strip() + if line == '\n' or line == '': + continue + + count += 1 + line = line.split() + if count <= (extra_head + head): + # Read head + if smdim == 1: + if count == 1: + # model + line = list(map(int, line)) + sg.model = line[0] + elif count == 2: + # initial twist/curvatures + line = list(map(float, line)) + sg.initial_twist = line[0] + sg.initial_curvature = [line[1], line[2]] + elif count == 3: + # oblique + line = list(map(float, line)) + sg.oblique = line + elif smdim == 2: + if count == 1: + # model + line = list(map(int, line)) + sg.model = line[0] + elif count == 2: + # initial twist/curvature + line = list(map(float, line)) + sg.initial_curvature = line + + if count == (extra_head + head - 1): + # analysis elem_flag trans_flag temp_flag + line = list(map(int, line)) + sg.analysis = line[0] + sg.degen_element = line[1] + sg.trans_element = line[2] + sg.nonuniform_temperature = line[3] + elif count == (extra_head + head): + # nsg nnode nelem nmate nslave nlayer + line = list(map(int, line)) + sg.sgdim = line[0] + num_nodes = line[1] + num_elements = line[2] + num_materials = line[3] + num_layertypes = line[4] + continue + elif count_node <= num_nodes: + # Read nodal coordinates + sg.nodes[int(line[0])] = list(map(float, line[1:])) + count_node += 1 + continue + elif count_element <= num_elements: + # Read element material/layer connectivities + line = list(map(int, line)) + eid = line[0] + lyrtp = line[1] + + sg.elementids.append(eid) + sg.elements[eid] = [i for i in line[2:] if i != 0] + sg.elementids2d.append(eid) + + if lyrtp not in sg.prop_elem.keys(): + sg.prop_elem[lyrtp] = [] + sg.prop_elem[lyrtp].append(eid) + sg.elem_prop[eid] = lyrtp + + count_element += 1 + continue + elif (sg.trans_element != 0) and (count_element <= (2 * num_elements)): + # Read element orientation (a b c) + eid = int(line[0]) + a = list(map(float, line[1:4])) + b = list(map(float, line[4:7])) + c = list(map(float, line[7:10])) + sg.elem_orient[eid] = [a, b, c] + + count_element += 1 + continue + elif count_layertype <= num_layertypes: + # Read layer types + ltid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + sg.mocombos[ltid] = [mid, theta3] + count_layertype += 1 + continue + elif count_material <= num_materials: + # Read materials + if line_material == 0: + line = list(map(int, line)) + mid = line[0] + mtype = line[1] + num_temp = line[2] + sg.materials[mid] = sgm.MaterialSection() + sg.materials[mid].eff_props[3]['type'] = mtype + line_material += 1 + continue + elif line_material == 1: + line = list(map(float, line)) + sg.materials[mid].eff_props[3]['density'] = dens + line_material += 1 + continue + + if mtype == 0: + # isotropic + if line_material == 2: + e = float(line[0]) + nu = float(line[1]) + sg.materials[mid].eff_props[3]['constants']['e'] = e + sg.materials[mid].eff_props[3]['constants']['nu'] = nu + line_material = 99 + continue + elif mtype == 1: + # orthotropic + if line_material == 2: + # pe = [] + e = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['e1'] = e[0] + sg.materials[mid].eff_props[3]['constants']['e2'] = e[1] + sg.materials[mid].eff_props[3]['constants']['e3'] = e[2] + # for ei in e: + # pe.append(ei) + # sg.materials[mid].property_elastic += e + line_material += 1 + continue + elif line_material == 3: + g = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['g12'] = g[0] + sg.materials[mid].eff_props[3]['constants']['g13'] = g[1] + sg.materials[mid].eff_props[3]['constants']['g23'] = g[2] + # for gi in g: + # pe.append(gi) + # self.materials[mid].property_elastic += g + line_material += 1 + continue + elif line_material == 4: + nu = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['nu12'] = nu[0] + sg.materials[mid].eff_props[3]['constants']['nu13'] = nu[1] + sg.materials[mid].eff_props[3]['constants']['nu23'] = nu[2] + # for nui in nu: + # pe.append(nui) + # sg.materials[mid].property_elastic.append(pe) + line_material = 99 + continue + elif mtype == 2: + # anisotropic + continue + + if line_material == 99: + line_material = 0 + count_material += 1 # next material + continue + continue + else: + # omega + sg.omega = float(line[0]) + + return sg + + +def readSCOutHomo(fn, smdim): + """Read SwiftComp homogenization results. + + :param fn: SwiftComp output file (e.g. example.sg.k) + :type fn: string + + :param smdim: Dimension of the structural model + :type smdim: int + """ + linesRead = [] + keywordsIndex = {} + + with open(fn, 'r') as fin: + for lineNumber, line in enumerate(fin): + linesRead.append(line) + if 'The Effective Stiffness Matrix' in line: + keywordsIndex['The Effective Stiffness Matrix'] = lineNumber + if 'The Effective Compliance Matrix' in line: + keywordsIndex['The Effective Compliance Matrix'] = lineNumber + if smdim == 1: + if 'Timoshenko Stiffness' in line: + keywordsIndex['Timoshenko Stiffness'] = lineNumber + if 'Timoshenko Compliance' in line: + keywordsIndex['Timoshenko Compliance'] = lineNumber + if 'Shear Center Location' in line: + keywordsIndex['Shear Center Location'] = lineNumber + elif smdim == 2: + if 'In-Plane Properties' in line: + keywordsIndex['In-Plane Properties'] = lineNumber + if 'Flexural Properties' in line: + keywordsIndex['Flexural Properties'] = lineNumber + elif smdim == 3: + if 'The Engineering Constants' in line: + keywordsIndex['The Engineering Constants'] = lineNumber + if 'Effective Density' in line: + keywordsIndex['Effective Density'] = lineNumber + + sm = sgm.MaterialSection(smdim) + + # ems = 6 # effective matrix size + + if smdim == 1: + ems = 4 + if 'The Effective Stiffness Matrix' in keywordsIndex.keys(): + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[1]['stiffness']['classical'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + # else: + # print('No effective stiffness matrix found.') + + if 'The Effective Compliance Matrix' in keywordsIndex.keys(): + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[1]['compliance']['classical'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + # else: + # print('No effective compliance matrix found.') + + if 'Timoshenko Stiffness' in keywordsIndex.keys(): + indexStiffness = keywordsIndex['Timoshenko Stiffness'] + sm.eff_props[1]['stiffness']['refined'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 8]) + # else: + # print('No effective timoshenko stiffness matrix found.') + + if 'Timoshenko Compliance' in keywordsIndex.keys(): + indexCompliance = keywordsIndex['Timoshenko Compliance'] + sm.eff_props[1]['compliance']['refined'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 8]) + # else: + # print('No effective timoshenko compliance matrix found.') + + if 'Shear Center Location' in keywordsIndex.keys(): + index_sc = keywordsIndex['Shear Center Location'] + sm.eff_props[1]['shearcenter'] = list(map( + float, linesRead[index_sc+2].strip().split() + )) + # else: + # print('No shear center location found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[1]['density'] = float(line[-1].strip()) + + elif smdim == 2: + ems = 6 + + try: + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[2]['stiffness']['classical'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + except KeyError: + print('No effective stiffness matrix found.') + + try: + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[2]['compliance']['classical'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + except KeyError: + print('No effective compliance matrix found.') + + try: + index = keywordsIndex['In-Plane Properties'] + for line in linesRead[index + 2:index + 8]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[2]['constants']['inplane'][label] = value + except KeyError: + print('No in-plane properties found.') + + try: + index = keywordsIndex['Flexural Properties'] + for line in linesRead[index + 2:index + 8]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[2]['constants']['flexural'][label] = value + except KeyError: + print('No flexural properties found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[2]['density'] = float(line[-1].strip()) + + elif smdim == 3: + ems = 6 + + try: + indexStiffness = keywordsIndex['The Effective Stiffness Matrix'] + sm.eff_props[3]['stiffness'] = utl.textToMatrix( + linesRead[indexStiffness + 2:indexStiffness + 2 + ems]) + except KeyError: + print('No effective stiffness matrix found.') + + try: + indexCompliance = keywordsIndex['The Effective Compliance Matrix'] + sm.eff_props[3]['compliance'] = utl.textToMatrix( + linesRead[indexCompliance + 2:indexCompliance + 2 + ems]) + except KeyError: + print('No effective compliance matrix found.') + + try: + indexConstant = keywordsIndex['The Engineering Constants'] + for line in linesRead[indexConstant + 2:indexConstant + 11]: + line = line.strip() + line = line.split('=') + label = line[0].strip().lower() + value = float(line[-1].strip()) + sm.eff_props[3]['constants'][label] = value + except KeyError: + print('No engineering constants found.') + + line = linesRead[keywordsIndex['Effective Density']] + line = line.strip().split('=') + sm.eff_props[3]['density'] = float(line[-1].strip()) + + return sm + + +def readSCOutFailure(fn_swiftcomp_in, failure_analysis): + fn_swiftcomp_in += '.fi' + lines = [] + kw_index = {} + results = [] + with open(fn_swiftcomp_in, 'r') as fin: + for i, line in enumerate(fin): + lines.append(line) + if failure_analysis == 'f': + # initial failue strength + if 'Initial Failure Strengths' in line: + kw_index['Initial Failure Strengths'] = i + elif failure_analysis == 'fe': + # initial failure envelope + pass + elif failure_analysis == 'fi': + # initial failure indices and strength ratios + pass + + if failure_analysis == 'f': + # initial failure strength + try: + index = kw_index['Initial Failure Strengths'] + for line in lines[index + 2:index + 8]: + line = line.strip().split() + results.append(map(float, line[:2])) + except KeyError: + print('No initial failure strength found.') + elif failure_analysis == 'fe': + # initial failure envelope + for line in lines: + line = line.strip().split() + results.append([int(line[0]), float(line[1]), float(line[2])]) + elif failure_analysis == 'fi': + # initial failure indices and strength ratios + for line in lines: + line = line.strip().split() + results.append([int(line[0]), float(line[1]), float(line[2])]) + + return results + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeSCNodes(sg, fobj, sfi, sff): + if sg.sgdim == 1: + nid = sg.elements[sg.elementids1d[0]][0] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + for eid in sg.elementids1d: + if len(sg.elements[eid]) > 2: + # node 3 + nid = sg.elements[eid][2] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + if len(sg.elements[eid]) == 5: + # node 5 + nid = sg.elements[eid][4] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + if len(sg.elements[eid]) > 3: + # node 4 + nid = sg.elements[eid][3] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + # node 2 + nid = sg.elements[eid][1] + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, sg.nodes[nid]) + else: + for nid, ncoord in sg.nodes.items(): + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, ncoord) + return + + +def writeSCElements(sg, fobj, sfi): + for eid in sg.elementids1d: + elem = np.zeros(7, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] # mid/cid + for i, nid in enumerate(sg.elements[eid]): + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + for eid in sg.elementids2d: + elem = np.zeros(11, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] + elem_type = 'quad' + if (len(sg.elements[eid]) == 3) or (len(sg.elements[eid]) == 6): + elem_type = 'tri' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 3) and (elem_type == 'tri'): + elem[i+3] = nid + else: + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + for eid in sg.elementids3d: + elem = np.zeros(22, dtype=int) + elem[0] = eid + elem[1] = sg.elem_prop[eid] + elem_type = 'hexa' + if (len(sg.elements[eid]) == 4) or (len(sg.elements[eid]) == 10): + elem_type = 'tetra' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 4) and (elem_type == 'tetra'): + elem[i+3] = nid + else: + elem[i+2] = nid + utl.writeFormatIntegers(fobj, elem) + + return + + +def writeSCElementOrientations(sg, fobj, sfi, sff): + for eid, orient in sg.elem_orient.items(): + fobj.write(sfi.format(eid)) + utl.writeFormatFloats(fobj, np.hstack(orient)) + return + + +def writeSCMOCombos(sg, fobj, sfi, sff): + for cid, combo in sg.mocombos.items(): + fobj.write((sfi + sfi + sff + '\n').format(cid, combo[0], combo[1])) + return + + +def writeSCMaterials(sg, fobj, sfi, sff): + for mid, m in sg.materials.items(): + # itype = 0 + # if m.type == 'orthotropic': + # itype = 1 + # elif m.type == 'anisotropic': + # itype = 2 + mp = m.eff_props[3] + utl.writeFormatIntegers(fobj, (mid, mp['type'], 1)) + utl.writeFormatFloats(fobj, (0, mp['density'])) + if mp['type'] == 0: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e'], mpc['nu']]) + elif mp['type'] == 1: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e1'], mpc['e2'], mpc['e3']]) + utl.writeFormatFloats(fobj, [mpc['g12'], mpc['g13'], mpc['g23']]) + utl.writeFormatFloats(fobj, [mpc['nu12'], mpc['nu13'], mpc['nu23']]) + elif mp['type'] == 2: + for i in range(6): + for j in range(i, 6): + fobj.write(sff.format(mp['stiffness'][i][j])) + fobj.write('\n') + fobj.write('\n') + return + + +def writeSCInH(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Extra inputs for dimensionally reducible structures + if (sg.smdim == 1) or (sg.smdim == 2): + # model (0: classical, 1: shear refined) + fobj.write((sfi + '\n').format(sg.model)) + if sg.smdim == 1: # beam + # initial twist/curvatures + fobj.write((sff * 3 + '\n').format(0., 0., 0.)) + # oblique cross section + fobj.write((sff * 2 + '\n').format(1., 0.)) + elif sg.smdim == 2: # shell + # initial twist/curvatures + fobj.write((sff * 2 + '\n').format( + sg.initial_curvature[0], sg.initial_curvature[1] + )) + fobj.write('\n') + + # Head + fobj.write((sfi * 4 + '\n').format( + sg.physics, + sg.degen_element, + sg.trans_element, + sg.nonuniform_temperature + )) + fobj.write((sfi * 6 + '\n').format( + sg.sgdim, + len(sg.nodes), + len(sg.elements), + len(sg.materials), + sg.num_slavenodes, + len(sg.mocombos) + )) + fobj.write('\n') + + # Nodal coordinates + writeSCNodes(sg, fobj, sfi, sff) + fobj.write('\n') + + # Element connectivities + writeSCElements(sg, fobj, sfi) + fobj.write('\n') + + if sg.trans_element != 0: + writeSCElementOrientations(sg, fobj, sfi, sff) + fobj.write('\n') + + # Material-orientation combinations + if len(sg.mocombos) > 0: + writeSCMOCombos(sg, fobj, sfi, sff) + fobj.write('\n') + + # Material + writeSCMaterials(sg, fobj, sfi, sff) + fobj.write('\n') + + # Omega + fobj.write((sff + '\n').format(sg.omega)) + fobj.write('\n') + + fobj.write('\n') + + +def writeSCInD(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + utl.writeFormatFloats(fobj, sg.global_displacements, sff[2:-1]) + utl.writeFormatFloatsMatrix(fobj, sg.global_rotations, sff[2:-1]) + fobj.write((sfi+'\n').format(sg.global_loads_type)) + utl.writeFormatFloats(fobj, sg.global_loads, sff[2:-1]) + + +def writeSCInF(sg, fn, sfi, sff): + with open(fn, 'w') as fobj: + # Addtional material properties + for i, m in sg.materials.items(): + utl.writeFormatIntegers( + fobj, + (m.strength['criterion'], len(m.strength['constants'])), + sfi[2:-1] + ) + fobj.write((sff+'\n').format(m.strength['chara_len'])) + utl.writeFormatFloats(fobj, m.strength['constants'], sff[2:-1]) + + # Load type + fobj.write((sfi+'\n').format(sg.global_loads_type)) + + # (fe) Axes + + # (fi) Loads + utl.writeFormatFloats(fobj, sg.global_loads, sff[2:-1]) + + +def writeSCIn(sg, fn, analysis='h'): + """ Write SG input files for SwiftComp + + :param analysis: Type of analysis. + - 'h': Homogenization + - 'd' or 'l': Dehomogenization + - 'f': Failure analysis + """ + + if not isinstance(sg, sgm.StructureGene): + return + + # string format + sfi = '{:8d}' + sff = '{:16.6E}' + + if 'h' in analysis: + writeSCInH(sg, fn, sfi, sff) + elif ('d' in analysis) or ('l' in analysis): + writeSCInD(sg, fn, sfi, sff) + + return fn diff --git a/io/iosc.pyc b/io/iosc.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d0d4a6c4d5262a36fb911e6394a4478d19e1387f GIT binary patch literal 14796 zcmc&*TWlQHc|J3{TyjY+?-aSZj%3+lY>2j~#HwXBc0pRU8jH4v4k!qe>E+Jsa>(6T za%Lzk?PcOvMuE0L5;RHEx`2Zgf$N79MS%wKLz@>DD2lc~iazuqZv|SQMH93!+M)>1 z^!xrZvlm`u3$+j}_w1Q-`Okm;^WV=X^@ukMk)i!>6Xf% za(bjPq@1i&hLzJRl@aCiNo7Guw}8($kzAY_1j63W@o%BXO+0{aAHH+= zhQHpURHW3pq1H{+Q0k7MmQA&bzf{tmQVo;syJ;24t{Q2=Z5ZlKO6_)WOS{-7>fKWI zNQteY1^YxVHi@#IQ_#3!s^ANRfEIsT8q?DFb<%>3peZUj6=jGL=yJzYU;R<{?trTQ zU3b*28Uw^<{^keW#C1^Z4yncv34BMzb&wFtgzvK*ROFIBQCsg5bVG2x1e7@^o6Nx@ z5X`BG4MPS0W42KJD?`s9jF*h4S4Dj)!q6dNMY71^*kHcVrB-rk>v5A4^p<`%nJ{qJ zgh4?uh1tL~t@33Y1Ko{W^^)7G);7PTfSB-Hpkx>R?3Im*)zd5OPCsbpCLbH)l!B;s%^voGJ-2gT9b?89a!RRdHys>;wqKG4^;xb&M(ZOM z{#B))x$=J`D|f1c=IhA(CV{x_FxEvs zwdaBTBf4Mfm;>gI?x79vBR#avedi#AaAIffg3uTMqW@5YD7olA6d}lI$En8%^@tI` z0cq>ge(E6m>AN5_0f6W~6roAoe<(uvJ*NkaiqlINhI|Qj0?CF1!}o~^j=CZK%_W$B zts-^V!CpA9l?U0qwB#IBFI~s;;?yy^Bws*b*>1tHoI+HvNHU191&geYoKcYaxA*V$g63BSL3!ASg9ni^pTRM1o>LkEf=O9 z00>QAvc2A4mHQ+AK@@n_Nge0!3!ZZLWV(1pv9=A=NC&B{kS;=3{g_ai=<1eAVaE-G z%?Fxo3DPpL;`#1MvqdOe$04n(L7)W%sgT7>U=^YQRw6@lqL2_gCwOSW^Qn0f6fSfIV@Bv%i&`)=fh~**;x)}&IfMj7J}m1 zau62RR*K}p^no-{^`aWiAIfBKEA!EXkA zdDJG2jMt^c-(&`0J0hY{M0=F zdNqJ40njH_`UwNMg1zJ5CUPO{Etm8?k^s-4?-7MR*y;h;YRC)n8%zMefd>f3d{GfUo@DN=gkf(VbQaoeLe(!cUd7kmkQvw*i`5{E(S#%x(!tdlzm;56t z4eh|f{1m7(1?N{|(AqN^J6N-L#b?T*>%$$Mt2K7boiv+XqP z%XF}mhIub6MVN>l+hMK<%qBbihTz$ca~f3j9bqR))J`<;b4J){^}etZQ3E>(E8WR- zu#(x$Z9niyoDo;y89{?|BQiTnyZRR>Brt1*Vr#W(4L-T|T7f>S1^;Hd5{1)y;r%htSyARy$TahQjp z62}xOmINPdkm4zDIHl3XK2&%T#mk5*pDC5-f!{)e8F{6W?}njut$9wQNAIL3iqwHB_q{zk0@C>-zXXib-%ML5HOHHYxP4Xo#DMIoYuo8c%XUiHtv zSSk3f^&$caxZNI5WBZkgTi*__&sMvHtYMz3aZ7#~-ZY{MziiDV;b+KztF_sA9ilEs z94n4@m~3C3ij_hbdc{H|4iA!jdvvoTun1hZ@UA0~4mJ+IRyHEe^)JhXgPp=ykM43v z23n|DZj;!})hu}PAhuOPEaFUu5AGd^kSQT~$4Cj#ZH5kl4Uz(+Rd__xy zNXe-D<(lR+`yz)NWx*VXP2RL8S@1{&w8m`XUY74d(zO`#v&}}7N3$%8@#m{!3gm}O1xkM zbJz|%Xk1d7OpPEi08W@=Y2Vn``g<(3^;07ag+a&j46G%Csb`-<<8}XU+A#gHQU7x; zI)C%|yRL$}A1ncH4U~%4gXL^Jz<0+3?k-Ro2(Xs?lB}W$@Kd9{^vSRCZLESLO~Y(6 zb+r&u3s09kF?*t&OJXUm>U5lSrQlVNB7@Zmpx>gkQ0o+(sVsS7k3{Q-X9%u_{Hx^g ztWrD?mRBamicl1@*u^RG;a-*4RBLUjZk6MXqmUV*Y?#eu#+T}U1(uuIHC}L~nG|F? zMiidtqK(@a2(1-cCvinY=h@c?pPJ6DqPQ=*hxj#hPfXGZBrp6dI+V^9hi zD2IMChujAh#&IZ`GZ19vXmY4Onf#hX8abx*xZOE?LjwvMBE2qI4T{qD4L(9bXy5ShtCNR0Zs0*8EW8&CuwYYZ2`G3HU;(<=Bu zTPvFGc^Rp^LJB9I|#z&U7#%<%jkNkfs+BqnOSQ@%$4Td$}k|B9;5 z2^^>`dG~;7Y#G@nyzn6Bao4j#S_+4>^_K<}=He7tVnJU~pr-FBJo1RXo1oZF(c0S& zv3MB8o?P|?oh7{;;EDj+3j&6;c_Ff26=Zu?-94C_Ww*u#7O5@CQ2((P^IWfA@D}#b|dMFn5>-d8l7MU>kF7 z+rS2pH0n~pi-tCnR8^=qw+BN$yLGIsLH2vh*uUSeHTT=MXTQjWea!v#AGqHfw-x?` zZ|Zas0@*2```-q0njd?=fQJ1(ZF0YZ0s!^S{r0sr$bS3G{ri0*mT{Vu1-WJvX>%AE zw+yVd&5Te)?zXlIj=sDkyNbjiwDK5}NdFo_ZLJ%q76c#!KPZTwV86~baJTaHE#GTy zPz%K>+tX;(Q4O_-)(?zC+w{X#1hHk@eb|q1#1~QQQ8pKyDGHKk3HwW|y~#q<0hI=C zd9hCiU5cpB{r`hnLM!?=t$6P4*%l>dy;{i5a;=?eS!y`=Jf=D|oQeqBo!t#7^rs?!>ud+FV?6 zQ2;|q8v$Uf7BpoaMb&zl=ikE*!ix!iS`bgUbQjM~I_CgCwRU z52}?wG*}^a3gRwzoanNvckRLX2!TWu{ zDn$sl1!A#Z;5|tF2?_iVoT8PE^OVAc2&GwH=D=!W#u$KOm$>Yru@{EBVJzuUj(!Dj zYhK4~GBasQxRNS-K~%d<4c&73B*i9Ku-4tQ&7-y{Kehy5Jbjx+!lqwjGfQF*upswn z6C%psIvzFW>%eV>W8PpP>8V#)`w9!{6R{q6`=ac&w2k&Q=?1nbc+n}0Jtg$4iiN4_ zmqn_ZMbY@uvwmL$g!VI5G)T3{3m4&sl#qwe0SU?Epk9tb+xH?L!}39VrcfJ&;&})j z;INh;>TwcIl+QPC2F#kk;lp7fYw|^t4xdO9R10_o&Sw0mGbFgEm9zx35oAIX5k@z| zvG%FKcapJ4qI9XTcs#^M|HTR`hKDh!))?@w#0X?C|5rLu_@_3ottnht+_ZdLgVbV2 z0zi@oR0V=3FNa&3$Tk<{JUBQg!SE<`$cW%Z&uUHS?UL#%>#?%Y*^k%Ol75h{vy)3edf}+zW!s+_5ndHe z#vTL4cTSEmz*Y($@6XIN`Tge3QlRuy}{XXHZO0=(LN1 z-NpW*sFN?T7vE3(-5aXPSY2H8pBn)A!soe;F=C!HPZ)#dS@SH)G2^5;lo~@mho2mt z@urM?GRCNxON|=2)DYSSp+<&M<5IE>iCUe_;m)|0sFSau(xFaXm$a;Q2(?~;PYAYv zD5mgkG&mS%tn)eKUYOvHPI3n%GFCN-R8n@B`Dq6WL)9^}jg^Sug8glj+CuO@ki^D8 z?1}KmtLWC&U<|c#rM8OK`y3B19Ttl6l3o~*p@gFCX#{=iAvwKxQc#Pu@7H2E$w0*2 zNIp2BgaQ&(w734Exn%wtxr8(Um%y0tBmOF_Fa=W1_}@_JIOBid8P~r~jiCSc4H$d= z=9{hPIrhA8{&*olyR7g%;s>2HfU~7f8FoV0f+sI-2h6|(HL;$mv(R$-@V;>9O`w27 z4nSHWhOl);EyYm4KQR=b*ekl1{L_ShR%q>x*=_to_e1lKzD{J0DZq8D_$qr@N)K7X zYrHGSE9r1l7_8-jyb5^YeMqV|NjYm}0~*HRjzl`|^6?K*w9-=c0SZJ1n{RU*zV}y^ zT=_highE19!k6AsV3Ps!wDDl-M3-pnwPFgr6gz%}X>OrHra>uE3&R{9q0g}hRe%@U zpbfx8I6ZIy?_TEVU1KyH6PyclTy z&$hjG9U{{owC?mfZu^tm{AKH9zFm#W`J`o9F5mg$^?B>u+PTXXZ|GVt5486>=PWF6 zu7W}HZs$7&?ZZmwPi`-gRiu9*Bn@A`{7cTR#F4LGzVP_;^d;1SqyoX}1k~zbUH(5Yz?+^dbb#Iy?W*x&3w5aTekxub4vmXvB zQkmwI2nzjNamn$j>gJ9P_UYrx+ojW9ZF)QMD*9%}tI0rV(5W2Csva1YwY?>$+Dk09 zSiHvKP+Q%klKowd`x6u`iIUnU{C4atUhvqixS!=g%6lBNy(7k8U-r+kzn`7STG^>= PcQ%tfmVGMQlO6v*nXE~o literal 0 HcmV?d00001 diff --git a/io/iovabs.py b/io/iovabs.py new file mode 100644 index 0000000..36bd573 --- /dev/null +++ b/io/iovabs.py @@ -0,0 +1,482 @@ +import os +import sys +import traceback +import numpy as np +import msgpi.io.utils as utl +import msgpi.sg as sgm + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readVABSIn(fn_vabs_in): + """ Read data from the VABS input file + + :param fn_vabs_in: File name of the VABS input file + :type fn_vabs_in: string + """ + + try: + fn_base, fn_extn = os.path.splitext(fn_vabs_in) + name = os.path.basename(fn_base) + sg = sgm.StructureGene(name, 2, 1) + + flag_curve = 0 + flag_oblique = 0 + + count = 0 + count_node = 1 + count_element = 1 + count_layertype = 1 + count_material = 1 + line_material = 0 + head = 3 # at least 3 lines in the head (flag) part + with open(fn_vabs_in, 'r') as fin: + for li, line in enumerate(fin): + line = line.strip() + if line == '\n' or line == '': + continue + + count += 1 + line = line.split() + if count == 1: + # format_flag nlayer + line = list(map(int, line)) + flag_format = line[0] + num_layertypes = line[1] + continue + elif count == 2: + # timoshenko_flag recover_flag thermal_flag + line = list(map(int, line)) + sg.model = line[0] + sg.analysis = line[1] + sg.physics = line[2] + if sg.physics > 0: + sg.physics == 1 + continue + elif count == 3: + # curve_flag oblique_flag trapeze_flag vlasov_flag + line = list(map(int, line)) + flag_curve = line[0] + flag_oblique = line[1] + if line[2] == 1: + sg.model = 3 # trapeze effect + if line[3] == 1: + sg.model = 2 # Vlasov model + if flag_curve == 1: + head += 1 # extra line in the head + if flag_oblique == 1: + head += 1 # extra line in the head + continue + elif flag_curve != 0 and count <= head: + pass + elif flag_oblique != 0 and count <= head: + pass + elif count == (head + 1): + # nnode nelem nmate + line = list(map(int, line)) + num_nodes = line[0] + num_elements = line[1] + num_materials = line[2] + continue + elif count_node <= num_nodes: + # Read nodal coordinates + sg.nodes[int(line[0])] = [ + float(line[1]), float(line[2])] + count_node += 1 + continue + elif count_element <= num_elements: + # Read element connectivities + eid = int(line[0]) + sg.elementids.append(eid) + temp = [int(i) for i in line[1:] if i != '0'] + sg.elements[eid] = temp + sg.elementids2d.append(eid) + count_element += 1 + continue + elif count_element <= (2 * num_elements): + # Read element theta1 + if flag_format == 1: + # new format + eid = int(line[0]) + lyrtp = int(line[1]) + theta1 = float(line[2]) + else: + # old format + eid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + theta1 = float(line[3]) + lyrtp = 0 + for k, v in sg.mocombos.items(): + if (mid == v[0]) and (theta3 == v[1]): + lyrtp = k + break + if lyrtp == 0: + lyrtp = len(sg.mocombos) + 1 + sg.mocombos[lyrtp] = [mid, theta3] + + if lyrtp not in sg.prop_elem.keys(): + sg.prop_elem[lyrtp] = [] + sg.prop_elem[lyrtp].append(eid) + sg.elem_prop[eid] = lyrtp + # self.elem_angle[eid] = theta1 + sg.elem_orient[eid] = [ + [1.0, 0.0, 0.0], + [0.0, 1.0, 0.0], + [0.0, 0.0, 0.0] + ] + sg.elem_orient[eid][1][1] = np.cos(theta1) + sg.elem_orient[eid][1][2] = np.sin(theta1) + count_element += 1 + continue + elif (flag_format == 1) and (count_layertype <= num_layertypes): + # Read layer types + ltid = int(line[0]) + mid = int(line[1]) + theta3 = float(line[2]) + sg.mocombos[ltid] = [mid, theta3] + count_layertype += 1 + continue + elif count_material <= num_materials: + # Read materials + if line_material == 0: + mid = int(line[0].split(',')[0]) + mtype = int(line[1]) + sg.materials[mid] = sgm.MaterialSection() + sg.materials[mid].eff_props[3]['type'] = mtype + line_material += 1 + continue + + if mtype == 0: + # isotropic + if line_material == 1: + e = float(line[0]) + nu = float(line[1]) + sg.materials[mid].eff_props[3]['constants']['e'] = e + sg.materials[mid].eff_props[3]['constants']['nu'] = nu + line_material = 99 + continue + elif mtype == 1: + # orthotropic + if line_material == 1: + # pe = [] + e = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['e1'] = e[0] + sg.materials[mid].eff_props[3]['constants']['e2'] = e[1] + sg.materials[mid].eff_props[3]['constants']['e3'] = e[2] + # for ei in e: + # pe.append(ei) + # sg.materials[mid].property_elastic += e + line_material += 1 + continue + elif line_material == 2: + g = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['g12'] = g[0] + sg.materials[mid].eff_props[3]['constants']['g13'] = g[1] + sg.materials[mid].eff_props[3]['constants']['g23'] = g[2] + # for gi in g: + # pe.append(gi) + # self.materials[mid].property_elastic += g + line_material += 1 + continue + elif line_material == 3: + nu = list(map(float, line)) + sg.materials[mid].eff_props[3]['constants']['nu12'] = nu[0] + sg.materials[mid].eff_props[3]['constants']['nu13'] = nu[1] + sg.materials[mid].eff_props[3]['constants']['nu23'] = nu[2] + # for nui in nu: + # pe.append(nui) + # sg.materials[mid].property_elastic.append(pe) + line_material = 99 + continue + elif mtype == 2: + # anisotropic + continue + + if line_material == 99: + dens = float(line[0]) + sg.materials[mid].eff_props[3]['density'] = dens + # self.mate_propt[mid].append(dens) + line_material = 0 + count_material += 1 # next material + continue + continue + except: + # e = sys.exc_info()[0] + e = traceback.format_exc() + print(e) + + return sg + + +def readVABSOutHomo(fn, scrnout=True): + """Read VABS homogenization results + + :param fn: VABS output file name (e.g. example.sg.k) + :type fn: string + """ + sm = sgm.MaterialSection(1) + + linesRead = [] + keywordsIndex = {} + + with open(fn, 'r') as fin: + for lineNumber, line in enumerate(fin): + linesRead.append(line) + if 'The Mass Center of the Cross Section' in line: + keywordsIndex['mc'] = lineNumber + if 'The Geometric Center of the Cross Section' in line: + keywordsIndex['gc'] = lineNumber + if 'The Neutral Axes (or Tension Center) of the Cross Section' in line: + keywordsIndex['tc'] = lineNumber + if 'The Generalized Shear Center of the Cross Section in the User Coordinate System' in line: + keywordsIndex['sc'] = lineNumber + if 'Mass Per Unit Span' in line: + sm.eff_props[1]['density'] = float(line.split('=')[-1].strip()) + if 'Classical Stiffness Matrix' in line: + keywordsIndex['csm'] = lineNumber + if 'Classical Flexibility Matrix' in line: + keywordsIndex['cfm'] = lineNumber + if 'Timoshenko Stiffness Matrix' in line: + keywordsIndex['tsm'] = lineNumber + if 'Timoshenko Flexibility Matrix' in line: + keywordsIndex['tfm'] = lineNumber + if 'Mass Matrix' in line: + if 'Mass Center' in line: + keywordsIndex['massatcenter'] = lineNumber + else: + keywordsIndex['mass'] = lineNumber + + ln = keywordsIndex['mc'] + mc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + ln = keywordsIndex['gc'] + gc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + ln = keywordsIndex['tc'] + tc = np.array([ + 0.0, + float(linesRead[ln + 3].split()[-1]), + float(linesRead[ln + 4].split()[-1]) + ]) + + if 'sc' in keywordsIndex.keys(): + ln = keywordsIndex['sc'] + sm.eff_props[1]['shearcenter'] = [ + float(linesRead[ln + 3].split()[-1]), float(linesRead[ln + 4].split()[-1]) + ] + + ln = keywordsIndex['mass'] + sm.eff_props[1]['mass']['origin'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + try: + ln = keywordsIndex['massatcenter'] + sm.eff_props[1]['mass']['masscenter'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + sm.eff_props[1]['mass']['masscenter'] = sm.eff_props[1]['mass']['origin'] + + ln = keywordsIndex['csm'] + sm.eff_props[1]['stiffness']['classical'] = utl.textToMatrix(linesRead[ln + 3:ln + 7]) + ln = keywordsIndex['cfm'] + sm.eff_props[1]['compliance']['classical'] = utl.textToMatrix(linesRead[ln + 3:ln + 7]) + + try: + ln = keywordsIndex['tsm'] + sm.eff_props[1]['stiffness']['refined'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + if scrnout: + print('No timoshenko stiffness matrix found.') + else: + pass + try: + ln = keywordsIndex['tfm'] + sm.eff_props[1]['compliance']['refined'] = utl.textToMatrix(linesRead[ln + 3:ln + 9]) + except KeyError: + if scrnout: + print('No timoshenko flexibility matrix found.') + else: + pass + + return sm + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeVABSNodes(sg, fobj, sfi, sff): + for nid, ncoord in sg.nodes.items(): + fobj.write(sfi.format(nid)) + utl.writeFormatFloats(fobj, ncoord) + return + + +def writeVABSElements(sg, fobj, sfi): + for eid in sg.elementids2d: + elem = np.zeros(10, dtype=int) + elem[0] = eid + elem_type = 'quad' + if (len(sg.elements[eid]) == 3) or (len(sg.elements[eid]) == 6): + elem_type = 'tri' + for i, nid in enumerate(sg.elements[eid]): + if (i >= 3) and (elem_type == 'tri'): + elem[i+2] = nid + else: + elem[i+1] = nid + utl.writeFormatIntegers(fobj, elem) + + return + + +def writeVABSElementOrientations(sg, fobj, sfi, sff, fmt): + for eid, orient in sg.elem_orient.items(): + fobj.write(sfi.format(eid)) + fobj.write(sfi.format(sg.elem_prop[eid])) + theta1 = np.arctan2(orient[1][2], orient[1][1]) + fobj.write(sff.format(theta1)) + fobj.write('\n') + # utl.writeFormatFloats(fobj, np.hstack(orient)) + return + + +def writeVABSMOCombos(sg, fobj, sfi, sff): + for cid, combo in sg.mocombos.items(): + fobj.write((sfi + sfi + sff + '\n').format(cid, combo[0], combo[1])) + return + + +def writeVABSMaterials(sg, fobj, sfi, sff): + for mid, m in sg.materials.items(): + # itype = 0 + # if (m.type == 'orthotropic') or (m.type == 1): + # itype = 1 + # elif (m.type == 'anisotropic') or (m.type == 2): + # itype = 2 + mp = m.eff_props[3] + utl.writeFormatIntegers(fobj, (mid, mp['type'])) + + if mp['type'] == 0: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e'], mpc['nu']]) + elif mp['type'] == 1: + mpc = mp['constants'] + utl.writeFormatFloats(fobj, [mpc['e1'], mpc['e2'], mpc['e3']]) + utl.writeFormatFloats(fobj, [mpc['g12'], mpc['g13'], mpc['g23']]) + utl.writeFormatFloats(fobj, [mpc['nu12'], mpc['nu13'], mpc['nu23']]) + elif mp['type'] == 2: + for i in range(6): + for j in range(i, 6): + fobj.write(sff.format(mp['stiffness'][i][j])) + fobj.write('\n') + + utl.writeFormatFloats(fobj, (0, mp['density'])) + fobj.write('\n') + return + + +def writeVABSMacroData(sg, fobj, sff): + utl.writeFormatFloats(fobj, sg.macro_displacements) + fobj.write('\n') + utl.writeFormatFloats(fobj, sg.macro_rotations[0]) + utl.writeFormatFloats(fobj, sg.macro_rotations[1]) + utl.writeFormatFloats(fobj, sg.macro_rotations[2]) + fobj.write('\n') + if sg.model == 0: + utl.writeFormatFloats(fobj, sg.macro_loads) + else: + utl.writeFormatFloats(fobj, [sg.macro_loads[i] for i in [0, 3, 4, 5]]) + utl.writeFormatFloats(fobj, [sg.macro_loads[i] for i in [1, 2]]) + fobj.write('\n') + utl.writeFormatFloats(fobj, sg.macro_loads_dist[0]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[1]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[2]) + utl.writeFormatFloats(fobj, sg.macro_loads_dist[3]) + return + + +def writeVABSIn(sg, fn_vabs_in='', fmt=1): + if not isinstance(sg, sgm.StructureGene): + return + + # string format + sfi = '{:8d}' + sff = '{:16.6E}' + + with open(fn_vabs_in, 'w') as fout: + utl.writeFormatIntegers(fout, [fmt, len(sg.mocombos)]) + model = 0 + trapeze = 0 + vlasov = 0 + if sg.model == 1: + model = 1 + elif sg.model == 2: + model = 1 + vlasov = 1 + elif sg.model == 3: + trapeze = 1 + physics = 0 + if sg.physics == 1: + physics = 3 + utl.writeFormatIntegers( + fout, [model, sg.analysis, physics]) + line = [0, 0, trapeze, vlasov] + if (sg.initial_twist != 0.0) or (sg.initial_curvature[0] != 0.0) or (sg.initial_curvature[1] != 0.0): + line[0] = 1 + if (sg.oblique[0] != 1.0) or (sg.oblique[1] != 0.0): + line[1] = 1 + utl.writeFormatIntegers(fout, line) + if line[0] == 1: + utl.writeFormatFloats( + fout, + [sg.initial_twist, sg.initial_curvature[0], sg.initial_curvature[1]] + ) + if line[1] == 1: + utl.writeFormatFloats(fout, sg.oblique) + utl.writeFormatIntegers( + fout, [len(sg.nodes), len(sg.elements), len(sg.materials)]) + fout.write('\n') + + # Nodes + writeVABSNodes(sg, fout, sfi, sff) + fout.write('\n') + + # Elements + writeVABSElements(sg, fout, sfi) + fout.write('\n') + + # Element, layer type, theta 1 + writeVABSElementOrientations(sg, fout, sfi, sff, fmt) + fout.write('\n') + + if fmt == 1: + # Layer types + writeVABSMOCombos(sg, fout, sfi, sff) + fout.write('\n') + + # Materials + writeVABSMaterials(sg, fout, sfi, sff) + fout.write('\n') + + # Recover + if sg.analysis == 1: + writeVABSMacroData(sg, fout, sff) + fout.write('\n') + + return fn_vabs_in diff --git a/io/iovabs.pyc b/io/iovabs.pyc new file mode 100644 index 0000000000000000000000000000000000000000..e7a36a977537cd38735a66e00599764f5300d93f GIT binary patch literal 11394 zcmcgyO^h4Kb*^r5I74cFI6s;Z^-KM%mTRoW((GD_WjW4fq}A@)8f|*yj8++W-0E&p z%^{nd?xr-OC5=tAP8=YBkzfrU5(K&A=9>>jatM$(x#W@32-*Qy)m%;CE;2V8{!o|O#8cN;k*9~>g zOxrW6o>lcRbuX8;jH~*Dx|dh3sXoT2f^subIi%dIR3?=>CY34W=A<&M+;OSQD0f0C zMdjwDGOJu*z}QLkv8nFOk(l+$6bSQwiN7_+K=l=T3;3?$8~xjaQn6BXL)A^yRqCms zHciztRL@M?GpcK_y_Zq3P~XiEFDQPRQN1iF?u|(au23H*s5c>{FxCT2XaQM(f}gDF zO{qA?>A>4lQ+@x=L~mN{zBdt%tL`-Cx_dwoQ|6|5lbI=#mm^?rf zc&NdaA>!FY(Nr#4uXV*JufR4nJKenME(o&U&1jQIdiEA&{G6(TmPJ~LQG;Y4PE(z_ ze*dKEEvfDjmmj&kD1(>S?OUeWxtQr4QC<9l0x+ub#SJz zA}@<)IH7tcRri?co@Bog)p7dj23AYq^L_X%6T{ zWXDsLLGUGYjb;b{#8y)Eqnv(wT!nv~r4pnu?UX*028v1s0x2X9zetZ*P&N20ROB(K zr|ozoQm}vl@d;H2f0$F&if-lV;#Z)j4u;aaM%CI^<4?!B1=XDsQ~T-Iu&JH;O#O@s z=oF?BmwE-BHL>JFs(V<_V;RE5!qGt7tGU5c@P^nJz!kMdJ;M&f3z9cJm&Lo6DQt@4 zH1%?=TOK@=>Swe&nM6D0pvRcx{M1%Co%8Q<8pioQnn+em{H9{JEqY(eUSzg2gr_biHz%an5&8lr7J8%dQ>UNmVXp^U$)Z9BA^Y)Hgwb&j8HFSNRzfF(1Q-WdTuG(r^)Oj(?wiEvY2!uOs|I6?4RSRGd~2#dZ)! zg2ju8Kj`=l`>^g^(myNu=Z%;%)-J6izsd8;8-g|H06_Gpq%8oAqs8@t$dC8dDrfN^ zxLT2PZQJp~81bWu+6_PUc4K*XU`HMWCDTN;X3X$&Jq|lg+zGv#Uf{{#R@)1N4KL_4 zz0i)mc#NnZu(qxdsdR1lQ7qV+c3Xb^AeJu8mg_ZSQai94dlCAP;PyjQok+S=8!cOR z4$wJ5SOvi}ye201-3URbQg$1pxPm5uMDUFJyRB$wv*om!55PPL1;foqsA(jOw8K_= zix`CHZEr8q{GnxwFf2QbbgLEm7$c$#G)E4wQauq^o5Alb$RqUa#=7UkeoG6gnY1Nn zt5v}oRT$3;Ujus;G2oH9NbAr~Ae4!6XdqEd4{O-9EvFMc@`UXwA6pL^{!T}$#A$Uv zECnaEEecYUsvi@X3E8@(VefgNnCWzXP)fWEZ1{n<)o-EIJ;cfs0aPIX(gxkLojpZv zLRGaDVpL45Fq*CYJX#`b(n_ZG_|#V?O^^-goLbnP@5<8AQ$e5F@EUtz+!o7#gNW@* zLRS-gX)Mz8h78-5@{y(!)}&36d|T{k(G;aiO}cOyRnVlQh;E{`dUf-gkrzgroxqRn z&B$&vHlvFx7oxc1?rlaZ7eX)c?9h3*8Ai^-Er%Ao*^Fv!f75T_pWbP?y(hW}AwfE` zcLUl_^g4aY_z`2yIBe#O9KZWMb4JOSG>)LXfZs{;uu;G#C$%YqpBbZ=JQGqxFMcNh z;Wvk|8Tl5Iu`)VyDZx@i?=PZXQJ#zDq{)`~3`>4;#&M&Jast@qjYZ%n8Vkl$hBK7R zlg3FIcSJ^ET!MAdIDvlim>07FCuu4I({W(RgZ4>t23YaQ3N_&2peNlI<1zLK`A#Y=1@iP0YN-T1!I}#u!QKFLaukG^A03bgQvNw4{GS&_s00`% zK1Ym_rT}Rj(L3@W7|ws1oTH%U_`k6ai4Am*d;)qQ&Z8RVq^jdEqc2qox2%>gfF$@Fo8p8I#P)H0vcm_6l7O^IV&HW0) z&go%$UlE!Yo;xyOZ)Cz()U}OZvCqP>ELd>QtDEnTg9PW|uYi-(jNohu&Q}whn18*( z{UKsBDB-FEs1IAsR?Q3iCpM#jGPZ?IBaTLb)vG#~ZFS-_n$-bo#apf|m%UxP*=~5t zQEhqqH67(9(dyOm{wS4x(85 zy%wUs(04v_P|b%5Uv@QMu9z%Xb)cf>%B9YtNrKcB0$;nFB zZ4hF)ia69=j$VDHXLX+(hx%3+3R#ORR#+TiQDVXAtTQYqX^ZN!-emDb6cWhUVQBAJ zuk$g3q7yenL`-w;v=X&4R`|NNcOwj2p`5ad)l_u^8#RS{Mz9fAvP8bf3LeAKqIUz= z+qG7Soq;twt#z6YywJMAHi>`(>jTs++$1RLcUU+CQT$MgBT=LjpKB;diKU_B?N0m- z9M8X?Df%7?m79@RlF_4%1Jig=|6b_Fhl6(1%q2vbImC;{jB}ZTYI8=xd=HS5=u<7KZ~kG&n8CkI7Hefd4NtuEKH|jdH{jE zE&0jYOgSSxobF9jR&70~iy=i-Kl!c7Z;&JwfkPMclC(OVQpnw8jcLRc@EA>?Ss{F0 zMDXooP-g?#7NPzEmEC_$gnIq{&m`&MJ7Z*+nt*8a zF%BYIoCZdL)%nr>)rqAj>6S6HW?@cH@0f)3vUJCo($E^1yzq1~sp64m zy)Xxj6mc`G6A90Lq%c(?>92SCku`lZC)+;Pf1$$d%vxsg8jB8#M=UO|kd1;}um|E4 zZskDjJpo)yrrM0PDrrRe_{X}9>Pylb;|#H|eVV&1hfWeeMz^U%It9J@Pr&IeLVbP5 zO1?yOcTm~AHKMv}2RK$_3rkq;OVG!XcCuWaRE*pgmT*=4tjP7jr%yUW`V}g|qD)3a zK5(Y{mPTK=PIC4lr*kh!qC;!p*6md}2)suzQlk{(n}&#wexxtHAECqUUmC;v(kF5< zp!X%1u_-o;V*(foja(cQ%^AY0%JG7cJi%<~F25sokJ~yXtr(-XFyv@h9!%Ac&f@-) zE-v_`bOoZsQ5FQ%0Yj0na%f)vAhl$I3=Bi~4{7M6{y~m?P^W1PozXu?@sdtK)Td=u z*ytI%-=?!iJ_pOVe*ZSp@!cD2KK+m#M;;!?RB)!f)19BO;+nX7i?%4HY0;4lj`YM*EjliK z7+qUDcC>hgWN~R(45h3h3cXp-PAm_FwtC4M^ycvfA>v~s1mzxK{T_?QC}d9x?V#oj zJWxgE9N8(j9IW>_g6>NvFgVf7@S*Kk+$;oAe@QEZnNr{9AV8M$>*(8DBHaJ5L~>Kc zQ5@+NjMLaU3dSq=R8YeKpm81>$5CTeMMTC%NPS-V&kfo zh2hRKTp}85Q<#`92K-nrsdRdd!h~5~eIMBW7t()-t|8j7U=QLH%?ELc=7iHxB-~cP zk;~ARcD*+0m$d=Yk|g`{xghj@C$=hP%AlrcJ7H_f^>H6$J95VoiP6g4fCh%GSZ_hr zDGZQ2QUee)<35JVW#2P##d=p5Si2ZH&H4YwkipV*Tp8l{S4Cyas$j5F#*%z^PQGYf zHZJyCmJEx_Pn?8=)btEouCkpLAwx(XW`U)sEje1`t6=M}WY`PziFg`o7D|T_p_T9d z+|&^@-8}vgcEg}xLK4y>ex6}BEDxn9na7&w0VgVfci3YE9fHegyN*Ka-d}Pbk9Egg8wZ-gLM;{ly5xfgK}cPdr2IU=leHNO|YP8FeFZgjLF&N zfW(4cQ-Z9N@bb5k3F3_{%p%mW)I=pUAWpY2Xb*!qlNQ8(3=f_X{QISn=xc`CQc}mO zXZ`vt>j`oYP{%uappn@rPS!CQXo5K)Wy!)t4Z-nlQ3^+NH3G!xdkSP?1gAMR(gke^ z3qLVqymqs%`!jQq9~|4is0)BRkk>-VW%!a6V?qIj#$Gm`M1emt)w$$czwp^|zm!+~@TW5)a8=fWG zx?~MKK7N5frN%UVta?-DT_va(*NBI_^EWguFqeXgf& nB<~vDX6Q6L=GQvG=4I%dF>@&Y_jxxzo}bMh%IESY@(ce9JQb{Y literal 0 HcmV?d00001 diff --git a/io/utils.py b/io/utils.py new file mode 100644 index 0000000..aa948ad --- /dev/null +++ b/io/utils.py @@ -0,0 +1,88 @@ +import numpy as np + + +def writeFormatIntegers(file, numbers, fmt='8d', newline=True): + ''' Write a list of integers into a file + + :param file: The file object for writing + :type file: file + + :param numbers: The list of numbers that is going to be written + :type numbers: list + + :param fmt: The desired format for each number + :type fmt: string + + :param newline: If append the character ``\\n`` after writting all numbers or not + :type newline: bool + ''' + fmt = '{0:' + fmt + '}' + for i in numbers: + file.write(fmt.format(i)) + if newline: + file.write('\n') + return + + +def writeFormatFloats(file, numbers, fmt='16.6E', newline=True): + ''' Write a list of floats into a file + + :param file: The file object for writing + :type file: file + + :param numbers: The list of numbers that is going to be written + :type numbers: list + + :param fmt: The desired format for each number + :type fmt: string + + :param newline: If append the character ``\\n`` after writting all numbers or not + :type newline: bool + ''' + fmt = '{0:' + fmt + '}' + for f in numbers: + file.write(fmt.format(f)) + if newline: + file.write('\n') + return + + +def writeFormatFloatsMatrix(fobj, matrix, fmt='16.6E'): + for row in matrix: + writeFormatFloats(fobj, row, fmt) + return + + +def writeFormatIntegersMatrix(fobj, matrix, fmt='8d'): + for row in matrix: + writeFormatIntegers(fobj, row, fmt) + return + + +def textToMatrix(textList): + ''' Convert the text of a block of numbers into a matrix + + :param textList: The block of text representing a matrix + :type textList: list + + :return: A matrix of numbers + :rtype: numpy.array + + :Example: + + >>> lines = ['1 2 3', '4 5 6', '7 8 9'] + >>> utilities.textToMatrix(lines) + array([[1., 2., 3.], + [4., 5., 6.], + [7., 8., 9.]]) + ''' + matrix = [] + for line in textList: + line = line.strip() + line = line.split() + row = list(map(float, line)) + # for j in line: + # row.append(float(j)) + matrix.append(row) + + return np.array(matrix) diff --git a/io/utils.pyc b/io/utils.pyc new file mode 100644 index 0000000000000000000000000000000000000000..5dffa6b1091b191985445f2d7757a2d8c0a51fca GIT binary patch literal 2946 zcmeHJZBH9V5T3J*ow_CjwW-w?Mf;%@u0*WhB!pP4m5Ky4YD=Um6%w3N&Ub73==tut zyDqi@QdRiJ{Q>RF>;Z@9N64q1Z`ZSL&(7>TJ9Fi~7wbvu!P@}xr;PtEu;@2fh4?xU z1ME2vmOz#uEW=p|_R0`exVga13N^zjoH?+!Kr_P{%T);Pz!@&D5zf!84-hl@2cKb?9||zYn{| zWl%rlhrmcB4i3EJ;6V68YLNA0zVPEX-$M)BGf7QhuM5mUn#SC(iC^y?bq#t0 zPBhi^?G};2ke3Knnkq`x4tAX(a*$_C4TbzH3<04nw$JL$piv74fGfR)(q`BDMdSH< z(+SVd?~&hs{ujT8aq64@**y^#Onfiwc8S+T za=XNEmnh4Oc8NAFvBV{=yTofI@9=iWyNvh+@@_FJKF7{$)D1{S(vFndq#0jfHG5z~ z@F#Fm27~bz-!etv(=rC%8tjIEVU~>E=Gy%b#TLQXyNh9yFEvIc70uXYzP~Cxx*da? zxb-!If8yOcIHc*+0q^U_&J7t%nT-v0DHynmRIzLC(g2p5k4T>)Xu-zwh-9tYz`>%G zm|t1Oe2pGMZ<5pTb%+U=mJb(*37mr$ zVHcfzSfg(2uzb+(C@gKJ$sbY~PAEpsIQ9Bs5U0U!l!eceVV;C7ORgj*TJ$64KAY<< z7BZHS8HMn~rX!KB<{Zvx`EqR`dn#$BD(Q;v^7#cdZ2nTTwo9sI^R};)KhG6-F$sB9U5rEB05`L`OVcSrsc!#8a_O!ZWcUp0D&7ljc%B!S4pC+oaY`YQ5EB zh_87T@73DhUu&<54wlF5{wfcUMfB`H!QQ7>)~}5`!`=-n&)faJMQhzbKXYoOl%H8C zLlyX9(t+?rocMN#^95g3`w76GX^9=Pit(jp(A*`=^&Bmc?(OlNE>`i-df3*q9V}Pz zDmUg`kf;;5VJNf! literal 0 HcmV?d00001 diff --git a/ms/__pycache__/analysis.cpython-37.pyc b/ms/__pycache__/analysis.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..86e02aefb4e41c54fae641bf26be71a40ff9b30a GIT binary patch literal 2710 zcmZ`5+in}Vb;!9$8u_9)S|{D23$_asD>N2j2MJmu=%#j(X3^BUsGTk#s0~I#N*ZS_ zT9R@kV+0gv_p$$AKlUs7AAIdo-iki%Q_tayoOFvFiR9toxgYY7N1L0W17H5xCx7m? z9p_&-S$_clPoSy4K*t?UQfE+&WI*hk4qQ8W1K*Cppkc>m)}jt^vT)EQ6~{)6*&J-u znC+}%X>JX60Mp_w_x?r(w|U5Y9stzkZQkHbfIi~n)ah-!11p@9vzX!tUkjQK?HHPx zLZ_W8QaaaUMHlx=dQPw0Yv-x+r=Jb!%3XQdyYfr#JmBOG()gv1G^P94`DXWPM>npT zB?U~tDXh~~pj%~7lE2bb1J=T_QBuG*%O;rOUhcEwNXjKE=9)!pk`|*VWg`(~EXt#F zsS?HJD#<78x6h8wLI&TZ-&d3Vn@g|HBi(BocP2-sF_Yp5xJ*-tSSLlUjGt(cDbvW3 z;vyO;vn@r$|M;7u(=X@xMUfT8%cD%#^@$j1(|#*~1J!(yn}(d{2%0u>Ap~5#)}V=| zyIj){~Yioj<5~??h@x~Cl3j5;if%x zN{0t`VC*$c1|A0KJHjq^A-imkLFjYLIHLRgerU+s9yNi|9D|3WkP^}lGf2O8UO;?Z zwN&zBFyHt#jm96Rb2SD0XlA?^0&2WsCUQdy*b-b^G3B-RD+OuhS0aJ0cEv}MO z<4VZgT|CEb8J3U}?te?6_r9fL`Xe}%UFVLo+9=7*3asPGYV!u?uc=pCfU~ozw{+Ve z+v5RP)PX#2c0ph1V9y)39p3!0Q)}NU!%ExxmD;MUbEpGdpsQAp8~I@7gM`+N1oYS~ zw^2?_-}zswAEWNA4|UfXw!cSTdsNY%gYF-&+3Lawm7Use8^1BNoZG;+#@|S^8arhR zD%&l9{qk1XE;q)Qv7K@U_IVhcK)Ge7cs!Y&3Z`_VlbFq9F-k?oiZPVGnM_oM=zA

5Bm>))w5~=e6-(VV_9U@yirlSvD7C;Ub)7Oi)>aXD7g#^H#?h(0`7s)l0EIS z7w|9mL{#=hQoto3c&+eo7D-@8MLyBfLv{wfrS%NMh{45@Y}Fm= z!H6`8@>qZ%V0^TA#NOi=Bts^P52dnL7E0~aHU0w~@3{^kc}gSGdi;OAi~UMpTEsSZ z?Cn&5muP!+IWV~R-a@v~Yb=YEBgj?>L)-+53w|~5?|-g=&rjaDfbY?FSleX}x_5Zn zuQ{I7z}+Xkz|eeVe7J<8r6Fg=<%Kr>r2scgZ~?nCe$%^cf-2f3eD9{Ai+w}K*bWAC z{J_u*;Oqd}14Et<8e$*zDi;I#d@-QEUl{TSLk4mW%rNA4hR&r4Dqp`>Vohi-v9De} zJ3f5~)avKZIo*$;!ZzWz<#x$0;YNcW4TC4-h~9S($z7Li?EV+T^TCn; literal 0 HcmV?d00001 diff --git a/ms/__pycache__/beam.cpython-37.pyc b/ms/__pycache__/beam.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..750cc904e331fe9a5a79522867895d235e1e7b45 GIT binary patch literal 2603 zcma)8UvJws5Eu2&vJ)rGmZpC;Y(mj(@vsDKFcd@34PDZr1GX9k(xHRuR!~IQQ7Xv} zDQB^pJhiX=4(;CdmH4`+euX{lca-BKZBUE`k9Wu8@kIXa=#6HxVW1VuUp)I|)-e7d z;kemg?x4#bK)AtKW&})+W^M&GGwvJQv6QMV`ejrmVRTbD<)5i6gSo69M_JevQJ#b;SKgGX#Y3@A zs%A0FL!R!ZGVK?##gr|Rtoyein%L>zdz-(@0XXIP|V)R9q^b>xk;6^F(%X*&C#w9Xh4e!AN;wm|}ex*5Oos0%9ah6w%w`?kzBwYqJJ~iH0Z>uxZM<*;SdPf)2vXxm3loJUN zjRZMWmK^3SCIH3i!}ZsTIhvZq3C6!AMX%i1+Pe&mBp(JFUp^i0X`^<_r7te%&_mS- z!@SRj8S%4WxHpWl$xJN_c|Q(AF%N%Ykq907DMA4V+Kz&hE{JE^rZ=y3F~n&jqdVyG z2FS;u&17L5x$b8&l6vkVi=c0ku)^WAMK>l zk9979$)mK$Mmj6Xoy5O?Z{v|qDEL@XaKxy=k30)7Kv~_K~|;ZNX=(y0qjo|f-0BtCRIVInL*qycuG!XFq5AM zWD{$(a&<<^puv%^BI`$72F)_flkTvnOoBN|9ZJ`T?f@xz9!xcct`Puy5M)=CRn=&_ zww7(3Z^*lGldk*2BqNwnVb)3Z7nlV!$S=otV4U3p*!d;Pc=^Zw*~YhZP4!<4QojB@KU()+9E`?KZaiJT zd297na#V>o&x&iuV#`z(76um9j29jx`8HDMZqlW@!~Kt+=9TJ2KJ^ZBN$#w7!9$CjGS)x(eP;u6%qtbqG4cI_=*XUXVHh_*p^N>9^#q6QAW$ev9#vcRJ zE!J+C`h&!`Sk|(uXg>q5SgKh3ev}Ki5B{*!_}Cb=Q49Wnt@J<)wx~lLFX;e3Fg!n6 z+VBz+@6pDRt!aDY?s%})VKZD^k z_6DLq75f>(cFPkNX~QKVmx$0qR=qnb#L=rE)`*<<0$stx=jf6!;n{eBzVa44SCGEt zDF;<@FcJ%pl&9bRRoPP}x+Qz5a{Jr62wT5dddYcvDDQUlO8!LuaH^+{%~}^O{RduR BWXu2n literal 0 HcmV?d00001 diff --git a/ms/__pycache__/iogebt.cpython-37.pyc b/ms/__pycache__/iogebt.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..9bcb3a3a51d3858af507c6bb51da57f754d0625e GIT binary patch literal 5907 zcmb7I&2JmW72nw}E|;I8WPMq76eo4fAU0yxZQRsOUB`Al+9*lw)&ZfIp*SmgC2~p4 zE*)CkMGs1WAV?7%il8WZkS;+B6hV#y?Hb5H*e-zspM+--RsJ`=g(`}zo;_)87N%E7Y0bU##vWuFqP{)L)VzrGaD95 zTkKx0;bi4}&sD93MiCU98}Dk35;wVpyv%K$Ltfzy&m*sLmlu#v@d_{T^1G~2(&@U0l@@ zGqGZ8&qPn;T;J5i^INC3PR05j(QR(*Gj1+u)MKA9ZR^w>?G9%2hW1nKO^ut09pxe? z%Ex+aqnEv}?dUr!$;CMu6|<&`lBG%ol+XvXV$3#NZbH_dW1CygYKa}3K`BQi)#jk> zpJT2&3}Z^;o{yC0#16?z^2i}AamT$!x1(}w5NEUtXKY5I>rR`)^Zht250z>}q3Yt-~3;^6ne4&0hyb1z)&g`I(qbk*;ZX6FVEq$|7@fAPk(+b;*?#o!v=--0@#`P%1MAG_*U8F~jEBQy4SDPJvR!KkYk%KH?(yO@M;dp4M6c z#WN@XAlT@1tid@>CVI>wEn)x+eb3leZclE2jm5THY%+o#wo3X{*wYq)brlZ`j4y z_rNoxEHR)Y+W?X!rYe~zS(uECFNe>GtPJ(az=b^ddV(oTuuD%0Tu+MlO1o?{hzsE2 z?yybxxcEGF6-MNgpO(3?sfQnMGcIr1$Sumn@8SZtVVJ^JW0wU_@|Q(H((A@dGo>WNU!Pwy6JWnGbwHi4jmjg|Z7UohEVo8;1Bq^>C0_^4Es>0$% z@Ch$vWB&p)AUP%2#EeS|+J>7<5ic4CYt%+>M}>Gw;r{JpI@U=$$Q-0H(X)I+UWLX5 zOOw{hk=9h3+E^2y;&p6}3a={dYNTCUiY?5)5!>ihN%e)4%wp(mV=GOq2wG|dQi68o zwsGakm2Knf*|UO>N}9I-6?I)&3K67D4(SPLy-olj2)G1oX7PRMY`>~n1??w6+eSQx zq+S$rR9-+Ljjk7n6I6YI67n~hr#-1~A>GEb14WIJ$0<3Eq&_1HV=TeO7T*F*me#wi zPIHWQKnsF3@idJ&OYJ%M|F8=fE2nVQ9M>!9EW$BCe4D!GfgjCOjdN5tLZF$^^SP)~ zwU~`o;3>XKBrEW`UQfCyyQcEEjfXvM>X66XiU7jx7RH9FNHk{S51`OZc%NwkZqgE9 zrwBlC&`PO`R)7%tQLT(VE&!;gXaE|0hPeP9hgH!=y>(X91)-lTRmBWKMh+r}v%(Kg=EYE_pY;JU$oRRNw~u*hotF%JAjuabcH*dvPu% zZy~?0U~eMc2;YP6R^Y?$ahr1S5xkS4j^fr`_%QrawSEHc$YJmP9X`Lmey96m%t$vo`59v880Ifa>qc<_SGM z<25x?U0jdnl%%=%*q##=qGD81@^W(;Y*S97NYa!n-qDiKzlj$W-^H(81#moA1>!bY2k`rvbvXVF*5UYM9rUfwVbukS>p(kxunw@&XYnj( z3q(^h%hrL)7;`M1g$B(m>QHt+j5yz1c89@TD>EdBdh+v)3q%#jD`-VL+`O0#s$OkHI|r6B5l z0#xTM9RNxS2nXwbm8cXy9|{T!hb(##L0+h?+?WC#CJqyCP5IdpT)C~2QO-7}B$$&dRej&y$f0YS9OoxK7m#eJqco+x=L0r9kr@(xpZ zCCjT>jyO%xf%Ph%&hnWoKbGaQJFJtZNLbCz;I5j{M8@bM$=I9Qr{Ml{cH5`c0(3G{ zjj~x=G#B6v?)OS-rQ_&LjIU(lU%_}AH+wtIZIsoWo^ra;=QifE`g~Sj0BsT9@s9Jb zc95}lM8-B`RH5G>ag^jDbaUv8oajd@=Wv{dyh8KKG zV(g<@pv7MgYk|}367-l#=eo<@I)mBf@s2>&B8^~=Y}MS2$7q&(8s7}P3Sjq`yLihO zeH=g3hFR6{B9QK@g0ir__1imI5RbN_02vzw9yLq)UsW2Rq0=v*e%=lqs6%NNPZt5;b&~Ur{0$q()T5G-CW#0v3y71 z1!0+P;LA~eIlAX9(|vi_4+aRv2^qG{b7+<9orZpir{^Nw9~2;M=YH|bbC-DHO6E57 zi+tNg{UT2k9l4v{9^7`nH|Pg=Y6u&~`R{RQtD`U!1u&{JX>s+Qp}KYbR?zk+!pB2` z$7R0VZG|C1>#jOvrP)JO+I``7{9vX2xHJPiTcp+3IZtuP(hq%`V65j(~I?%2MR_y!By_^)LH)_n>E&Ry4$kDJ{&X z7mCM7;xZ-b`g95_{D-5qA7WIVJi1bHm8*hEn(Y(o#k777eCXAN87|u43vWPES(O;m{GR4rBE=<6Q3b f`LOFf4+Ci-y}dZl8nH~lx$cym8_qAainH)P3twX3 literal 0 HcmV?d00001 diff --git a/ms/__pycache__/prebeam.cpython-37.pyc b/ms/__pycache__/prebeam.cpython-37.pyc new file mode 100644 index 0000000000000000000000000000000000000000..fd2132797a53ceb69bfab8dda548430f2ef01526 GIT binary patch literal 7564 zcmcgw$!{E2datGTjTG6WxQS}Zv?Nj{MN+mFirOejvSlr_sgdkAl-tFA)y*b*rM9no5gF-3fQh=%%~&U{?g^so>s%k!q*+YL zXeY%D?2Wl;83SzY>zkQ3I%wGt6FI=-LQM1k(<9BzERSs&VlSk{4`{tgb^@~FJn1E5 zymNx5uu94#Pj4l~ufUOfqbo{^XFy50eQv*-m#NK6XJB9J4DM@h4D3j|Cp$wjy{~tM zc^3Q9GQpEu#$HtXs;YSzu6x<`3gnOQ+`hi7ybW=ei9$=@?yKJ@akPig2w)W;Cx0hvei zPC(P390b%yaUd_elX75lv~vntP7WjXI^y0L@2BzeMaqOkg|eV4(N`EIbPfN0o)-}Chmby}W=m**7zWEhKxc->ul1zA8YkRQq>AKo(e zK&6TdE2;o_$j|Y*LWJH^RZBI6{8PVC_PsO5u*nenzI07l55P4gC*%kFhIdw;4fmW?d$0oSEI17B9N-JscTS!yqm#Y!SbwN`X0pBHPV?9s z>W|CbRD0%Nb;i4ZT{AK*r!W)F%L}m6kdtysUML&rhRF8i7j!?5nD(oh@0b0E*hTCy zH)m!4qQ*zP4^_Oz4&wJWs3nRtFx^_ho)5t{CkMbcFZ)!dZ9po`5FK%iikskygX`20 zu4#oCNBmFl)9`5qy)!ecc^Bn0I7j8hVGQD7Ew6dAe3XyM>5A!1b4)m{1DeCGxr%{# zc8Tk9R^f0Dsq=DH&MWFiZC7=8K~0ddcS&A?9jc4N*^yS5_qvuHOs9{+I-vb${)GB5 zZ2VYte-EPjk^Gpy|5wP2uI)v)A0n>F#sa3<1bQ?l{+*}2%eb!Kx~jbX0RCQ+*HnMT z#GN!vs<>G<7#gKY6jQjbn{pgDu z%7RI7+>|$=cM8xgc`M|&UD3QdxR%_@{K6aZM{oZv_Koll)r{0%{-t*pGjv(r1@ClM zu3YBRs25${K-8D{3}z6~XCV8Yyr=Anivqvs-Iw=}?^)m<$OoalzmpGuzt1lsTG$oF z{a#>EHy>+;_lftZ{6yZ!Yw}aNcQv6lmib($aSoCn%7;pFLL~US_nCYMjhBFbBp-#A zd?_CR{~4cG8ebtw4+CwOj7nDHJ5WCY6*9n|;va+mv4SoGdZM5!fSzLVJnr&(E1~UA zj@#~ij$Ka=c75Sq>0Fg7@^kd#GQYa0p--<-)Pjl+Ad|bmW^d3CbHo*V@bO4Z>npTtz&myj8qY+>~*g zuLjQ6TxXH@@ayg^IrT<|UT^?j#d)-NQj@u2W-qb*44N+U8&nynl4rPnh6?$YuG~HS z^;1=Q5)VH8yCe&Nv1%C15kOsb__zRMMMCUVp8?n8!sQrLW#9Q$`KcX|a{aon` z_?H&7hW}2<6x+ZXitYEkiNLB~9{5MPbcH$eFS-58F1hm5mbS9Dq_o?Z2cLmjmSsSn0;Cf8|^9ZjPS?m+!^z#vWh8$H&%HYRhAaNmQ`D%B%`18a?@rXO8^9Wb4YhB|&5(Xo(eo5t^dRXJJ;bHsmoj;JOOr60>vX8XH$yAwV4kw# zpMNkJ&TLj@O8{(_<)>pwhm0FF``p61v zTH#GFElvECerDakujs;xuv?;8vK^->3e9ri!M&wtsZ;@V%i zvkMKYZqqh*!zNbRU8q~tMxokh?YM<9?KwZwuGc058SZY&J`~}&qS~k&6JKv~d!g_c zcM>ah8V%Y(`j36)3e$y>ur1dvtUf55FNvn(OgnbTtu`BzN0d8%$~q5L$OYg3LcDY+ zw^D1aTeZU6a)JCkBJ+qRfwk3!t8jw;hVTa8%RmYEzw|=&6wTo5P78V)jCmTEU zqT|{v=PmtV!E&^=s};NP%BtbASlZN04(X+^;QjbIxM+vIW`Ff zPtDrhY0=YeRNM`9nkdfbll# zX;7SN@SPX!1}8~&8LC(m zMa>|eEUjD*7H|I>A4VC*8W~J=Y0Xw3@mqm-L!eVRQobbdb%c8Ib+TC7Z8-$`dXbq} z`1-60(@u*SC5OeUjZ$rg+i(AiG$Dgw+7J*Yu)){oSX?cjG*NUa4JG*uIn-5u%pX`1 zjD@g{{%~vZu5Z?zbrxAdpZSLEvM5H4gEFkxu4p!0-)N$lS{5>$a9TC=QX7R$1;RECr z^NmIebs&Vb%i;n(Vzpozsk3S|7E{gVurwLJ8unhyYPC>C>;&yvaql)!9!tAVs6S&R zS~vk|H5ga{$v08Wa4^8ZHm6c$J&#fQ0wor16*+RjEIAI#k{^d97DH=TwF-++DYK|* zb{2buswr);=xZAj(vPZY_f6iULRO=s%2`Z7j_4=R(68o)Bk~ln0m7wQ2dxVkEzE3$ zYsamz+{!@<%rCF5pcJbOG*!uFX>effAVIF}$J}NRVh7Pq+OV63q0LgOJL~sqwi;I^ zLheU4G0FTW)eX3WU2nd!CmRTMaOU_Cnz!N_Z1N*iR4n2O8_8%3btbWN8}lrjzATAF zFiTksr%}1?0%a9LFa>HjOd&s*sw~oK%u`Bh*1BC|32af>MLnw78bchAJQiJdXkUEY z2{y*d>j4l0hf}FDW8G2QLCYGg1~XbE76lbO6)QQxq%H;QRQyd}N3#((^|67{iD0oV z`Y8NHYZliWq|PyOeQqAWC0z4eyVp?<)lFZo`ub*&MB`gXd*uXHCjbQI8Lu2>e(RQ* z`N}PeLCWvJxMe8>8GTG$#nKdT2M0iRfyodoGBE3KBT;QukT`yg;WC)Y2w3Q{sEtwx)==@ zyNIr(=C2YMjB66DCbq)ia6)|0&uWyTPch0J77ftGRtx#Y!Ck34Vua+!f|Zz{1qF{{ zN2wlUh#po5jUS>ucscz5t`YwyAe4uDtifj}T@Sv&< zM>6X@gahG>%C7_qj5g*BU|1BS60Jg$SydoPQ_-5Gf6&BJ82_K2++A4vnS(QOZKqLn ztu+U8ZOxgTpLX0GzPskkPh;xgI4^Ci38%DykfE*C>Q1Fq#eJTPmBOyYsI-Kn>B#btsA@P0FcHwoLG2~V46bZSHu zilvX>#Ch7ERi!;scV^a6%f5k*e*nBVN5n)xbSlhnDq?|9{Q>IO>$S)dCE`@jwCRAp r<5p{EWc3g?Ngx@lTJXqi1vo#kh=XpYX5XYA6b>Evlo(2=3WWJTjq_~} literal 0 HcmV?d00001 diff --git a/ms/analysis.py b/ms/analysis.py new file mode 100644 index 0000000..e8887e0 --- /dev/null +++ b/ms/analysis.py @@ -0,0 +1,105 @@ +import os +import sys +import datetime as dt +import subprocess as sbp +import numpy as np +import msgpi.ms.prebeam as prb +import msgpi.ms.iogebt as miogebt +import msgpi.io.iovabs as miovabs + + +def solveGEBT(beam_xml, scrnout=True): + """ Carry out a global beam analysis using GEBT + """ + + # Preprocess + beam = prb.preBeam(beam_xml) + + for sid, mso in beam.sections.items(): + section = miovabs.readVABSOutHomo(mso.name+'.sg.k') + beam.sections[sid] = section + + fn_gebt_in = miogebt.writeGEBTIn(beam, beam.name+'.dat') + + # Solve + fn_gebt_out = runGEBT(fn_gebt_in, scrnout) + + # Parse results + results = miogebt.readGEBTOut(fn_gebt_out, beam) + + return results + + +def runGEBT(fn_input, scrnout=True): + cmd = ['gebt', fn_input] + cmd = ' '.join(cmd) + + try: + sys.stdout.write(' - running gebt...\n') + sys.stdout.flush() + if scrnout: + sbp.call(cmd) + else: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + except: + # sys.stdout.write('failed\n') + # sys.stdout.flush() + return + + return fn_input + '.out' + + +def solvePLECS(length, compliance, x1, f1=0, f2=0, f3=0, m1=0, m2=0, m3=0): + """ Solve the static problem of a prismatic, linearly elastic, + cantilever beam. Find the three displacements and three rotations + of a point x1 given loads f1, f2, f3, m1, m2, m3 applied at the + tip. + + Equations (5.59) and (5.61) from the book + Nonlinear Composite Beam Theory by D. H. Hodges + are used. + + :param length: Total length of the beam + :type length: float + + :param compliance: The 6x6 compliance matrix of the beam cross-section + :type compliance: list of list of float + + :param x1: The location where the result is wanted. + :type x1: float + """ + + F = np.array([[f1, f2, f3, m1, m2, m3]]).T + e1tilde = np.array([ + [0.0, 0.0, 0.0], + [0.0, 0.0, -1.0], + [0.0, 1.0, 0.0] + ]) + + Fx = np.array([[f1, f2, f3]]).T + Mx = np.array([[m1, m2, m3]]).T + (length - x1) * np.dot(e1tilde, F[:3, :]) + + R = compliance[:3, :3] + Z = compliance[:3, 3:] + T = compliance[3:, 3:] + + K = np.zeros((6, 6)) + + # Equation (5.61) + K[:3, :3] = ( + x1 * R + + (x1 * length - x1**2 / 2.0) * np.dot(Z, e1tilde) - + x1**2 / 2.0 * np.dot(e1tilde, Z.T) - + (x1**2 * length / 2.0 - x1**3 / 6.0) * np.dot(e1tilde, np.dot(T, e1tilde)) + ) + K[:3, 3:] = x1 * Z - x1**2 / 2.0 * np.dot(e1tilde, T) + + # Equation (5.59) + K[3:, :3] = x1 * Z.T + (x1 * length - x1**2 / 2.0) * np.dot(T, e1tilde) + K[3:, 3:] = x1 * T + + ur = np.dot(K, F) + result = np.vstack((ur, Fx, Mx)) + + return result diff --git a/ms/beam.py b/ms/beam.py new file mode 100644 index 0000000..80ed51f --- /dev/null +++ b/ms/beam.py @@ -0,0 +1,98 @@ +import numpy as np +import msgpi.sg as sgi + + +class BeamSegment(object): + """ Class for a beam segment. + """ + + def __init__(self): + self.points = [] #: Point labels [start, stop] + self.coords = [] #: Coordinates [[x1, x2, x3], [x1, x2, x3]] + self.css = [] #: Cross-section labels [start, stop] + + self.rotate_a1 = 0.0 + self.twist = 0.0 #: Twist + # self.dihedral = 0. #: Dihedral + # self.sweep = 0. #: Sweep + self.local_frame_id = 0 #: Local frame + self.frame_id = 0 + self.curv_id = 0 + + self.num_divisions = 1 #: Number of divisions for meshing + + def calcLengthSq(self): + return ((np.array(self.coords[0]) - np.array(self.coords[1]))**2).sum() + + +class Beam(object): + """ A slender beam-like structure + + This class is mainly for the GEBT code. + """ + + def __init__(self): + self.name = '' #: Name of the beam + + # Analysis + self.analysis_type = 0 #: Analysis type (GEBT) + self.max_iteration = 1 #: Max iteration + self.num_steps = 1 #: Number of analysis steps + self.num_eigens = 0 #: Number of eigen analysis resutls + + # Design + self.angular_velocity = [] #: Angular velocity of the rotating beam + self.linear_velocity = [] #: Linear velocity of the first key point + + #: Points + #: {ptid: [x1, x2, x3], ...} + self.points = {} + + #: Beam segments + #: {bsid: BeamSegment object, ...} + self.segments = {} + + # self.divisions = [] #: Divisions of members + self.pconditions = [] #: Point conditions (B.C. and loads) + self.mconditions = [] #: Member conditions (B.C. and loads) + + #: Effective properties of cross-sections + #: {sid: MaterialSection object, ...} + self.sections = {} + + self.frames = {} #: Local self.frames + self.distrloads = [] + self.timefunctions = [] + self.initcurvatures = [] + + # Results + # self.results = None #: Results of GEBT analysis + + def echo(self): + print('') + print('Key Point Coordinates') + print('='*40) + for pid, coords in self.points.items(): + print('point {pid:4d}: {xyz[0]:16.6e} {xyz[1]:16.6e} {xyz[2]:16.6e}'.format(pid=pid, xyz=coords)) + print('') + print('Member Definition') + print('='*40) + for mid, bs in self.segments.items(): + print( + 'member {mid:4d}: {pids[0]:4d} {pids[1]:4d} {cs[0]:4d} {cs[1]:4d} {frm:4d} {ndiv:4d} {curv:4d}'.format( + mid=mid, pids=bs.points, cs=bs.css, frm=bs.frame_id, ndiv=bs.num_divisions, curv=bs.curv_id + ) + ) + print('') + + def findPtCoordByName(self, name): + for i, c in self.points.items(): + if i == name: + return c + return None + + def findSectionByName(self, name): + for i, s in self.sections.items(): + if s.name == name: + return i + return 0 diff --git a/ms/iogebt.py b/ms/iogebt.py new file mode 100644 index 0000000..0d7103d --- /dev/null +++ b/ms/iogebt.py @@ -0,0 +1,419 @@ +import numpy as np +import msgpi.ms.beam as msb +import msgpi.sg as sgi +import msgpi.io.utils as utl + + +# ==================================================================== +# +# +# READERS +# +# +# ==================================================================== + +def readGEBTIn(fn_gebt_in): + beam = msb.Beam() + + # results = {} + lines = [] + + with open(fn_gebt_in, 'r') as fin: + for ln, line in enumerate(fin): + line = line.strip() + if line == '': + continue + else: + lines.append(line.split('#')[0].strip()) + + li = 0 + + line1 = list(map(int, lines[li].split())) + beam.analysis_type = line1[0] + beam.max_iteration = line1[1] + beam.num_steps = line1[2] + if beam.analysis_type == 0: + li += 1 + else: + li += 5 + if beam.analysis_type == 3: + beam.num_eigens = int(lines[li]) + li += 1 + + line2 = list(map(int, lines[li].split())) + nkp = line2[0] + nmemb = line2[1] + # results['ncondpt'] = line2[2] + # results['nmate'] = line2[3] + # results['nframe'] = line2[4] + # results['ncondmb'] = line2[5] + # results['ndistr'] = line2[6] + # results['ntimefun'] = line2[7] + # results['ncurv'] = line2[8] + li += 1 + + # Read keypoints + for i in range(li, li + nkp): + p = lines[i].strip().split()[:4] + beam.points[int(p[0])] = list(map(float, p[1:])) + li += nkp + + # Read members + for i in range(li, li + nmemb): + data = list(map(int, lines[i].strip().split()[:8])) + bs = msb.BeamSegment() + bs.points = data[1:3] + bs.css = data[3:5] + bs.frame_id = data[5] + bs.num_divisions = data[6] + bs.curv_id = data[7] + beam.segments[data[0]] = bs + li += nmemb + + return beam + + + + + + + +# -------------------------------------------------------------------- + +def readGEBTOutNode(lines): + out = [] + for l in lines: + out += list(map(float, l.strip().split())) + return out + + + + + + + + + +def readGEBTOutStatic(fn_gebt_out, beam): + flag_analysis = beam.analysis_type + nstep = beam.num_steps + nkp = len(beam.points) + nmemb = len(beam.segments) + + results = [] # for all steps + points_results = [] + members_results = [] + + # lines = [] + + with open(fn_gebt_out, 'r') as fin: + all_lines = fin.readlines() + + i = 0 + sid = 0 # at least one step + mid = 0 + while i < len(all_lines): + line = all_lines[i].strip() + if (line == '') or ('===' in line) or ('---' in line): + i += 1 + continue + elif ('Step' in line) or ((nstep == 1) and (sid == 0)): + sid += 1 # step number + results_step = {'point': [], 'member': []} # point results, member results + mid = 0 + i += 1 + continue + elif 'Point' in line: + # Read point result + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + i += 5 + results_step['point'].append(out) + # points_results.append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + i += 2 + member_out = [] + for j in range(nelem): + if flag_analysis == 0: + # Static + out = readGEBTOutNode(all_lines[i:i + 3]) + i += 4 + else: + # Dynamic + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + member_out.append(out) + results_step['member'].append(member_out) + # members_results.append(member_out) + if mid == nmemb: + results.append(results_step) + continue + i += 1 + + return results + + + + + + + +# -------------------------------------------------------------------- + +def readGEBTOutEigen(fn_gebt_out, beam): + + # if len(gebtin) == 0: + # gebtin = readGEBTIn(gebtin_name) + # gebtout_name = gebtin_name + '.out' + + nstep = beam.num_steps + nkp = len(beam.points) + nmemb = len(beam.segments) + members_in = beam.segments + # members_in = gebtin['members'] + # print 'nmemb =', nmemb + + # evei = [p1, m11, m12, ..., m1n, p2] + # p1/m1i = [x1, x2, x3, t1, t2, t3] + + with open(fn_gebt_out, 'r') as fin: + all_lines = fin.readlines() + + + + + # Read steady state results + + results_steady = [] # for all steps + points_results = [] + members_results = [] + + end_of_steady = False + + i = 0 + sid = 0 # at least one step + mid = 0 + # while i < len(all_lines): + while not end_of_steady: + line = all_lines[i].strip() + if (line == '') or ('===' in line) or ('---' in line): + i += 1 + continue + elif ('Step' in line) or ((nstep == 1) and (sid == 0)): + sid += 1 # step number + results_step = {'point': [], 'member': []} # point results, member results + mid = 0 + i += 1 + continue + elif 'Point' in line: + # Read point result + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + i += 5 + results_step['point'].append(out) + # points_results.append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + i += 2 + member_out = [] + for j in range(nelem): + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + member_out.append(out) + results_step['member'].append(member_out) + # members_results.append(member_out) + if mid == nmemb: + results_steady.append(results_step) + end_of_steady = True + continue + i += 1 + + + + + # Read eigen results + all_lines = all_lines[i:] + + eva = [] # Eigenvalues: [[eva11, eva12], [eva21, eva22], ...] + eve = [] # Eigenvectors: [eve1, eve2, ...] + + eig = 0 + i = 0 + block = '' + while i < len(all_lines): + # print i + line = all_lines[i].strip() + if line == '': + i += 1 + continue + if 'Eigenvalue' in line: + eig += 1 + # print 'eigenvalue:', eig + eva.append(list(map(float, all_lines[i + 1].split()))) + evei = [[], []] + # store all Point results for an Eigenvalue temporarily + # evei_p = [] # = [[x11, x12, x13, u11, u12, u13, t11, t12, t13], ...] + # store all Member results for an Eigenvalue temporarily + # evei_m = [] # = [[[x11, x12, x13, u11, u12, u13, t11, t12, t13], [], ...], ...] + mid = 0 # Member id + pid = 0 + i += 2 + continue + if eig > 0: + if 'Point' in line: + pid += 1 + out = readGEBTOutNode(all_lines[i + 2:i + 5]) + # pr = utl.textToMatrix(all_lines[i + 2:i + 5]) # point results + # x, u, t = pr[0], pr[1][:3], pr[1][3:] + # [x11, x12, x13, u11, u12, u13, t11, t12, t13] + # eve_p.append(pr[0] + pr[1]) + i += 5 + evei[0].append(out) + continue + elif 'Member' in line: + # Read member result + mid += 1 + nelem = beam.segments[mid].num_divisions # number of divisions + # nelem = members_in[mid - 1][-2] # number of divisions + i += 2 + evei_m = [] + for j in range(nelem): + out = readGEBTOutNode(all_lines[i:i + 4]) + i += 5 + evei_m.append(out) + evei[1].append(evei_m) + # members_results.append(member_out) + if mid == nmemb: + eve.append(evei) + continue + i += 1 + + results_eigen = [eva, eve] + + return results_steady, results_eigen + + + + + + + + + +def readGEBTOut(fn_gebt_out, beam): + flag_analysis = beam.analysis_type + if flag_analysis <= 1: + return readGEBTOutStatic(fn_gebt_out, beam) + elif flag_analysis == 3: + return readGEBTOutEigen(fn_gebt_out, beam) + + + + + + + + + + + +# ==================================================================== +# +# +# WRITERS +# +# +# ==================================================================== + +def writeGEBTIn(beam, fn_gebt_in=''): + """ Write data to the GEBT input""" + + if fn_gebt_in == '': + fn_gebt_in = beam.name + '.dat' + + with open(fn_gebt_in, 'w') as fout: + utl.writeFormatIntegers( + fout, [beam.analysis_type, beam.max_iteration, beam.num_steps], '8d' + ) + if beam.analysis_type != 0: + utl.writeFormatFloats(fout, beam.angular_velocity[0]) + utl.writeFormatIntegers(fout, beam.angular_velocity[1]) + utl.writeFormatFloats(fout, beam.linear_velocity[0]) + utl.writeFormatIntegers(fout, beam.linear_velocity[1]) + if beam.analysis_type == 3: + fout.write('{0:8d}\n'.format(beam.num_eigens)) + fout.write('\n') + + line = [ + len(beam.points), len(beam.segments), len(beam.pconditions), + len(beam.sections), len(beam.frames), len(beam.mconditions), + len(beam.distrloads), len(beam.timefunctions), len(beam.initcurvatures) + ] + utl.writeFormatIntegers(fout, line, '4d') + fout.write('\n') + + # Write the block of points + for pid, coord in beam.points.items(): + fout.write('{0:4d}'.format(pid)) + utl.writeFormatFloats(fout, coord) + fout.write('\n') + + # Write the block of members + # for i, member in enumerate(members): + for bsid, bs in beam.segments.items(): + line = [ + bsid, bs.points[0], bs.points[1], bs.css[0], bs.css[1], + bs.frame_id, bs.num_divisions, bs.curv_id + ] + utl.writeFormatIntegers(fout, line, '4d') + fout.write('\n') + + # Write the block of point conditions + for condition in beam.pconditions: + fout.write('{0:4d}\n'.format(condition['point'])) + utl.writeFormatIntegers(fout, condition['components'], '16d') + utl.writeFormatFloats(fout, condition['values']) + utl.writeFormatIntegers( + fout, np.zeros(6, dtype=np.int8), '16d') + utl.writeFormatIntegers( + fout, np.zeros(6, dtype=np.int8), '16d') + fout.write('\n') + + # Write the block of cross sectional properties + # for i, member in enumerate(members): + for sid, ms in beam.sections.items(): + fout.write('{0:4d}\n'.format(sid)) + # tcm = st['tcm'] + # for row in tcm: + # utl.writeFormatFloats(fout, row) + if len(ms.eff_props[1]['compliance']['refined']) > 0: + for row in ms.eff_props[1]['compliance']['refined']: + utl.writeFormatFloats(fout, row) + else: + for i, row in enumerate(ms.eff_props[1]['compliance']['classical']): + utl.writeFormatFloats(fout, (row[0], 0, 0, row[1], row[2], row[3])) + if i == 0: + utl.writeFormatFloats(fout, np.zeros(6)) + utl.writeFormatFloats(fout, np.zeros(6)) + fout.write('\n') + if beam.analysis_type != 0: + # mass = st['mass'] + # for row in mass: + # utl.writeFormatFloats(fout, row) + for row in ms.eff_props[1]['mass']['origin']: + utl.writeFormatFloats(fout, row) + fout.write('\n') + + # Write the block of self.frames + for i, cij in beam.frames.items(): + fout.write('{0:4d}\n'.format(i)) + for row in cij: + utl.writeFormatFloats(fout, row) + fout.write('\n') + + return fn_gebt_in diff --git a/ms/prebeam.py b/ms/prebeam.py new file mode 100644 index 0000000..8d08dd8 --- /dev/null +++ b/ms/prebeam.py @@ -0,0 +1,475 @@ +import os +import numpy as np +import xml.etree.ElementTree as et +import msgpi.ms.beam as msb +import msgpi.sg as sgi +import msgpi.cross_section as sgcs +import msgpi.utils as utl + + +def preBeam(fn_beam, mode=1, sections=[]): + """ Preprocessor of GEBT + + :param fn_beam: File name of the PreGEBT main input file (.xml) + :type fn_beam: string + + :param mode: Mode of running + 1 - create SG (cross-section) input files + 2 - run SGs + 3 - create global 1D beam input file + :type mode: int + """ + + beam = msb.Beam() + + + # ---------------------------------------------------------------- + # Parse XML parametric inputs + tree = et.parse(fn_beam) + xr_beam = tree.getroot() + # beam.name = xr_beam.get('name') + beam.name = os.path.splitext(os.path.basename(fn_beam))[0] + + # self.fn_gebt_in = self.name + '_gebt.dat' + + # Read analysis settings + xe_analysis = xr_beam.find('analysis') + beam.analysis_type = int(xe_analysis.find('type').text) + beam.max_iteration = int(xe_analysis.find('max_iteration').text) + beam.num_steps = int(xe_analysis.find('num_steps').text) + if beam.analysis_type == 3: + beam.num_eigens = int(xe_analysis.find('num_eigenvalue').text) + + + # ---------------------------------------------------------------- + # Read design + xe_design = xr_beam.find('design') + inpfmt = 1 + xa_inpfmt = xe_design.get('format') + if xa_inpfmt is not None: + inpfmt = int(xa_inpfmt) + + # Read global frame a + xe_frame_a = xe_design.find('global_frame_a') + frame_a = list(map(float, xe_frame_a.text.strip().split())) + frame_a = np.array(frame_a).reshape((3, 3)) + # print(frame_a) + + if inpfmt == 1: + # Read points + pn2l = {} # point name to label dictionary + plabel = 0 + for xe_point in xe_design.findall('point'): + plabel += 1 + pname = xe_point.get('name') + pcoord = list(map(float, xe_point.text.strip().split())) + pn2l[pname] = plabel + beam.points[plabel] = pcoord + + # Read beam segments + bsn2l = {} # beam segment name to label dictionary + bslabel = 0 + fblabel = 0 # local frame b label + for xe_sgm in xe_design.findall('segment'): + bslabel += 1 + bs = msb.BeamSegment() + + # Start station + xe_start = xe_sgm.find('start') + pn = xe_start.find('location').text.strip() + bs.points.append(pn2l[pn]) + pc = beam.findPtCoordByName(pn2l[pn]) + bs.coords.append(pc) + xe_cs = xe_start.find('cross_section') + csname = xe_cs.get('name') + cslabel = beam.findSectionByName(csname) + if cslabel == 0: + cslabel = len(beam.sections) + 1 + mso = sgi.MaterialSection(1) + mso.name = csname + for s in sections: + if s.name == csname: + mso = s + beam.sections[cslabel] = mso + if xe_start.find('rotate_a1') is not None: + bs.rotate_a1 = float(xe_start.find('rotate_a1').text.strip()) + bs.css.append(cslabel) + # print(beam.sections) + + # Stop station + xe_stop = xe_sgm.find('stop') + pn = xe_stop.find('location').text.strip() + bs.points.append(pn2l[pn]) + pc = beam.findPtCoordByName(pn2l[pn]) + bs.coords.append(pc) + xe_cs = xe_stop.find('cross_section') + csname = xe_cs.get('name') + # print(csname) + cslabel = beam.findSectionByName(csname) + # print(cslabel) + if cslabel == 0: + cslabel = len(beam.sections) + 1 + mso = sgi.MaterialSection(1) + mso.name = csname + for s in sections: + if s.name == csname: + mso = s + beam.sections[cslabel] = mso + bs.css.append(cslabel) + # print(bs.css) + + nd = int(xe_sgm.find('mesh_size').text.strip()) + bs.num_divisions = nd + + # Local frame + if xe_sgm.find('local_b') is not None: + xe_b = xe_sgm.find('local_b') + if xe_b.get('method').strip() != 'global': + p12 = np.array(list(map(float, xe_b.find('p12').text.strip().split()))) + p0 = np.array(beam.points[bs.points[0]]) + p1 = np.array(beam.points[bs.points[1]]) + b1 = p1 - p0 + b1 = b1 / np.linalg.norm(b1) + b12 = p12 - p0 + b3 = np.cross(b1, b12) + b3 = b3 / np.linalg.norm(b3) + b2 = np.cross(b3, b1) + frame_b = np.vstack((b1, b2, b3)) + # print(frame_b) + # Calculate C_ij^ab + cij = np.zeros((3, 3)) + for i in range(3): + for j in range(3): + cij[i, j] = np.dot(frame_a[i], frame_b[j]) + # print(cij) + fblabel += 1 + bs.frame_id = fblabel + beam.frames[fblabel] = cij + # for i in range(3): + # a_temp = cij[i, 0] * frame_b[0] + cij[i, 1] * frame_b[1] + cij[i, 2] * frame_b[2] + # print(a_temp) + + # Twist + twist = 0.0 + if xe_sgm.find('twist') is not None: + twist = float(xe_sgm.find('twist').text.strip()) + bs.twist = twist + + beam.segments[bslabel] = bs + + if beam.analysis_type != 0: + ws = float(xe_design.find('angular_velocity').text.strip()) + vspname = xe_design.find('linear_velocity').get('at') + vsplabel = pn2l[vspname] + vs = ws * beam.points[vsplabel][0] + ws = [0., 0., ws] + wtf = [0, 0, 0] + beam.angular_velocity = [ws, wtf] + vs = [0, vs, 0] + vtf = [0, 0, 0] + beam.linear_velocity = [vs, vtf] + + # Read point conditions + # conditions = [] + for condition in xe_design.findall('condition'): + pname = condition.find('location').text + components = list(map(int, condition.find( + 'components').text.strip().split())) + values = list(map(float, condition.find('values').text.strip().split())) + beam.pconditions.append({ + 'point': pn2l[pname], + 'components': components, + 'values': values + }) + elif inpfmt == 2: + fn_cs_base = None + fn_bsl_base = None + fn_lyp_base = None + + xe_templates = xe_design.find('templates') + if xe_templates is not None: + xe_template_cs = xe_templates.find('cross_section') + if xe_template_cs is not None: + fn_cs_base = xe_template_cs.text.strip() + + xe_template_bsl = xe_templates.find('baselines') + if xe_template_bsl is not None: + fn_bsl_base = xe_template_bsl.text.strip() + + xe_template_lyp = xe_templates.find('layups') + if xe_template_lyp is not None: + fn_lyp_base = xe_template_lyp.text.strip() + + length = float(xe_design.find('length').text.strip()) + + xe_layers = xe_design.find('layers') + + xe_domain = xe_layers.find('domain') + domain_name = xe_domain.text.strip() + + layers = [] + kps = [] + + print(' - reading layers...') + for xe_layer in xe_layers.findall('layer'): + layer = {} + layer['lamina'] = xe_layer.find('lamina').text.strip() + + layer['spanrange'] = [0, length] + xe_spanrange = xe_layer.find('spanrange') + if xe_spanrange is not None: + defby = 'normalized_ends' + xa_defby = xe_spanrange.get('define') + if xa_defby is not None: + defby = xa_defby + if defby == 'normalized_ends': + [start, stop] = list(map(float, xe_spanrange.text.strip().split())) + layer['spanrange'] = [start*length, stop*length] + + layer['angle'] = 0.0 + xe_angle = xe_layer.find('angle') + if xe_angle is not None: + func = 'constant' + xa_func = xe_angle.get('function') + if xa_func is not None: + func = xa_func + if func == 'constant': + layer['angle'] = float(xe_angle.text.strip()) + elif func == 'polynomial': + xe_order = xe_angle.find('order') + xe_coeff = xe_angle.find('coefficients') + xe_ndivs = xe_angle.find('divisions') + + porder = int(xe_order.text.strip()) + coeffs = list(map(float, xe_coeff.text.strip().split())) + ndivs = int(xe_ndivs.text.strip()) + + layer['angle_func'] = utl.Polynomial2DSP(porder, coeffs) + + layer['kps'] = np.linspace( + layer['spanrange'][0], + layer['spanrange'][1], + ndivs+1 + ) + pass + + # kps.append(layer['spanrange'][0]) + # kps.append(layer['spanrange'][1]) + kps.append(layer['kps']) + + # layer['sectionrange'] = [] + # xe_sectionrange = xe_layer.find('sectionrange') + # if xe_sectionrange is not None: + # defby = 'domain_list' + # xa_defby = xe_sectionrange.get('define') + # if xa_defby is not None: + # defby = xa_defby + # if defby == 'domain_list': + # # + + layers.append(layer) + + # for layer in layers: + # print(layer) + + kps = np.concatenate(kps).tolist() + kps = list(set(kps)) + print('kps:', kps) + + for i, kp in enumerate(kps): + xe_point = et.SubElement(xe_design, 'point') + xe_point.set('name', 'p'+str(i)) + xe_point.text = str(kp) + ' 0 0' + + # Generate beam segments + print(' - creating beam segments...') + bm_sgms = [] + layups = [] + for i, kp in enumerate(kps): + if i == 0: + # first station + bm_sgm = {} + bm_sgm['loc_start'] = 'p' + str(i) + layup = [] + for j, layer in enumerate(layers): + if kp == layer['spanrange'][0]: + layup.append(layer) + bm_sgm['layup_start'] = layup + bm_sgms.append(bm_sgm) + elif i == len(kps) - 1: + # last station + layup = [] + bm_sgms[-1]['loc_stop'] = 'p' + str(i) + for j, layer in enumerate(layers): + if kp == layer['spanrange'][1]: + layup.append(layer) + bm_sgms[-1]['layup_stop'] = layup + else: + # inner stations + layup = [] + bm_sgms[-1]['loc_stop'] = 'p' + str(i) + for j, layer in enumerate(layers): + if (kp > layer['spanrange'][0]) and (kp <= layer['spanrange'][1]): + layup.append(layer) + bm_sgms[-1]['layup_stop'] = layup + + bm_sgm = {} + bm_sgm['loc_start'] = 'p' + str(i) + layup = [] + for j, layer in enumerate(layers): + if (kp >= layer['spanrange'][0]) and (kp < layer['spanrange'][1]): + layup.append(layer) + bm_sgm['layup_start'] = layup + bm_sgms.append(bm_sgm) + + # for i, bs in enumerate(bm_sgms): + # print('') + # print('- beam segment {0}'.format(i)) + # print(' - start layup') + # for j, layer in enumerate(bs['layup_start']): + # print(j, layer) + # print(' - stop layup') + # for j, layer in enumerate(bs['layup_stop']): + # print(j, layer) + + # Write PreVABS input files + print(' - writing sg input files...') + fn_lyp_tmp = fn_lyp_base + '.xml' + fn_lyp = beam.name + '_lyps.xml' + + xt_layups = et.parse(fn_lyp_tmp) + xr_layups = xt_layups.getroot() + for i, bs in enumerate(bm_sgms): + # start layup + layup_name = '_'.join(('lyp', str(i+1), '0')) + xe_layup = et.SubElement(xr_layups, 'layup') + xe_layup.set('name', layup_name) + for j, lyr in enumerate(bs['layup_start']): + xe_layer = et.SubElement(xe_layup, 'layer') + xe_layer.set('lamina', lyr['lamina']) + + f = lyr['angle_func'] + x0, x1 = lyr['spanrange'] + y = (kps[i] - x0) / (x1 - x0) + # xe_layer.text = str(lyr['angle']) + xe_layer.text = str(f([y, 0])) + + # start cross-section + fn_cs = '_'.join((beam.name, 'cs', str(i+1), '0')) + xt_cs = et.parse(fn_cs_base+'.xml') + xr_cs = xt_cs.getroot() + xr_cs.set('name', fn_cs) + xe_include = xr_cs.find('include') + xe_include_layup = xe_include.find('layup') + xe_include_layup.text = fn_lyp[:-4] + + xe_cmp = None + for xe in xr_cs.findall('component'): + if xe.get('name') == domain_name: + xe_cmp = xe + break + if xe_cmp is None: + xe_cmp = et.SubElement(xr_cs, 'component') + xe_sgm = xe_cmp.find('segment') + if xe_sgm is None: + xe_sgm = et.SubElement(xe_cmp, 'segment') + xe_lyp = xe_sgm.find('layup') + if xe_lyp is None: + xe_lyp = et.SubElement(xe_sgm, 'layup') + xe_lyp.text = layup_name + + bs['cs_start'] = fn_cs + xt_cs.write(fn_cs+'.xml') + + # stop layup + layup_name = '_'.join(('lyp', str(i+1), '1')) + xe_layup = et.SubElement(xr_layups, 'layup') + xe_layup.set('name', layup_name) + for j, lyr in enumerate(bs['layup_stop']): + xe_layer = et.SubElement(xe_layup, 'layer') + xe_layer.set('lamina', lyr['lamina']) + + f = lyr['angle_func'] + x0, x1 = lyr['spanrange'] + y = (kps[i+1] - x0) / (x1 - x0) + # xe_layer.text = str(lyr['angle']) + xe_layer.text = str(f([y, 0])) + + # stop cross-section + fn_cs = '_'.join((beam.name, 'cs', str(i+1), '1')) + xt_cs = et.parse(fn_cs_base+'.xml') + xr_cs = xt_cs.getroot() + xr_cs.set('name', fn_cs) + xe_include = xr_cs.find('include') + xe_include_layup = xe_include.find('layup') + xe_include_layup.text = fn_lyp[:-4] + + xe_cmp = None + for xe in xr_cs.findall('component'): + if xe.get('name') == domain_name: + xe_cmp = xe + break + if xe_cmp is None: + xe_cmp = et.SubElement(xr_cs, 'component') + xe_sgm = xe_cmp.find('segment') + if xe_sgm is None: + xe_sgm = et.SubElement(xe_cmp, 'segment') + xe_lyp = xe_sgm.find('layup') + if xe_lyp is None: + xe_lyp = et.SubElement(xe_sgm, 'layup') + xe_lyp.text = layup_name + + bs['cs_stop'] = fn_cs + xt_cs.write(fn_cs+'.xml') + + # for i, bs in enumerate(bm_sgms): + # print('') + # print('- beam segment {0}'.format(i)) + # print(bs) + + xt_layups.write(fn_lyp) + + # Write beam segments + fn_beam_sgn = fn_beam[:-4] + '_sgn.xml' + print(' - writing beam segments in file ' + fn_beam_sgn) + for i, bs in enumerate(bm_sgms): + xe_bs = et.SubElement(xe_design, 'segment') + xe_start = et.SubElement(xe_bs, 'start') + xe_loc = et.SubElement(xe_start, 'location') + xe_loc.text = bs['loc_start'] + xe_cs = et.SubElement(xe_start, 'cross_section') + xe_cs.set('name', bs['cs_start']) + + xe_stop = et.SubElement(xe_bs, 'stop') + xe_loc = et.SubElement(xe_stop, 'location') + xe_loc.text = bs['loc_stop'] + xe_cs = et.SubElement(xe_stop, 'cross_section') + xe_cs.set('name', bs['cs_stop']) + + xe_design.set('format', '1') + xe_design.remove(xe_templates) + xe_design.remove(xe_layers) + tree.write(fn_beam_sgn) + + # Calculate C_ab + # self.frames = [] + frame_a = [ + np.array([1., 0., 0.]), + np.array([0., 1., 0.]), + np.array([0., 0., 1.]) + ] + # for member in members: + # for i in range(1, self.num_points): + # frame_b = self.stations[i]['frame_b'] + # Cab = utl.calcCab(frame_a, frame_b) + # self.frames.append(Cab) + # self.num_frames = len(self.frames) + + # ms = float(design.find('mesh_size').text) + # # divisions = [] + # for i in range(1, self.num_points): + # dx1 = self.stations[i]['coordinates'][0] - \ + # self.stations[i - 1]['coordinates'][0] + # self.divisions.append(int(np.ceil(dx1 / ms))) + + return beam diff --git a/presg.py b/presg.py new file mode 100644 index 0000000..47309aa --- /dev/null +++ b/presg.py @@ -0,0 +1,329 @@ +import os +import sys +import subprocess as sbp +import numpy as np +import xml.etree.ElementTree as et +import msgpi.sg as sgm +import msgpi.io.iosc as miosc +import msgpi.io.iovabs as miovabs +import msgpi.utils as utl + + +def readMaterialFromXMLElement(xem): + type_flag = { + 'isotropic': 0, + 'orthotropic': 1, + 'anisotropic': 2 + } + elastic_label = { + 0: ('e', 'nu'), + 1: ( + 'e1', 'e2', 'e3', + 'g12', 'g13', 'g23', + 'nu12', 'nu13', 'nu23' + )#, + # 'anisotropic': ( + # ('c11', 'c12', 'c13', 'c14', 'c15', 'c16'), + # ('c12', 'c22', 'c23', 'c24', 'c25', 'c26'), + # ('c13', 'c23', 'c33', 'c34', 'c35', 'c36'), + # ('c14', 'c24', 'c34', 'c44', 'c45', 'c46'), + # ('c15', 'c25', 'c35', 'c45', 'c55', 'c56'), + # ('c16', 'c26', 'c36', 'c46', 'c56', 'c66'), + # ) + } + + m = sgm.MaterialSection() + m.name = xem.get('name').strip() + + m.eff_props[3]['type'] = type_flag[xem.get('type').strip()] + m.eff_props[3]['density'] = float(xem.find('density').text.strip()) + + ep = xem.find('elastic') + if (m.eff_props[3]['type'] == 0) or (m.eff_props[3]['type'] == 1): + ep_labels = elastic_label[m.eff_props[3]['type']] + for tag in ep_labels: + m.eff_props[3]['constants'][tag] = float( + ep.find(tag).text.strip() + ) + else: + stf = np.zeros((6, 6)) + for i in range(6): + rowtag = 'row' + str(i) + row = list(map(float, ep.find(rowtag).text.strip().split())) + stf[i, i:] = row + for i in range(1, 6): + for j in range(0, i): + stf[i, j] = stf[j, i] + m.eff_props[3]['stiffness'] = stf + + if xem.find('strength') is not None: + xe_str = xem.find('strength') + if xe_str.find('criterion') is not None: + m.strength['criterion'] = int(xe_str.find('criterion').text.strip()) + m.strength['constants'] = list(map(float, xe_str.find('constants').text.strip().split())) + + return m + + +def preSG1D(xr_sg, smdim): + # debug = True + # print('running presc...') + + # ==================================================================== + # Parse XML parametric inputs + # tree = et.parse(fn_xml) + # xr_sg = sg_xml_tree.getroot() + name = xr_sg.get('name') + # sgdim = int(xr_sg.get('dim').text) + # smdim = int(xr_sg.get('model_dim').text) + + sg = sgm.StructureGene(name, 1, smdim) + + # Head + sg.analysis = int(xr_sg.find('analysis').text) + sg.model = int(xr_sg.find('model').text) # model (0: classical, 1: shear refined) + sg.degen_element = int(xr_sg.find('elem_flag').text) + sg.trans_element = int(xr_sg.find('trans_flag').text) + sg.uniform_temperature = int(xr_sg.find('temp_flag').text) + + if smdim == 2: + if xr_sg.find('init_curv') is not None: + ic = xr_sg.find('init_curv').text.strip() + sg.initial_curvature = list(map(float, ic.split())) + + # Materials and Laminae + materials = {} + d_lamina = {} + + if xr_sg.find('include') is not None: + xe_include = xr_sg.find('include') + if xe_include.find('material') is not None: + fn_material = xe_include.find('material').text.strip() + '.xml' + xt = et.parse(fn_material) + xe_materials = xt.getroot() + for xem in xe_materials.findall('material'): + m = readMaterialFromXMLElement(xem) + materials[m.name] = m + for lmn in xe_materials.findall('lamina'): + ln = lmn.get('name').strip() + lm = lmn.find('material').text.strip() + lt = float(lmn.find('thickness').text.strip()) + d_lamina[ln] = {'material': lm, 'thickness': lt} + + if xr_sg.find('materials') is not None: + xe_materials = xr_sg.find('materials') + for xem in xe_materials.findall('material'): + m = readMaterialFromXMLElement(xem) + materials[m.name] = m + for lmn in xe_materials.findall('lamina'): + ln = lmn.get('name').strip() + lm = lmn.find('material').text.strip() + lt = float(lmn.find('thickness').text.strip()) + d_lamina[ln] = {'material': lm, 'thickness': lt} + + # Layup + # List of material, thickness, angle + # [{'material': 'name', 'thickness': t, 'angle': a, 'combo': id}, {}, ...] + ld_layer = [] + l_material_used = [] + layup = xr_sg.find('layup') + method = 'll' # ll: layer list, ss: stacking sequence, pp: ply percentage + if layup.get('method') is not None: + method = layup.get('method') + if method == 'll': + # layer list + for lyr in layup.findall('layer'): + ln = lyr.get('lamina').strip() + lm = d_lamina[ln]['material'] + # if lm not in l_material_used: + # l_material_used.append(lm) + # d_material[lm]['id'] = len(l_material_used) + lt = d_lamina[ln]['thickness'] + nl = 1 + lon = lyr.text.strip().split(':') + lo = float(lon[0].strip()) + if len(lon) > 1: + nl = int(lon[1].strip()) + for ni in range(nl): + ld_layer.append( + {'material': lm, 'thickness': lt, 'angle': lo} + ) + elif method == 'ss': + # stacking sequence + ln = layup.find('lamina').text.strip() + lm = d_lamina[ln]['material'] + # l_material_used.append(lm) + # d_material[lm]['id'] = len(l_material_used) + lt = d_lamina[ln]['thickness'] + layup_code = layup.find('code').text.strip() + code = utl.parseLayupCode(layup_code) + for lo in code: + ld_layer.append({'material': lm, 'thickness': lt, 'angle': lo}) + elif method == 'pp': + # ply percentage + pass + # print code + + # Record used materials and create material-orientation combinations + mid = 0 + cid = 0 + tt = 0.0 # total thickness + for layer in ld_layer: + # Check used material-orientation combination first + cid = sg.findComboByMaterialOrientation(layer['material'], layer['angle']) + if cid == 0: + # Not found + # Check used material + mid = sg.findMaterialByName(layer['material']) + if mid == 0: + # Not found + mid = len(sg.materials) + 1 + sg.materials[mid] = materials[layer['material']] + cid = len(sg.mocombos) + 1 + sg.mocombos[cid] = [mid, layer['angle']] + sg.prop_elem[cid] = [] + layer['mocombo'] = cid + tt = tt + layer['thickness'] + + # ---------------------------------------------------------------- + # Meshing + pan = 0.0 + if smdim == 3: + ens = 2 # by default, use 2-node element for 3D structure + sg.omega = tt + elif smdim == 2: + ens = 5 # by default, use 5-node element for 2D structure + + if xr_sg.find('translate') is not None: + pan = float(xr_sg.find('translate').text.strip()) + + gms = float(xr_sg.find('global_mesh_size').text.strip()) + + if xr_sg.find('element_nodes') is not None: + ens = int(xr_sg.find('element_nodes').text.strip()) + + nply = len(ld_layer) + # tthk = thickness * nply + ht = tt / 2.0 + # print tthk + + # nodes_major = np.array([-ht + pan, ]) + # nid = 1 + nid1 = 1 + eid = 0 + yprev = -ht + pan + sg.nodes[nid1] = [yprev, ] + for lyr in ld_layer: + t = lyr['thickness'] + ne = int(round(t / gms)) # number of element for this layer + # lyr['nelem'] = ne + ns = np.linspace(yprev, yprev+t, ne+1) + for i in range(ne): + # nid = nid + ens - 1 + nid2 = nid1 + ens - 1 + sg.nodes[nid2] = [ns[i+1], ] + eid = eid + 1 + sg.elements[eid] = [nid1, nid2] + sg.elementids1d.append(eid) + sg.elem_prop[eid] = lyr['mocombo'] + sg.prop_elem[lyr['mocombo']].append(eid) + nid1 = nid2 + yprev = yprev + t + + # Change the order of each element + if ens > 2: + for eid in sg.elements.keys(): + nid1 = sg.elements[eid][0] + nid2 = sg.elements[eid][1] + nq0y3 = sg.nodes[nid1][0] + nq4y3 = sg.nodes[nid2][0] + nid3 = nid1 + 1 + sg.elements[eid].append(nid3) + if ens == 4: + pass + else: + nq2y3 = (nq0y3 + nq4y3) / 2.0 + nq1y3 = (nq0y3 + nq2y3) / 2.0 + nq3y3 = (nq2y3 + nq4y3) / 2.0 + if ens == 5: + # nid = nid + 1 + nid5 = nid3 + 1 + nid4 = nid5 + 1 + sg.nodes[nid3] = [nq1y3, ] + sg.nodes[nid5] = [nq2y3, ] + sg.nodes[nid4] = [nq3y3, ] + # sg.elements[eid].append(nid3) + sg.elements[eid].append(nid4) + sg.elements[eid].append(nid5) + else: + sg.nodes[nid3] = [nq2y3, ] + + return sg + + +def preSG(sg_xml, analysis, solver='swiftcomp', write_input=True, scrnout=True): + """ Preprocessor of a structure gene + + Parameters + ---------- + sg_xml : str + File name of SG design parameters (XML format) + analysis : str + The analysis that will be carried out after the generation of + the input file. + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + Format of the generated input file ('vabs' or 'swiftcomp') + """ + if scrnout: + print(' - preprocessing structure gene...') + # fn_head, fn_tail = os.path.split(xml_sg) + fn_base, fn_extn = os.path.splitext(sg_xml) + fn_sg_in = fn_base + '.sg' + + xt = et.parse(sg_xml) + xr_sg = xt.getroot() + + if xr_sg.tag == 'cross_section': + smdim = 1 + cmd = ['prevabs', '-i', sg_xml] + if 'vabs' in solver: + cmd.append('-vabs') + elif 'swiftcomp' in solver: + cmd.append('-sc') + cmd.append('-' + analysis) + + cmd = ' '.join(cmd) + try: + FNULL = open(os.devnull, 'w') + sbp.call(cmd, stdout=FNULL, stderr=sbp.STDOUT) + return fn_sg_in, smdim + except: + e = sys.exc_info()[0] + print('[ERROR] ', e) + return + else: + try: + sgdim = int(xr_sg.get('sgdim')) + except TypeError: + print('[ERROR] Unable to read structure gene dimension (sgdim):', xr_sg.get('sgdim')) + return + + try: + smdim = int(xr_sg.get('smdim')) + except TypeError: + print('[ERROR] Unable to read structural model dimension (smdim):', xr_sg.get('smdim')) + return + + if sgdim == 1: + sg = preSG1D(xr_sg, smdim) + + if write_input: + miosc.writeSCIn(sg, fn_sg_in) + return fn_sg_in, smdim + else: + return sg diff --git a/presg.pyc b/presg.pyc new file mode 100644 index 0000000000000000000000000000000000000000..581cae78a2f87aff686c04559a0ade00d4460489 GIT binary patch literal 8074 zcmc&(O>7%ka;|PtltfYg)(>Tk-S*fPHDk$^b};|<2^S+qRRtRgMbVmGO# zo844*(*u?_55~+ z_6Gs7h$ub2S5>cGy{dZks-AK`n=bm58(*(c`jf-=d-#jKhsMUQMjp}0ps&%1mcey8 z$z}aKos4Du0yW2}IYB3r8Kg+fDLN@-!)a>HP;-_}<}%2QgmvF?sQp*XC^P=k^wVS#pnCzn%xQ;O3+YsZ&=ryVnJsLi#3s;{#g zSoq~6h2Pdw{R=R`o=Rel^${9RkLRzz;o!<0z=u~@AbXu2`!u2k#bxD}sv(`WBGWK9 zsu*zC=PnvnKioD7;xK5tHJ7~_`X>s)_#y+su@%tr2b3{WIx(aCj(j;g>T6eh?Zr9Z z8aun`cRJqPi}NxEh_sVIGL+SF1Qqw%4h|I79!Fz4e&ohI={uek#coYxuLXV-TYen5 zd^lyXFnBE}(A4X`6Gia^P*EH@ej~n6vqCp^LO1X$yaSV0FD4e?IisFvl zJB@a?LML*puy%18MzxD-EwG*Fv=uel?r9sS#&)|GuYphKSazoU55u7K=~te8;5jYF zM>v6ahiFL4Y2)}5b$sijkyb2~v{n6%wyDo*+uD*g4T_!IthT7lXtLSMiiZcG-y4-6xw;bvb)e%1t;n5=I=Q!p=O-z>O`v7E!=<1MH z5lRQB{8iNW2ukROvnEDRDgX}{B5Dp8j$dVn0Eq}SAwhmZOXZQ6YOAny6^P}UCAqWm zSqPeys6RK1yNEH>dAgj?5T;z*Ze!A_OA4DUvH&Ub1BrosQwuRhkjtKmfah z>i~ORpyr~4Ybbh>`nXT7iF`%sKIHh13tq81JiE!kJYJ^0uvZZn*B8fOM90fN#jTKJ zQ>Iyh!!wtXY3d`D$c&lpenxRLi=mfdgJSEqt zzbHzU3@G)EujNL8pkeh~V|q?WKFW9hNMiyRSdyc|vwx&yiTX>zNDMf%mZ-T5#g}0t zrSvlNUxrO^CKEUXre)6E1Q1>X=~l#^$P1UKJhbB~#jDhpD;JYr<|&-Y^9chTGK3-B zGE3aur?^ag$+`~DDwJS=OmvOY#W#EbY3DNSV_ofXHGh2B$sPNDEjwP39sg&Zw+qiw zd&{Iv;Ecawd+8j>sRlbt)+oMC&2{Rp$v&@B_kHSLhsQ%e?CsYJ__z6EADVO=C5TcBh5+x1O1Zl`TH?WSZ-ug zYy?OtOxfyxFvh9}g`|v7_X?wNh!s;wjj>RgBq>ru&JvSiAyBJ@8#s#KWANIK;hsFW zB?ERdCFXTxa!M4q$_ln~IX8C98^bId&vlWQ3W%Qxrh*co?>6~x!a-17Dj@L(mYJ$( ze?!Dk(eKjnm2K6_np|Vd9-=)XhpDL67M|#@i7*+U+iSJOqLjBvc1nGGL zn-t%rCb$U~;8yyZyvz#UfbUTE&*W%qGVS5nA5!uzFVMe5E1+3T?oc1-!%0m(E z!zlLzn1D0r?^9|^ta1-X7_!0#3|@Hc=r2Lv#e0+&Y!?^AM@F89WdSzZW> z&0qdo;TQt(4T)pEKG}UwV@4a)e21@@!?XJuOCcbvSLQ$gCL;hu1l|cKM16{4D+An5 zr)=|1#GN0EOlH|0s7*kI+v0A&H*xvn@je)ALrrWrTyq#bOh-H19pH9Zn{YrrI%AXI zB}?p4e^2!uQGaK|x(}eeha=-XjCV)IkCb*r&!V_TO1Q!HS>Qc^v5=#zz#1UF3ZROk z%G=14o=fhlst(SwE#LBb5fWEUj$6ph$}8O{kM-913Y137pH$jLobBcTT26**NE2JJTDjrtK8TR|;ood;6xxI9$} z!#njLay$o6!wb$WuiA2=i)!S4=A@aqB(znZli-Z!doC|lc^8dY=4J(rkRQcirxtfY z=c(g6<_6>5;^qc7oWq*?XoS>u8jfFeRKhF}rT7fOnLHBodJwj%Yy&5>n6Y_}xp4Yz zvW3k%+^lnR7Y&MAWgZknVj-EMY6(>mwXHC6L?@_c!XSvH&lM2L0P_(unqfvcG@j#| zw|Kyrx7BVtD3ci8i9Ng!5DBkmY_%t>c-|y4dcvxHzn7H;ABQeV0!xasg_I=2z2AF| za@727fR+C!J4Iq1F4aUqV(g&hG%V?fEuA{mMr9_>%_j%dR(2S;h^_n8 z%)web?SmgkZvFt%s`dm7LwQxK7_SBJ&~h=Cg!3yF%+Wb}`I110Y`wA(vMz3&AO`H$rZFhTEVWP zjxC|4Eot*;8GA#!Ex4<&n`4**+5%`6w35!WGx-^eS21rvE2Ba#Xd8M-U)5%EGl0&b zcHRVKF3&Ry+O)o=m2!)KECRczt)QYVY74-YfGy!OiFI%3H>54;MOkqQRMVha*5A?= zb>>@?T6|u5GdY%hKK&H6Yv7Ut6@S11v~&E?uW2j#iY_NBdhOO@MB8?aZO6VksM~)6 zjo}sa>6g{REEnv=SJ^2t)I24~F_C(67Af^P@^FPtt5y`-C~@Sav031fp0C`uH9P~O zEL9J|pl=}W!F%U8zXiN3Q1}h0H}NPd@0*&BJQqc*!u@52lZu(EX)b3$i?>L;g>q_x zA{Z}(0Am2Yy~liWu=rR5MFIpz!8^qzJ7|!nq4Kp+q5{GtgsySjFv%~ON5_{x*6@SC zsh&{nIixLkOPA+J5Ip@d9;c5!^?yh`Jjl**#oe4w3ol-L8qc}+D35adWsyRgd7*?F zr5ZB1!a?4wnZx824&Mc)K%_XXhyVcP*C3eH7Q& zBv2!J1?-4-Q{js-dK;0%_-g1Np4aeN7=%VpH!K4WQ|hiXa1lDiqJf{US|LAA;&Dob zTZ11NMvW>`3gdvM3V^44=z5OfV-|BgdTQXFb{oEdyD8;lRPZom@I7k9n=&m+gWllu z;=(ZosQAK)jn}T{8Rw2svqD^Awh?qU6(BL*#tjY}Dij#zpi4`W=v4azzs+aAC# zj4k6LXa$(%ekQX9C>to-xr%(?1vSe{DIbJREqLV&H(kdJ&cblb8+FU|kQlI16dRU`HKt5V%G zQu9&~-(p*?+zT@13u%jc=AWP^oVKEx+)i848q&U&oF|HHm^|w`a01Msh*Rm`{lGLo zHqVUbDp6=4QHZ`Y0{_If&f%DGVDK~P8+$N7&F}5-oJzz~HxJAPb3XX>9AI(D+^_K9 zJZ1BsvW)3+730;Kqi@^7dpt`kZ_{|(Hcbr7BW_-Bv(3#h8X<26c(CSq@DE?|Lfd)e zcaV*K$Rdb&=*31impd6;y5Zr=YIs?@QJpD zm5=;?VA~Dw3%l~Ea|3Y~9-H4|DcHOtWPAw4x^5s{KDF=)>2vz5NIS9Xsd`m<=BL1@ zGDTd8z4b-dam>F#XJ|X={h$4_V* [{0:%H}:{0:%M}] eval {1:6d} {2:s}'.format( + dt.datetime.now(), evid, ' '.join(argv) + )) + + +def interface(params, results, tmp_map, res_map, sg_xml, analysis, solver, scrnout=True): + """ Create a workflow for a design problem using Dakota + + Parameters + ---------- + params : Parameters + Dakota Parameters object + results : Results + Dakota Results object + tmp_map : dict + Dictionary of file names that will be used in dprepro + {'template name': 'true name', ...} + res_map : dict + Disctionary of keywords for results of Dakota + {'dakota_response_keyword': 'sg_result_keyword', ...} + sg_xml : str + Main file name for preprocessor + analysis : str + The analysis that will be carried out after the generation of + the input file. + h - homogenization + d - dehomogenization/localization/recover + f - initial failure strength + fe - initial failure envelope + fi - initial failure indices and strength ratios + solver : str + The name of the executable ('vabs' or 'swiftcomp') + """ + evid = int(params.eval_id) + + # try: + print(' - eval {evid:d}: writing input files...'.format(evid=evid)) + for fn_t, fn_i in tmp_map.items(): + di.dprepro(template=fn_t, parameters=params, output=fn_i) + + # Solve homogenization + print(' - eval {evid:d}: running solver...'.format(evid=evid)) + sm = sga.solve(sg_xml, analysis, solver, scrnout) + + # Extract results needed + print(' - eval {evid:d}: extracting results...'.format(evid=evid)) + for k, v in res_map.items(): + tmp = sm.eff_props + v = v.strip().split('.') + for i, tag in enumerate(v): + if '#' in tag: + tag = int(tag.replace('#', '')) + if i > 0: + tag = tag - 1 + tmp = tmp[tag] + results[k].function = tmp + + # except: + # e = sys.exc_info()[0] + # print('Error:', e) + # results.fail() + # with open('results.out', 'w') as fout: + # results.write(stream=fout) + # # printInfo(evid, '>>> FAIL <<<') + # return + + # else: + # results.write() + # # printInfo(evid, '>>> FINISH <<<') + diff --git a/sg.py b/sg.py new file mode 100644 index 0000000..f0780ee --- /dev/null +++ b/sg.py @@ -0,0 +1,241 @@ +import numpy as np + + +# class Material(object): +# """ Material class +# """ + +# def __init__(self): +# self.name = '' #: Name of the material +# self.density = 0. #: Density of the material +# self.type = None #: Type of the material (isotropic, orthotropic, anisotropic) + +# # Elastic properties +# # isotropic: [[e, nu]] +# # orthotropic: [[e1, e2, e3, g12, g13, g23, nu12, nu13, nu23]] +# # aniostropic: 6x6 matrix +# self.property_elastic = [] + + +class MaterialSection(object): + """ A macroscopic structure model + + 3D constitutive + 2D plate/shell + 1D beam + """ + def __init__(self, smdim=3): + self.smdim = smdim + self.name = '' + + self.eff_props = { + 3: { + 'density': 0, + 'type': 0, # 0: isotropic, 1: orthotropic, 2: anisotropic + 'stiffness': [], + 'compliance': [], + 'constants': {} + }, + 2: { + 'density': 0, + 'stiffness': { + 'classical': [], + 'refined': [] + }, + 'compliance': { + 'classical': [], + 'refined': [] + }, + 'constants': { + 'inplane': {}, + 'flexural': {} + } + }, + 1: { + 'density': 0, + 'mass': { + 'origin': [], + 'masscenter': [] + }, + 'stiffness': { + 'classical': [], + 'refined': [] + }, + 'compliance': { + 'classical': [], + 'refined': [] + }, + 'shearcenter': [] + } + } + + self.strength = { + 'criterion': 0, + 'chara_len': 0, + 'constants': [] + } + + def __str__(self): + s = '\n' + s = s + 'Effective properties of the SG\n' + s = s + 'Structure model dimension: {0}\n'.format(self.smdim) + if self.smdim == 3: + pass + elif self.smdim == 2: + pass + elif self.smdim == 1: + pass + + return s + + def summary(self): + print('') + print('Effective properties of the SG') + print('Structure model dimension: {0}'.format(self.smdim)) + + ep = self.eff_props[self.smdim] + if self.smdim == 3: + pass + elif self.smdim == 2: + pass + elif self.smdim == 1: + stf = ep['stiffness'] + print('The Effective Stiffness Matrix') + for row in stf['classical']: + print(row) + if len(stf['refined']) > 0: + print('Generalized Timoshenko Stiffness') + for row in stf['refined']: + print(row) + + +class StructureGene(object): + """ An FE level structure gene model in the theory of MSG + """ + + def __init__(self, name, sgdim, smdim): + self.name = name #: Name of the SG + self.sgdim = sgdim #: Microscopic dimension of the SG (1, 2, 3) + #: Macroscopic structure (1: beam, 2: plate/shell, 3: solid/block) + self.smdim = smdim + + self.fn_gmsh_msh = self.name + '.msh' #: File name of the Gmsh mesh file + + # Analysis configurations + #: What analysis + #: 0 - homogenization (default) + #: 1 - dehomogenization/localization/recover + #: 2 - failure (SwiftComp only) + self.analysis = 0 + + #: Physics included in the analysis + #: 0 - elastic (default) + #: 1 - thermoelastic + #: 2 - conduction + #: 3 - piezoelectric/piezomagnetic + #: 4 - thermopiezoelectric/thermopiezomagnetic + #: 5 - piezoelectromagnetic + #: 6 - thermopiezoelectromagnetic + self.physics = 0 + + #: Macroscopic structural model + #: 0 - classical (default) + #: 1 - refined (e.g. generalized Timoshenko) + #: 2 - Vlasov model (beam only) + #: 3 - trapeze effect (beam only) + self.model = 0 + + self.trans_element = 0 #: Flag of transformation of elements + self.nonuniform_temperature = 0 #: Flag of uniform temperature + + self.initial_twist = 0.0 #: Initial twist (beam only) + self.initial_curvature = [0.0, 0.0] #: Initial curvature + self.oblique = [1.0, 0.0] #: Oblique (beam only) + + # Materials + # self.num_materials = 0 #: Number of materials + self.materials = {} #: Materials {mid: MaterialSection object, ...} + + # self.num_mocombos = 0 #: Number of material-orientation combinations + #: Material-orientation combinations + #: {cid: [mid, orientation], ...} + self.mocombos = {} + + # Discretization + # self.global_mesh_size = 0 #: Global mesh size + self.degen_element = 0 #: Flag of the type of elements (SC) + # self.num_nodes = 0 #: Number of nodes + # self.num_elements = 0 #: Number of elements + self.num_slavenodes = 0 #: Number of slave nodes + + #: Nodal coordinates + #: 3D SG: {nid: [y1, y2, y3], ...} + #: 2D SG: {nid: [y2, y3], ...} + #: 1D SG: {nid: [y3], ...} + self.nodes = {} + #: Elemental connectivities {eid: [nid1, nid2, ...], ...}, no zeros + self.elements = {} + self.elementids = [] + self.elementids1d = [] + self.elementids2d = [] + self.elementids3d = [] + + #: Assignment of property to each element + self.elem_prop = {} #: {eid: mid/cid, ...} + self.prop_elem = {} #: {mid/cid: [eid, ...], ...} + + #: Element local orientations + self.elem_orient = {} #: {eid: [[a1, a2, a3], [b1, b2, b3], [c1, c2, c3]], ...} + + self.omega = 1 + + # Dehomogenization + self.global_displacements = [] #: [u1, u2, u3] + self.global_rotations = [] #: [[C11, C12, C13], [C21, C22, C23], [C31, C32, C33]] + self.global_loads_type = 0 #: 0 - generalized stresses, 1 - generalized strains + + #: Global loads + #: 3D structures + #: generalized stresses = [s11, s22, s33, s23, s13, s12] + #: generalized strains = [e11, e22, e33, e23, e13, e12] + #: Kirchhoff-Love plate/shell model + #: generalized stresses = [N11, N22, N12, M11, M22, M12] + #: generalized strains = [e11, e22, 2e12, k11, k22, 2k12] + #: Reissner-Mindlin plate/shell model + #: generalized stresses = [N11, N22, N12, M11, M22, M12, N13, N23] + #: generalized strains = [e11, e22, 2e12, k11, k22, 2k12, g13, g23] + #: Euler-Bernoulli beam model + #: generalized stresses = [F1, M1, M2, M3] + #: generalized strains = [e11, k11, k12, k13] + #: Timoshenko beam model + #: generalized stresses = [F1, F2, F3, M1, M2, M3] + #: generalized strains = [e11, g12, g13, k11, k12, k13] + self.global_loads = [] + + self.global_loads_dist = [] #: Distributed loads for Timoshenko beam model (VABS only) + + def summary(self): + print('') + print('SUMMARY OF THE SG') + print('') + print('Structure gene dimension: {0}'.format(self.sgdim)) + print('Structure model dimension: {0}'.format(self.smdim)) + print('') + print('Number of nodes: {0}'.format(len(self.nodes))) + print('Number of elements: {0}'.format(len(self.elements))) + print('') + print('Number of materials: {0}'.format(len(self.materials))) + print('Number of material-orientation combinations: {0}'.format(len(self.mocombos))) + print('') + + def findMaterialByName(self, name): + for i, m in self.materials.items(): + if m.name == name: + return i + return 0 + + def findComboByMaterialOrientation(self, name, angle): + for i, mo in self.mocombos.items(): + if (self.materials[mo[0]].name == name) and (mo[1] == angle): + return i + return 0 diff --git a/sg.pyc b/sg.pyc new file mode 100644 index 0000000000000000000000000000000000000000..7a81f685d1df8d02cc3d73c5a9155003ea596d0d GIT binary patch literal 4843 zcmcIo-EJGl6`mz2N=r%Bj~yv;5@(zCLMUq2mI1eL)71YqT37~F4ve@=x2xrlT6=%! zouMMN)|K7&=nM2Ee9?F46ZEP;fj&U{eP_6&q_{yY6v}({%sDe>X3qb*`Jcs3v3cna zi5mZ!c;CWfJ_qRGpQkdVo_ng-P`aV|K$_~UQg4*{>UmT3S~YcMLTz!;&So0|(*NSO zZI^&Pj;*fJINQ;QP0M1kVFSE7c+AfLwA%`>K_S@lt*5w+?D-s$TioXPj0$%)XSnJa zKyW+Acc2wygFMB%cbOQXY-wTqmFgCR#o9=6&ZQToo8IJXHv!}C6 zNZ0NKyE@JVPpON9=i|fO?&2|Xo>H$p^;)S&)I-Emn4He=z&Pp+zFO+W0Q?#EYO5P5 zMzi3}0mVh!%(S66%XhOnU}M3ziCo@LlrZlEEnp^vEn+y?cD(MK6Sf5Q^81ASza4mx z6)aj+H+U!P2gT!<=Rw2kf)k5f$puTn$R9u>+ZKR_Q_9W)^mJiTduUrE91gYhfk5;9 zexZ%A9gLDPA7*J>B#tSTi3@8)hUSe{6j5XB@go08p zJoUy?LX`o0|7Al}zf{Ff#%y>sDZf?4YE5YirK<{wrkzcC&bHy+9rG#R{eGXBVqXVT zF>IJ>6O{eH?&)CXp_{XFGSUNBh`v`AH-kT4J?e-ziED_Ronh!C{j$nq>rBI&hIC`g zz0Fyhir0!F_(l|MQ~m`2<+Z(zx8^OuHshxuTBp-YwElMtp!J835PSe9mHqG$_1+jO zbAtvsuKHiA;x9GDss4#7OigjBe<#+QRU1bL?m296u&!PM_0`6`ujkcKQ@w72`szSE zv%taWs14W|Dd=TORqs@9fr~EItP_j=rk=feN89$8B^Gg&?NIc43Fvp5)^N&qTlSP; zLkTu6dkW{6O5dq-PJrT5r3WF65i)N`vN8d%hq};AAb--m;AxtdNDjrzGMFrt1Q4!s zAtO0F4*+#2RVN4MJMVc~|8(KZG-xQJ%-u;qS4b`3yjLGx%F+EZ44tcZU?D zgEe_eB_8Jt=B_A`^cY1bkE4u3=cK(diK1`~BEt;=Ds*}S!OkRog^;Tl2Zp9-x7*&r zb$`*)NhW5eSrQQBN5{E?$7}(dY7yR|8;~+Z1ZCAV3Mh%LQD{#ZgBlwYlN!0?@&I%1wxJfUPH8bI+Z zTg8Qmbf!@utfi9|Wr6z3QV?04qr$}O8D%mnD_Ygaegm17abHuCQMK>pF=oo$EPXX{ z%>W9D^tY}V$V;?1yQL9vJq-_=?0t@VkMhW5@xCrt+vTVZV?>M(={lwHAnh5uFy&!* zt#^vBb((O!C%p_MNSh)f7E?KjC{lQ6nWz)$xy9U% zX=C>BM$|K&ot1IVMAGx|u*a$2o(9p@Wmi^Tv32-kV$}CE@0?A_p~2avnlAucr{p)5 zXZ?;p>o?dp`AC~fEq}&u`AxqqH1<;-MtmCb_WV`fRcfXrRsqQbe+6=^QcKM~U9J(i z<7TPBW(fsdzznG@ugDfHG90gcj65?X%aYk+y!MngcZ@tcCCd`lU^e078cYT6eEInC z?eMGM$%EkOqxFM4mb?@s>Rkxwgskyd{ou4NS^fKG8 zomb21860XyTHI|(tecc4x0J}n*G+CD$WFzaw7sdo{EQV~nZNPj)5FAzTyZAgLLfKvFu z0qr!8X@~nb?YJVt{lo{tTg1dAcxL&p1#yXK35ttBro-FdTDAardDhCJw@_JkSvc+- z);fD4*;$1B!GJF6qJ*r@5lTz7y4)K@{K)cay z_}V|`w?%iZqK$HVt5R9T(0!ZlzQ89PZL31A{vnPf{lMq0v;BgRE$zm_`L*u<0DW{7 AM*si- literal 0 HcmV?d00001 diff --git a/utils.py b/utils.py new file mode 100644 index 0000000..e815fef --- /dev/null +++ b/utils.py @@ -0,0 +1,281 @@ +import os +import math +import numpy as np +import pandas as pd +import xml.etree.ElementTree as et +import msgpi.sg as sgm + + +def tilde(v): + """ + ~v_{ij} = -e_{ijk} * v_k + + :param v: 3x1 vector + :type v: list or array-like + """ + vtilde = np.zeros((3, 3)) + vtilde[0, 1] = -v[2, 0] + vtilde[0, 2] = v[1, 0] + vtilde[1, 0] = v[2, 0] + vtilde[1, 2] = -v[0, 0] + vtilde[2, 0] = -v[1, 0] + vtilde[2, 1] = v[0, 0] + + return vtilde + + +def calcRotationTensorFromParameters(rp, method=''): + """ Calculate rotation tensor from rotation parameters + + :param rp: Rotation parameters (3x1) + :type rp: list-like + + :param method: Rotation parameter type + er - Eular-Rodrigues + wm - Wiener-Milenkovic + :type method: string + """ + C = np.zeros((3, 3)) + if method == 'er': + # Eular-Rodrigues + t = np.dot(rp.T, rp)[0, 0] + nmr = (1 - t / 4.0) * np.eye(3) + np.dot(rp, rp.T) / 2.0 - tilde(rp) + dnm = 1 + t / 4.0 + C = nmr / dnm + elif method == 'wm': + # Wiener-Milenkovic + t = np.dot(rp.T, rp)[0, 0] + c0 = 2 - t / 8.0 + nmr = (c0**2 - t) * np.eye(3) - 2.0 * c0 * tilde(rp) + 2.0 * np.dot(rp, rp.T) + dnm = (4 - c0)**2 + C = nmr / dnm + + return C + + +def listToString(flist, delimiter='', fmt=''): + sfmt = '{0:' + fmt + '}' + s = '' + for i, n in enumerate(flist): + if i > 0: + s = s + delimiter + ' ' + s += sfmt.format(n) + return s + + +def angleToCosine2D(angle): + angle = math.radians(angle) + c = [ + [math.cos(angle), -math.sin(angle)], + [math.sin(angle), math.cos(angle)] + ] + return c + + +def calcBasicRotation3D(angle, axis=1): + c = np.array([ + [1., 0., 0.], + [0., np.cos(angle), -np.sin(angle)], + [0., np.sin(angle), np.cos(angle)] + ]) + if axis == 2: + p = np.array([ + [0., 1., 0.], + [0., 0., 1.], + [1., 0., 0.] + ]) + c = np.dot(p.T, np.dot(c, p)) + elif axis == 3: + p = np.array([ + [0., 0., 1.], + [1., 0., 0.], + [0., 1., 0.] + ]) + c = np.dot(p.T, np.dot(c, p)) + + return c + + +def calcGeneralRotation3D(angles, order=[1, 2, 3]): + ''' Calculate the general rotation matrix + + :param angles: Three rotation angles in radians + :type angles: list + + :param order: The order of axis of the rotation operation + :type order: list + + :return: The general rotation matrix + :rtype: numpy.array + ''' + c = np.eye(3) + for angle, axis in zip(angles, order): + ci = calcBasicRotation3D(angle, axis) + c = np.dot(ci, c) + return c + + +def rotateVectorByAngle2D(v2d, angle): + c = angleToCosine2D(angle) + vp = [ + c[0][0] * v2d[0] + c[0][1] * v2d[1], + c[1][0] * v2d[0] + c[1][1] * v2d[1] + ] + return vp + + +def calcCab(a, b): + '''Calculate the direction cosine matrix between frame a and b + + :math:`C_{ij} = a_i\ \cdot\ b_j` + + :param a: list of three a basis (a_1, a_2, a_3) + :type a: list + + :param b: list of three b basis (b_1, b_2, b_3) + :type b: list + + :return: 3a3 matrix of the direction cosine + :rtype: numpy.array + + :Example: + + >>> a = [ + ... [1., 0., 0.], + ... [0., 1., 0.], + ... [0., 0., 1.] + ... ] + >>> b = [ + ... [], + ... [], + ... [] + ... ] + >>> utilities.calcCab(a, b) + ''' + a = np.array(a) + b = np.array(b) + Cab = np.zeros((3, 3)) + for i in range(3): + for j in range(3): + Cab[i][j] = np.dot(a[i], b[j]) + return Cab + + +def distance(p1, p2): + s2 = (p1[0] - p2[0])**2 + (p1[1] - p2[1])**2 + return math.sqrt(s2) + + +def ss(v1, v2, scale=None): + if scale is None: + scale = [1., ] * len(v1) + assert len(scale) == len(v1) + sumsqr = 0.0 + for i in range(len(v1)): + sumsqr = sumsqr + (v1[i] * scale[i] - v2[i] * scale[i])**2 + return sumsqr + + +def parseLayupCode(s_code): + l_code = [] + ilsb = s_code.find('[') + irsb = s_code.find(']') + # print ilsb, irsb + angles = s_code[ilsb + 1:irsb] + # print angles + ns = 0 + if irsb + 1 < len(s_code): + ns = s_code[irsb + 1:].strip('s') + if len(ns) == 0: + ns = '1' + ns = int(ns) + # print ns + + # ilrbs = [] + # irrbs = [] + islash_level1 = [] # index + is_in_brackets = False + for i, c in enumerate(angles): + if c == '(': + is_in_brackets = True + # ilrbs.append(i) + if c == ')': + is_in_brackets = False + # irrbs.append(i) + if (c == '/') and (not is_in_brackets): + islash_level1.append(i) + # print ilrbs + # print irrbs + # print islash_level1 + + islash_level1.append(len(angles)) + list_angles = [] + i = 0 + for j in islash_level1: + list_angles.append(angles[i:j]) + i = j + 1 + + # print list_angles + + for i, a in enumerate(list_angles): + temp = a.split(':') + temp[0] = temp[0].strip('(').strip(')') + # print temp[0] + repeat = 1 + if len(temp) > 1: + repeat = int(temp[-1]) + for k in range(repeat): + if '/' in temp[0]: + temp2 = temp[0].split('/') + for l in temp2: + l_code.append(l) + else: + l_code.append(temp[0]) + + l_code = list(map(float, l_code)) + # print l_code + + for s in range(ns): + temp = l_code[:] + temp.reverse() + # print l_code + # print temp + l_code = l_code + temp + + # for angle in l_code: + # print angle + + return l_code + + +def updateXMLElement(parent, tag, value): + """ Update/Create an element with given tag and value. + + :param parent: Parent element + :type parent: xml.etree.ElementTree.Element + + :param tag: Tag of the element + :type tag: string + + :param value: Value of the element + :type value: string + """ + _se = parent.find(tag) + if _se is None: + _se = et.SubElement(parent, tag) + _se.text = value + + +class Polynomial2DSP(object): + def __init__(self, order, coeffs): + self.order = order + self.coeffs = np.asarray(coeffs) + + def __call__(self, x): + basis = [] + for i in range(self.order+1): + for j in range(i+1): + basis.append(x[0]**(i-j) * x[1]**j) + basis = np.asarray(basis) + + return np.dot(self.coeffs, basis) diff --git a/utils.pyc b/utils.pyc new file mode 100644 index 0000000000000000000000000000000000000000..77e1dc0bdeaa7cefe5207037836d900d5d0f3901 GIT binary patch literal 8526 zcmcIp%~Kr574MnlYhn3@0J1EN964Uw0t-}SrzlC51yWLuMJXd}VhOUI-DzN8znC7- zVzWslQst6U@;Bs?@-fH!1*ys@2cMJwAe9_*NPfT9v&$~Rhg6w`>FIv``gQl~_tkGQ z{~jK!=CAy`q|#po?|XR6zfgqucT`!ar|r6`d;Mx}Kt0X2 zAcJaeNIlJ|&{Z!%ZCHgE*rQ&8+KANA>{Ty8a8!kTf-|PVeiaS~&bZV;ZBXq1b4Y*_ zD$EJyqzZ?nGNr;1sZ4X(;i#a^sBlbxvnm{y%A5)(l)j|4fip?S3#C4P`WnYP$xmV1 z(SPwU`(Hy8#)23U)^SNK+Un0wXH$h8nUfi ziaKnYBWZKgHpeXPI9goF1n1mDG3epn*k&8W zHqRh)m=Ngb<%n~~o`+0&7Fl{v$DTnxZyBd!rI#awRKjMpp7i4NrH*UHLfL!)lZSnMRJFF;37w>Wn%g_>;E_IZ~>`1@`0N3OTpAhQ*~KHY6hE zDhiXU_FWbKR=NOz3}#rqaH#`xO;M&3(V^*C4)38bLqez6kO6iK2lT1w87z>)9W~no zYL80C&NKl>t5;1oT~+{BOwRDnrQ>|P1xMb}h& z=BZgg+BDnCo0{Q@28laMI^8Q$>coX+jF461V3oHLluL)@Akki2OM)b-RlP)44di?~ zu2oKfq8>WIv3BW8TwnJ5*1%I7FAptxyF-77M+LHa)2VS4o$S`a^Qd^#JVB-UcKmeg z&3g~9>UiF-g>kfVsLgqJ&nv*&h;&uQ^N*skuI|^4qEg4IT63Euj;cG-A?hr0uj5mN z-T}Yx8@rcQ)IC9YUJ=Dlf!@7{md+~S8Rh5bv+7M)o?7Iy=2NQ@4hw5Z!p6F^^s)AZ zEX7=Ajq9R&cJYF_B4r}@)k0vD;NKWhM*F|I=LSi1%jEbDR1{`w< zg&G=la_*dS8@f9zwGn5Pb+^Rra$ju^q0xZmK}ctyvolo8jDkQnKBlN!m0R7qtu|c+ z!KSvrY4)i79u)3?G{8 zvi8zZLpyF)yKRAhhc5intJD2*mbEV1IB01D-dx|f}6u6S{AE?hnry@JQQfh}1!F+74Rz>QjAEC#`jxzphp8Dgcqw^avC5O?`c zXH=3Ypet~MSHu~$=)6kS&Skgve@u~j7oOqT4hsDK`C(dtW>MtCzMfh=ar6>V`Pb0v zzrlhh;ZvPNc>M39CK41hBI93Yi=tJ(0NI!XUD~z#fr;A2?bhl~05)8jf)?eRtIk_) zPHH)4+Pw}4p;C8GH!0#n^WN67bu`-IV^6(tmo^=C4abjjTHTr#(qrl-!j{pUld)z5OzCjv$*%Tx==y^4sTYS1#!;iogNw46Wp8~q)*V5B1$$A| zOVzc*ZMKHc)O69R#i5Q#N+VE4)vIlLJT$(!o_2^)tJ9hy*yx~|4(=ct>*O%5T6*8M z6@~@NF(ZX#4@SCvydcta>dW(#M(be*$$wh>FQdBuHo=rM5k|4=c3&crvKWyy5kLPe z4kepQ(FKUzGzmlOeq_(9b9@Dut*(tSCX&l8B+0 zgI3Q$u=tvSl<*CO77$)jGlZdX7LA}zG#Ydf#3y;OZY@HQ-g2mg{+lSQk!SbOQfL__ zYD$J3U4)*JVIoZWGr6hVKmLe5V`=pW9jPVHxo zWyvfv2L_!6;M9;Xim@R^;qZ*`O`g?Quc(vfT2~R2Bf9beFQ|rIG1Waf*5xl&TDRMv z5N&#!B}9>%Ua_$EMVE#Kt$QpLk_sPqMc6vy<%7cF4KFAxvAES4x~D_B$P~|#DYnQI zNv23L#cnd`kd#bIwQmKtT3bwQ#90nl3vrHy+l~i~pi(dEWlQ(&-MiS)9q*ZRE-Wm_ zt>M|?!VT}HJX<%qdkD3?=k?g$t@aRG<0QpbCpk;tmG(Jgc;AhZNSg(!|4LBI2Vk7# zzV(z;A)8{L#B><5AE0UP43zE|wo6-c#-=#PK+2+6JdpM;w>tJNz=SbIO84+fRBhCG znOF5uXA*z!Si{9lwEf`+Xrva3(GLvNEKJWNK~{(O63`3+m{e!s{1{x4NNS6OMYKiY z2mPz4*|7T{PQ+%3;jS;**OI)<(gpB<-Xw%i45}sl2h@sUBgheqo`9M@uh2_O+v;q;N^ z`_Yp4RRbBZmnd^sF$b{-vtK@nDL6)!0$Y{2&;l6Z%}D!~xaJsuN^&ZHDQu3Q-JD zL5RQsq5h9_dkT@Aq*NBWanb90j2Mf+d`wQv2N(tbu#g#GA3Sq6(0Y>R4So&0SL&IRTXGYX zn9FBEMSzP^=B0dF%J-yPmb|+p_rw5hN7c~3%TW|`Cb{aOw^20_$itv)w9i;d(AVQb ztC~T*uB)LPREK{OjTRJ1bkLYp0-F`yF4uy@-s<%cl7l)=y^?6!NM={6Wm}~=M1Hp@ zq@vgq{db8>T}3)zn7yVPnB77dw91QOTOw14s)b@4l=gLE#9Y`b;1JNcoPv$o{rA!K zZ=n$CiLTfzS*+_CD8;l6F}*dI82~1d|~bYV;?;@nL-h`RzXeVm?E`3nW~k zJMQ+uAk8=vPSzRi?GwY4bF#QrX5C3=s&^8QK4ieb>w62{Bv+(;4*SG9_#@Dxx zPO&=D!fQ72`UGtC$y*Uxoj~=KdU9Ay7o3V-BeBksR}Q;5Q9G?r2!?f}2K^gAm+QG$k-*6H zU93i_rnk3^jo>KUMwhoOdJU}$VB7%67GQ3oq)=csC5_**HgxTZ-+lvBdPP??p2_7~ zlBZzK43uqk>!Dl{??{iRodv%uaT#W%r}5c{*d`=r*1;h(qgeKuYQVdoBn``E1O#r{&7Ek9fcUfhKQ!z z1+9^#J-sj=p{5WQ9Vd79N(-^bat=$8Dfv zd9u*&>f@td0qx<|`S&^5p8RD%MCEJ7-Y>{JNjF*1wVGZH10z`VFrimH*PscX#mpgS zP^!b!`A{{LaOd(dbiwQhF}!%hW~U@J{95Wx9xe|r7_Q%