From 84a6fb53999d205bced5757f7c26a6bcdefc80f1 Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Fri, 22 Jul 2022 10:17:30 -0700 Subject: [PATCH] Refine the tutorial --- examples/multi_material_tutorial/diagram2.png | Bin 0 -> 14028 bytes .../optimism_tutorial.ipynb | 360 ++++++++++-------- 2 files changed, 192 insertions(+), 168 deletions(-) create mode 100644 examples/multi_material_tutorial/diagram2.png diff --git a/examples/multi_material_tutorial/diagram2.png b/examples/multi_material_tutorial/diagram2.png new file mode 100644 index 0000000000000000000000000000000000000000..9cbafb49a5e9bd01412fee45e5714ae85f8353f6 GIT binary patch literal 14028 zcmeHuWmH^EwKp@QF z{l0tW&RTc=&diT#R*wg*f`qf6YNyfDrbmzt6 z@n~rT66BW`v7UN6bCD?W>7zC16@+1WgfMgHdJ`LQT+DJ`d;jM6 zV&Q5-E8axo%~OIkE%V&>(A=jbD0u zxB6x#BbT715+&L8Mb|mTi&(C)-74HKB`!lSauXezI-@>rl&%=%j?*0J%-RsMQY32Y zM1at?pWiLLka@x5V`QxjO-12;tDr(i+&*?OUv5Jr;0=~NO6&GL4#RKyGwtDVLYn7+ zc+dOklN_pOae3y-=TCtk3MwSeBXuEt<)iXn`#wW|F!}OOVo%sRZ1zVO0)n68d<)*b zv%@*1MR}Ti(Zclp2USdI1lK(VIfrYX9=hxdFb;tZ#wv~)8c3Xg4MsvkCPG35EM(vh z37HBB?H?NnNfnv;KW$xP_J5Z_K|+diK?41|j45z@{J?=P;Qj9-N*T(3miSnP`X4s3 z>SIy9NN+sgfQf*c_#z=;6F+{DkqU|^kdRQ&UG$Cpj5XB7?YumoHuhe&4$wdk#G@CI zWS}@;dN}ymKmtA7J$=Okr5OJyAr9D&r+FD6{}l0alVUX1(1yr+`8Ys?p?pw2MrkYv z1S08U?EsAOITh1nTAE#LF)xCdSJrz$+lY1C-$Ledg(B z6UgJ~%k-~G{-YiR2VXlM7lfaSmnY=0UK?94e?KWk#>a;K^YgEM`nfp%&z3xW|9xA) z4e~x-;pK<&@&0#i4uLNJ8@II4nJRCfI|M3wM;FIM2Z!i6~t|tG*l~444as4ls{_QHs`}m~(%ai-p zu>5lt7#C?QN#6e)V`;1*{}WFnBmpNC1zG(-4;%c39(KsP^|kTN3z1lQWRk$tcBs>;V39?@!m?Kc3eJ$ z2J8z9Y4~f2=vi1;_-+<82VEUCZ3bI?`SUZ_Ad`si#WDvc?J@wszJ9J>t%yU2djgFgN}`jO^>Tp zP>^`&!~O5&HvwL~Uvtd5-P3fZz6*tSN@r6ePwJ0FMxxMNFW~s4VMl)&d@kev{bM>^ z*tJd2=`hCO<>8VJAzF0gWxd6lr(b_cdGsF;*0zE&$!v#2oabTCN^xGxzP#N~psXLz z#AK0ly_{eMV;T&kTg*+JQQb#hMZUg08t1Y%mHyqT!x6Xt>Vw($$@}y3v_C~K{KW7$ zLPVCeHvY^XOtOppca^9Zc-pJY7GmE9{`{=;8obXr(1(~jyHTi2R5b}Q1QEr+V!-Gs zL`kUL)&w^{-<$+)#JnmL@v6b26-NkblWr5RedrS_rz1g5Wk6VhP<2ivtV9BSn%^EZ zwA5-83cG(44Z37^@;@jJSp50rjY$y4@OyF#jRO|Y`y}f$zU%Rd8k6ePXiQ=X)V$B> za|5A&zKvk>?zR5jd|o=&w)cYt)(K4vCylsZ{MhS%{4IvuvP>nDb)n5}QgqOnr@4MNk{S*&rm2I+Ug5F2QgerkYy>i`G#I6B zgj{HN^8Sp;6ld&j{h3`@Xog-|Z_sd#N5CzR6ns?oGSGrsw}i-OE-^fpFg}P@y$Vym z*R<1hp-%6m^H@CL`s<1cMpZR6?9ZAE8?QY+f7>tgs&RZ*aMqk+&|pZoZ-?mG43lVy zvUx)a#*v2?i~9#Q#|f-(piIZZ^oh9jGu*%QX0P=8_}rz}VAxc_8sX1H0>0Tv;%!$2 zZw*}UU>_~l5!)!o6L1o_WD>eJ`F@p$RDg)H35c9WV{tFH|F!J^8Ar9!+{Tqu?bW#8kw>KM@fG@5icuXlS3I=!b654#in z5_H^jc|OCd=jg~v%~`<$0&`WN^^B$3mOH+s%;5SYDVF{?q>KbP*0y=ZgJDv*wBlK1 zc(yI&v*W-xOlY8}RIs4##3XSoa_iEh+Sabo!#Bg$gR+c?Wh&>P*K~JNv?4C|En2}UIAH_&O7U6@*dcj}R#95xf z#a$C5PY7nzh7L!=%MYyaj&wjI5v|~lRAM6&Pp@Os$kr0A-s6#1uW_JY%9fd9l6f?h z$-TWbsS>R4v7IY-Pn8Qf8x^Pg!Mce!OMI7$qwbI6qA}4HvI5(A3DM!kT!?JP!^gVI zKByRO?32742nHht&xaD5)s?~Mi)C8~RB3aAQ)pgyQ4!5A+vj$Q7Mv5E$sS_Lg&dPz zjszYf7fEbnpD)y!>iS;5VlY|v#Gjor|M=9&bzQIW3>Sm4%YBU~IrhtKV8OxHuZpv% z3bC(7jUvomIS-G;HnyFAQwR$`$BJJ}M!``b>zB9BWEa}}{!@lB8`M7Y+1;h4=#na{ zgz9E#!WYLUIiDFA(# zAmI9RTRo>(BH-M*W8jHnE}qui-Cf8?fjO`asqK4K-&!9gY&`U({WLRE$X$h(6Ynb? zI=5+PaMUb0)oW)?xYhbM#LfQhH>0R0p{$4ya8Tj)&qRn(i)lWKkoPmY7}UwE_rM|h z_(S(;J+^H*Z6LeU{pFJ1-W1^kY}tO=Y>reODbTRNFhxXht6_n{cOhBhjV^!naF^ZB z{@<2IOE2-Uqnux*ViDtLq3ieBjK*HhSp8jW_9PrVMSpURwDx{AM!Faarw;EcG3^|q zxb0n|fKs~H&e!}0E1j;|yz`{N>YWO4)@;HU^%iykJ{{fJMO(`)2r0E~jPD}I9pR&y|wcK00Et-p}7OdN4J48Lq zg!{{{!+*hrcbj#K#nke%v+m=eHLOuSDQZIcF?}zL1!^(biOc=BqSW64% zuFoldEt?h#$wZsA$lqt+ZMXGC6LWvSVpF^O_=eHMbpJfe?U#BLUJCJ$hcc}Y8F2Cqqm%C2)(dCQS?$@=)hC1qfp%n1M0<{o3!!$T!Kdcc*f$a3M6=`lf)9{rb3b+Y?d-@($QfW`-@5$n8qzYk@RzMaUZADa|mX|8zL)2FIu! zwehND9n73%4avMyu$2wEziN42jY}=eH+@QK=${B{XG|lAjL?TyVztV@(3`5JU(Tewn+0Ie0uFc{l#}q z>XrAd(hI%8uuv+=Kt44nTtJ81Uo!_sUhNnRGEcWYx8mM^S2+0qZcu9w@wh{sG2pNyozn+UK$8)8J8+LhG)EhQY zZ7NY!d2IBFSv*~Hph3SP%E0mE&$cgTKk8g1h?f|)55;~kA>d>zx_S`Cgg~&TqDAAe z*~UQmuy)|d>B-6>`=4xBtSBOD&Xp)Zeu7X{w%R*7@|jO}pYM-MOI^z{rV(C18LG$9 z#fs8!hzNY!fmgbHSe}J4P4_?%w4=!`^Q|K3x+pt)JQ|yY=NtSg%i`3&P0)iA_H;K`H2)aQ3h_qVHq z`~9(eGhUkm#qTZ(g*_?&0KCYQDk+CcCAei8DQHNpDlj|U+EibdpBZ6)*eC7y+Qx+V zu`>{q=DfoN=`=Lv28%bs_cwFTYk=%v>GtaQa3}&btFUn9*!o)Q6!HkEN91gk+qI-` zqQ(83y~rxtb}O4WUWPAmEsTuRikQm{2}PDm$qCek0qAp=9!}0#cS^s*B*rz($0qSc z^TUgQwVrHOLpQs>#eHj=U{q{lLtCeWtm+Lt&`u=Ggq9nsnDvgtK60g{wKdHP++!U^ z?pC&Pts{K3#>F*lO%u@z4p^G|ys++2y7bVHugHdt`z+f$If;%a3TCfE7$+`kYm!0^ zA)d-ePbcKYUF8?c(PzQhK086tcLvXn>C*DEQV6wQfV)-sbTG~e-c~3kfeS+}E-s!~ ze2n#k^Wfhifz6;+@Eno|crNklP2zgAyUwRLq6QW(YC!YK=x=onphU zSBPnsj+sUGxy74YQsFB$BxEBnYxEFaa_~ui?#8nN+M~4&r_+AkSJa$4ZY}*v#?d*? z@@J4Gs>&^*$c!AodpHqBNAj~LPU(Uu3Z}P-bk<|0ATT` zgb?&qLjcg@9bjLGl7KM}sX`k}+4n})vE6Kxo~jT@aiL%;4v}(IQ1ImNTD?z^0>D_N z+6W;yAv(tCd(C8|UoElEcAwj7a~a?!cl$~|l<*0ONWUgGiZ>et7#YJ3-BLLsDjikx z=UI$cIO3xQJ=F_0xM)#18PvQxEz|lp#8=?n5b7caA10H%txcap!Nl~A6 z0>T?YK$gtkMZ~rjxwcVHhSg)tvv(htAqOO4T?B}1YgBZcvItZRB57JAF2-dkQ3A}( z1LHSOdA=Y74eD8d*Mzx&3`I zXQ%+gG5nM82y@^Ji`2#JlRGTz#fZQ+0VjA=g3;#U2Poo48OX>+%oO1S8QV?vH>^90 z)oRapa5$MPC&tLDe0+Sa?W*8Wiig=;Mk5WrP$gOX>DX{L5GtRX0KDC;xhM8{o>k1x z67m(^ETB5o*W1-s{)&pPMsUtySEwdi(dOMsqDB6+>dSKLZjp~ zgfa{yBTakhbG_4y?FBn%Vp43>D3}a+;r137nS4cUn&_2r+>FQKxumpLY{V#-EP1rr z+GzBFTgZQ3Au$(a%;2%$s*!=rOQW(eaP$>Jx={D9*tXx!$e#fDpzM$@(8)AgRcH$i zVUw~V!71}_BYJA0pym=MN`ev)2h{}GO$QjaNPyt{GHp481eqELv%+Q0Y5)J+|Ji47 z36qY14WV|yABMX(0Fi@hSU04RN%t4B$NJ{TJDJ`T5h&WK~0N4_^xU?yDC_1VH0F zKD|oNES1&(fXeRivi=X($?#MLCPgXq);=;J9-b6v->IA-F+*K%z=@_+=oP8QT6+ut zBa~OIo=X{-tc)i9rg;?$a4&|F-Y>}XSd7Qp0;Di>egpB`1t3QBTP&qmsrKy=+N5!K zw5px1D?7*xXBI8BsaVnI;TQ!)5iWwvYFex~RkC4}e1V*X3jV~BbijM_Q%Q%?`%8@m zHLH>0xmHaZVQMdmRkBd+CxFMGB*tNu02r)Cj*d#ibCx+tkSsR)MWUe$pV? ze6if1UcVf`VvlSN1s=SSl*@3wOg8G;bcf}Mp=GTEo$le`;c}Okl`xGBCg_vAAhZsg zlbbu~F`@@TB^k8$fsoPRQCxo$l%$+HuwfXhv+NhY`nJFXmR+egs0%pR2;tG%+|sxj zOGZMVzo4X~ROe3n-g8h=38bhxkL=Q&idFQCBnLQFY2;mjP*O?AxymE2h_8*KyN`R0+-YEm)mlKo0}U!&vnt)IX06c09O~ei1a%tQO#g7n)e2l zt(d}ooh9#7_7|IsyUfh1KfAf_T=m0oHR>F>a#dGNwm_3j`kk({ks-bf;p4Ks>Gq{p zsI-!SoN>Z_FJ4%$S(XM*J7NnlaR2M`+5^B&b&NLxu{X=JDk>dK-(0^!JdtEO-jhxh zLn={SOYP-dU&yWB%}YTbMJtdF4QgctkT}lKxwg#mGi+qZ5WfAh@Sxvc0x&0-YGV-# zu3-MFBVf`RS0@PjzTA?N=dY6tx#n|9aN{nsedh$wr8MXj7GdM0%VWk*+Z*gDONwij zK`QN8mq=gt^s9HCAe*)Z()9~|B`pwcvd7i~a?GA3VXqDKA!1z$e4JVu&t0r7*AI=; z#DVBK_znc&HR#!4Cma2Ob`v;3C;h*#rAI1Bh4~5H7pqh*r4J`yIgBtsRDP#h zvYuN2F1)GB$jlri6Kq6_a9%kUdP!b-j~bd8O4$kbQg=$!oA!8traQ~H*H&=3jX}TA z{?J@GICw6(O6h1}%^HkKC>s3By>~>CWIn;O-W4M7Rsx`X-j&Lf76`db)Vsfmilx@a z#5*4vZKe?AJ{mRogB6xADI$*tbm{lcFGVz;eTcqVCrbP5fn{X*U@)lhA_SR_u>Y(k zec69(clqg7sGoS6JO}a<-b_rLJ}wWxMUu*Dx{q1P#Lu58w+gB!2YmlHyr_y|V;g7S z2S7xkth6VHs1Ce>X#Fv}7}GQ3c6-dZKMr5i=oEbSXUjdm^qvjOBXoIKBXKDOq6$IO<_-0IcFyU*^gRx$SvcPCTzfE6byDV+UeL7Dz#p%<0x!{^P|u(q1K z<|bW*Nt6lZpY{%nUmeCA%GhHHe$7+~y~$Z@CXGA%dDe9%Ez5ko-t)xktzs~&8EVlx zuOae9!NkO*c_*Q_4DZ@zRjSA=bUKD*82Mv3{2QBb@Fu_yXT*f)(<|W{cg~pkCMQ@M zS%vz`l}AJ%fkvS*@^N?zADYa1q>-!vMr*lFk8MCx!un|@wNyR^E;lXoD8W%4w^-aa zi>lLA{Y!IyJhb<_Z-I;#atggniLmGT)YT~TH^3D0Z!4ux8Px2G%Aw$$>ubUmO1O#a zFlJ-{Vt4nJ$7-O6{(CFPZX;wf^!M3>>lQi_csWFZoOt-V0D=iei&)UOq}^8Few6<2 z5Va})pxE&&iXTUvx`5Gij5tLP%%v{hJ>w-##qAZ z{DiXw^BRMKScN>}2^(B{nTVE>>L&`eU-$PH4U8;E8$F-Vk|x*EtrDb3WApIk4EVeo zA~T&fZmHIRkd;8TS+m>rsavkMzjFpcfF`Fi`7HhCLx@vhx^O-&bxXFK&gS&HX~t(v zow6s{sjx+O>v>@XKMC>)LmBA8veY>xBpRS9l#A;$&gcN_baIvOVM-a`972}!Kau?Z zLjP|)g1`*x@$qpIFen7nAzahkt`$6TH{gtGZj zx=DnMfYEoN{~|F~tpb2bt7YJp(=O4C#iSjNx2fwGG`fxzjqEM=T&@bBI2NYTnGSz+ zZ5n1AQfRn>7QZw~O`&z-4^}FtE&_tlwRwXfkL%~xJO(Y6#SOq>NGWXva0TiU$1Q^; z)9i;|8ar~T!TKQNRpz{G8v3wn23$JHyr5s-RUR`SI;o&}0HJBpKfmMeVv*{P{!J_K zw2IoZI|4{TC$`1e8sJ85pA;S$m5z+#t5vI)R(Ct(>mMdksP&#c75;i~u+;7_ljrxu zaTda0+#hsR8~W#$x&SkQ3M;20Bhr!pX@v1(wnZ)MZtc_Qyb!7KD&WMDlbhTA`AiOv z3lJ;(u2zgX?B$t&H(d=NAj&e?v`XX!xLoX9S|k31pk+PjpX06u9rThn`0FQ%c z!uD+jdZ`V_P1_qsD-B|n05P-DVpR%2z4yNSwMkRhzE6s>?Qt^cB$Hw52i9^ZY0dyI zg)@G(J9(2FyaS^c8)sW!Te`R<&U`!|WsBBIWAQU6Z-$Afv z&v%?s-yvV{Ca`%kOv3~y6c-m4^X$7ZRARfu**cJQqa-696V`fjB_;G>POwpvsgU%H z6a%LI7^gw~d^(fz^#F}RhK(k$%6F9=(*2IHXIX4jAkFrBR76u?>jnyf+WNylD(VxA zPa_I@q;>>kX=mTQtpMO^M%%erzA|k(O{B_=@3;yFB&IKY?LBovw zLfJ~^3Gg`yXqD*j*KVuOiJku68sy@~hQ7)_dRbpJL0`z^n6SrTAY(?DZr%X`P> zrVFTgn6dy9qjDUJPp=6)qFuE8lfj3QUq3!8i`s2J+a8Iz+@*iW<2#7#P~}tcVxc6f zWw>dTW{t(qOvFl>MYoY}?_ecSX>w8=vHKV@CJ*bbHi&2uTMZ z-rg?Qh$Zhu+Y7L9adX>d46{E*V;#p>7sIoBqR2=GC!u%nwIeCTA#ad56NM|D?a=U1 zZf2%D07!S-#NvhKT@|x}%bmEUdw+k^JanmZtW~2 zBVY51r-*^kJdjy3+P(*udhrmCmYt;o{XXy`qBro)o$AvwaIt8y(Zn7-2XFsN@E9tx z5hg2hWOO7|{0%u-u8m@}fp#+00+fzPrm+wZ-U)?&oUFf2|9kgRxmGoG<8qP~ZMjPo z>ua|97o!IjkMrmw}SRK$A=b z=$B=7lP?zbrZaE40?%wqO1-1`;@bBOxvNaSZIJ zGW0hw;ao&Gob*+TK*B@IOiq^6=8L*F(gQfHJLVnA-0Ttm`fNhnv^RoXcpWx*eYShZ zFW=0aiLaC1##nIm1g`#G5kC>O^Y|Jhyp6l!75&bUMSIaL4N%Q<<)^N1ftc*_zDNu# zJqcmlL*bBZF53&b)M0h80V-gv$@ddKvqic}Lx6Um!!RrYoJIJ)F$7#hpmq0A4lN2J ztXy--K;VE=2Yd1b?MwzX5rM95GIP1hl_?>d1YnKw1RJN|eJ zD1{2&!K3(K;=|x&DIp;u;U;`!1Tt3!s90hSw7TVH3klD%?=^8_!=Xhn?jMWOb92c- za#Jj>4EyE~@_q9(b4adgGsOJpjN*}*`owi?nrrdYr1i%CXxYXm90tK`O!x?mNg2Av z!zfC_esoJSuO~xi(AOe@Y+(&6ByKnDJ#HBQ7|$C6(nrc@i1H-KM)^LAL+_fF*E%j> zG}@fcIo!O0>{DHaN_dxHTL!~c0V&}kSmDia^PB58)aMQyTpM~XH645&5m40Gk1=94 z&DPF&Rj2`WQ*GTQFo)&%kdS;5-Ucjbh!9#1toQnJ6&{pBZIzWPl~O3U@}ste$nq|- z%`NHVOuOLVZY4f_IXJq^Ktg=!N5 zga|4?uUU$B4NDA111vP0D^3(mCP38z#3wDnRzTnB2rsu(9C87;C{3pN8x9mQz#`*Y zZUBp>0IZ=(BSf^WFsn&!`o6 z>LIJ+?{Pz8W1}@wM@Pq6X-U}ADYXLYD$&6AIb?COLTZXbxS(iuV2=&!9gLg61rPy; zcVzxv+kNWsWh)l*jm5&oJ_I(4)QI~Zjop_@aSlgXl$F>k!oac5#ke0f%|7XNx*8ag zo(%l)#H=X&d-+l@@N{d-Us6gOBZ;lIpo2(hxFV3< zFL!&}9Rbq1Se$9e-mLGDCN5D0EOc|<-A-=_~T<+RM=I6l}?Lgzvye~Jvmv`*2VV=dgn|K zPOBs`9%mqfUWXWC0Yv!wzd*1)@56DkArjeWGK?Ijwt2TH1cFi3h9Q9Zuyg)Aa#`yQ zfC@#kPu8QIN5uUR_kGa%8rzn8@qHRKcbqmU!}#A{U4#E*bS7OJL0ohF>sglU9|w=E z9iuNGlMFHZ(BFt0lGEaTP7xcvku|}68>;LM&w~Wc2qjK(WO!S0GLGd5u;`OSm_P2> z@*^Fl9vffhA;bKp7D;MXRES(+c&nD;)5>Jods9vIw5rO!fE%)PCZR_2vl<+yFQBEZ zhkRGz?Wys)vM$ML2VBEA9`yU?yEhW^(#M8%7aqLvWIWR$ihM4CPn8OWp|DNTJS&@U z0bF0zTop!JtcD-cC>6Jg<9@w?No&aS^8_oPEHysV?i((=ycjk4*4ajN)fx7c0j4x+ ztxv$Lco9IQ@Y#ltx=bVrY|S06E+30y<3lSQl81$b#c5v_-^vhd4j`reX&Z|UPFtE#mq%m*XX;zP{BK33?n75 zM{=tIY!#72XI|ibmnE4dbCZrw?fC+Tc$j9*k76bMxbwLhKdP?M+D+DRf}nu4PQKIM zbMl~5k-_aC9@BFUi@VVQzna!DZfq*WueztbsEZ0OXscW}q7UmR1^?`(eoGuST1;X3 zo60KnU{r6`MQ8SUUbY->UOw@r3|p>nKi|z7`v;N?8$wA&Nxl(+R5M@o)G-|W_+eR3 z%VrisJmB!X)3s5PxlYiaC3eH-w zh`r62kiJ)uE_WKsDaYXXPi0PC-GTB$xtX>%uhodDSXbM!cw90J6<=i>=R1{@ zW6+1Bg`GyxYaf4TdKQd>elZXl?DFfo*t%gF&Ijfy;%r%%XlX}yO5Y>>3Q)fl-y1nE zE9tLZt}FHQMNxnD&EYYdXwg>B7Z|Ti@*k~FI?%(=*bN0`=ti!kA3yHkdsp1)w(96c z2u>^f9Nz^0hOVxDhEhTy$GIrBDM$jCPupNa<$%o>u9xw{#_n_l;`06 zuUcjk#iMGM{y5koAvaCWr%pHJM>cmDw+;n^(YlhY3fr80I+{J?$=j%^HdQFIA|;T8 zfK}me?eAN$3Zxldo__XFZD;rXrS__dF^{&VRnAKp>PIjYZk5TZ@t%Pa%=o+V$1Gy6 zcPEjn79dg8q0`tm*wc_{Xzw@1+mq)H4H@c4^~?uV`B5g*SLAky4=McP0f&n*fchQ$ zVz2!dEj(iQ#pXtTFaW163NzQzv58!_h9YErR7}dhjt+U>Ub@`g9=9AWX~|SSszHCx zQlEb-n#Ib}x>sQs{yr76Iyl2beml!qR|b9CMIpB`4!?I(hR%L1BQs?oe~aZ0jsg+~ zoFm)KfzZR^fQ{{d%!)Q>?EcHu>6q{XOh*>-#+Wn!v~2ysF86Il%+CyA=$v|>!)Zmi zrnO%XOc$g^GvK$8?MT1gnkt&1+I-K;95ps-_Bsy`hh48*2`_J}$Oe4+42 z=V^x@&@{9t+Cf00pI`6sy?H81JdkaEaWIEZMC{wVzP>Kt3w0CQzMIloR^;Dq=%v8p zc^(~{WBqh5(v;$T8KGK0xHTdVco4I$Z!~rR>k22Jy0QhFJZ=^p3P)aUG(Hnly=+%?Rjrsq}9a0ZA^@ zLX}Ok2>cVpyMrig8A!yWsqH}cA~UP)MT&r$CJ}fFd2UC;+qKOi2B(6TlUoKxScYC1 z2=f(5KFfGPFZ14&iP-{3t2uaidD*}4i836X76KbJq;B>yC7n}gB_ki21VG=(Y4A;{2}hHZ20aj4y(}CNay4-_4`fh66le zrRz#P+F%ME=GC9ghmzG1m}6B`HsN9;0A~G0aYYam=dk`|77mmeRH?7CLmf?XG7mRe z?l$;&osHo-#~;0f=gx0Rf`4LCw6E0^N(Bo~@}o-Sy$_j48V&o^QHnN8Gb?Q2yp_SC zhBDDTe3mF zXk^Xv$ls@*JavTE2eR5&2NLh z3jXa0h!^6|i%G*~({Xk}@G(N7)^Xqi6Bw-9vs(uz=GF&jGy~6&+0N4F!e-0225yhI zq?w9^yg(z#4|f-H7om-X!wh&~wb=W9v3oSo8*Ay98-^p?-<@kwV76yW#SAtG&#Q5A ztJ6@iGPven9D zwKE;c(73Dvj{gIW5+cGg)M>~RhKIJPBRLrj2o>LE;zvv|@kabeACB#0`!1E0(ne9# z&Gx&=z}>+Tn-D|I>Ra4jmf|YXB>E3MTQG_qwws}I42(lewTcy$pvLN9n?G5i#?ohr zSVebSdedfL-kWU=Cs77Et8c^P;J_TvnOa`%qw3+)rK2BbZ!GO9*#l;e5+qQO8Ht}{xWTZMIdY2v# zhhDmYmZXl86}Efx{oBGBh^Vn|<+k2n-2dYk#Yat=W|)i(ciLBpr&1pU8r(y3>SQc<4k6%Eu*;M> zmz7GJqxR|>_SSJyae5keDzz>m@f)?3xq$XM0@%WLE;M3OZ**VoM3GlM&Sz$5MB>LGFLYCw*HpLG=y~fA29} zIg5mvuc+$3p>(tmdObP$+?6^ru$Q5B-OpE>_#wzgDd-dzMMajcjMjx{ttLcf<6EM literal 0 HcmV?d00001 diff --git a/examples/multi_material_tutorial/optimism_tutorial.ipynb b/examples/multi_material_tutorial/optimism_tutorial.ipynb index 54af4951..6f375ab2 100644 --- a/examples/multi_material_tutorial/optimism_tutorial.ipynb +++ b/examples/multi_material_tutorial/optimism_tutorial.ipynb @@ -2,16 +2,42 @@ "cells": [ { "cell_type": "markdown", - "id": "b5d338a3", + "id": "6356bf65", "metadata": {}, "source": [ - "# Tutorial: a multi-material simulation with optimism" + "# Tutorial: a multi-material problem" + ] + }, + { + "cell_type": "markdown", + "id": "3a2ebec7", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "In this problem, we'll subject a rectangular slab to plane strain uniaxial tension. As shown in the figure, the slab is made of two equal-size layers made of different materials. One material will use a hyperelastic material model, while the other will use a $J_2$ plasticity model. The slab will be loaded by applying an extensional displacement $d$ to the top edge." + ] + }, + { + "cell_type": "markdown", + "id": "b60dcbd2", + "metadata": {}, + "source": [ + "![diagram](./diagram2.png)" + ] + }, + { + "cell_type": "markdown", + "id": "045899d7", + "metadata": {}, + "source": [ + "## Imports\n" ] }, { "cell_type": "code", "execution_count": 1, - "id": "24fd5985", + "id": "3f49e265", "metadata": {}, "outputs": [ { @@ -41,25 +67,7 @@ }, { "cell_type": "markdown", - "id": "9c7de108", - "metadata": {}, - "source": [ - "## Overview of the problem\n", - "\n", - "In this problem, we'll subject a rectangular slab to plane strain uniaxial tension. As shown in the figure, the slab is made of two equal-size layers made of different materials. One material will use a hyperelastic material model, while the other will use a $J_2$ plasticity model." - ] - }, - { - "cell_type": "markdown", - "id": "37c03013", - "metadata": {}, - "source": [ - "*Figure here*" - ] - }, - { - "cell_type": "markdown", - "id": "b1c62590", + "id": "15441e66", "metadata": {}, "source": [ "## Set up the mesh" @@ -67,23 +75,23 @@ }, { "cell_type": "markdown", - "id": "a8b6f7ba", + "id": "da642525", "metadata": {}, "source": [ - "As in any finite element code, first we need to set up the geometry of the problem with a mesh. Optimism can read in meshes from the Exodus II format, but it also provides utilities for generating structured meshes for simple problems like this one. Let's generate a rectangular mesh with width `w` and length `L`. We'll also specify that it has 5 nodes in the $x$-direction (hence 4 elements) and 15 nodes in the $y$-direction." + "As in any finite element code, first we need to discretize the geometry of the body with a mesh. Optimism can read in meshes from the Exodus II format, but it also provides utilities for generating structured meshes. The structured mesh route is often best for for simple problems like this one. Let's generate a rectangular mesh with width $w = 0.1$ and length $L = 1.0$. We'll also specify that it has 5 nodes in the $x$-direction (hence 4 elements) and 15 nodes in the $y$-direction." ] }, { "cell_type": "code", "execution_count": 2, - "id": "b7750f13", + "id": "4938c5c6", "metadata": {}, "outputs": [], "source": [ "# set up the mesh\n", "w = 0.1\n", "L = 1.0\n", - "nodesInX = 5 # must be odd\n", + "nodesInX = 5 # must be odd for element boundaries to contain material boundary\n", "nodesInY = 15\n", "xRange = [0.0, w]\n", "yRange = [0.0, L]\n", @@ -92,16 +100,16 @@ }, { "cell_type": "markdown", - "id": "f18bd7ea", + "id": "619b23e4", "metadata": {}, "source": [ - "The mesh object has an attribute called `blocks` which lets us group elements together to specify different materials on them. Let's see what our mesh has:" + "The mesh object has an attribute called `blocks` which groups elements together to specify different materials on them. Let's see what blocks our mesh has:" ] }, { "cell_type": "code", "execution_count": 3, - "id": "bb1bd192", + "id": "219257f7", "metadata": {}, "outputs": [ { @@ -127,7 +135,7 @@ }, { "cell_type": "markdown", - "id": "2bea405c", + "id": "5afe9f12", "metadata": {}, "source": [ "`blocks` is a standard Python dictionary object (i.e. a `dict`), where the keys are strings that lets us give meaningful names to the element blocks. The values are Jax-numpy arrays that contain the indices of the elements in the block. The `construct_structured_mesh` function returns a mesh with just one block. We want two equal-sized blocks for our problem like in the figure, so let's modify the mesh.\n", @@ -138,7 +146,7 @@ { "cell_type": "code", "execution_count": 4, - "id": "ec126372", + "id": "77813b01", "metadata": {}, "outputs": [], "source": [ @@ -148,7 +156,7 @@ }, { "cell_type": "markdown", - "id": "d88821e7", + "id": "450d807f", "metadata": {}, "source": [ "We'll loop over the element node connectivity table, which is the `mesh.conns` object, extract the coordinates of its vertices, and use the `triangle_centroid` object on them to determine if the element is in the left or right layer. We will store the results in a list called `blockIds`:" @@ -156,13 +164,14 @@ }, { "cell_type": "code", - "execution_count": 26, - "id": "a35d7cdd", + "execution_count": 33, + "id": "70c7dced", "metadata": {}, "outputs": [], "source": [ "# initialize the block IDs of all elements to a dummy value (-1)\n", "blockIds = -1*onp.ones(Mesh.num_elements(mesh), dtype=onp.int64)\n", + "\n", "for e, t in enumerate(mesh.conns):\n", " vertices = mesh.coords[t,:]\n", " if triangle_centroid(vertices)[0] < w/2:\n", @@ -175,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "b87eed4b", + "id": "960c5019", "metadata": {}, "source": [ "This will mark an element as block 0 if it's centroid is on the left hand side of the slab, and as block 1 if it's on the other side. Now, let's make the `dict` that we want to attach to the mesh object." @@ -184,7 +193,7 @@ { "cell_type": "code", "execution_count": 6, - "id": "e245933d", + "id": "2f8ce025", "metadata": {}, "outputs": [], "source": [ @@ -194,16 +203,16 @@ }, { "cell_type": "markdown", - "id": "61954ea8", + "id": "24f443e0", "metadata": {}, "source": [ - "Now we can make use of a function that takes in the original mesh (with one block) and the block IDS we just created, and returns a new mesh that is the same as the old one except that the blocks have been updated." + "Now we can make use of a function that takes in the original mesh (with one block) and the block IDs we just created, and returns a new mesh that is the same as the old one except that the blocks have been updated." ] }, { "cell_type": "code", "execution_count": 7, - "id": "e2a1c548", + "id": "7b5c94f6", "metadata": {}, "outputs": [], "source": [ @@ -212,7 +221,7 @@ }, { "cell_type": "markdown", - "id": "bcb47738", + "id": "35f85f8a", "metadata": {}, "source": [ "Let's check to make sure this process worked. To do this, we'll make use of optimism's output facilities. Optimism provides a class called `VTKWriter`, that, as the name suggests, writes data to VTK files which can be visualized in ParaView (and several other visualization tools). First, we instantiate a VTKWriter object, giving it the base of the filename (the name that will be suffixed with the \".vtk\" extension)." @@ -221,7 +230,7 @@ { "cell_type": "code", "execution_count": 8, - "id": "ea13862b", + "id": "da768ca5", "metadata": {}, "outputs": [], "source": [ @@ -230,7 +239,7 @@ }, { "cell_type": "markdown", - "id": "619fe562", + "id": "067f340c", "metadata": {}, "source": [ "Next, we can start adding fields. In this case, we only have one field - the block ID numbers - which is a scalar field of integer type. Finally, once our data is loaded into the writer, we call the `write()` method to tell it to write the VTK file to disk." @@ -239,7 +248,7 @@ { "cell_type": "code", "execution_count": 9, - "id": "42104d3f", + "id": "19b50965", "metadata": {}, "outputs": [], "source": [ @@ -251,15 +260,15 @@ }, { "cell_type": "markdown", - "id": "607566b1", + "id": "cdee4b9d", "metadata": {}, "source": [ - "This is what we get when we open Paraview and visualize the `block_id` field. If you try it, you should see something similar." + "This is what we get when we open the file `check_problem_setup.vtk` in ParaView and visualize the `block_id` field. If you try it, you should see something similar." ] }, { "cell_type": "markdown", - "id": "50aa988b", + "id": "bdc4fb70", "metadata": {}, "source": [ "![paraview shwoing blocks in mesh](blocks_elements.png)" @@ -267,25 +276,25 @@ }, { "cell_type": "markdown", - "id": "f0637278", + "id": "da15ed2a", "metadata": {}, "source": [ - "Success! This is what were shooting for." + "The left layer is blue (element block 0) and the right layer is red (element block 1). Success! This is what were shooting for. (The particular colors may be different in your ParaView installation depending on the color palate chosen for contour plots, but the spatial layout of the blocks should be the same)." ] }, { "cell_type": "markdown", - "id": "bddeee2a", + "id": "6ddc1328", "metadata": {}, "source": [ "## Essential boundary conditions\n", - "We're going to make one more modification to the mesh. Looking again at the problem setup figure, we can see that we need to apply boundary conditions on the bottom, left, and top boundary edges of the slab. Similar to the `blocks` attribute, `mesh` has a `nodeSets` dictionary that maps names to sets of nodes. We'll make the nodesets we need by performing range queries over the nodal coordinates:" + "We're going to make one more modification to the mesh. Looking again at the problem setup figure, we can see that we need to apply essential boundary conditions on the bottom, left, and top boundary edges of the slab. We must create node sets that group the nodes on these edges together so that we can set the boundary condtions on them later. Similar to the `blocks` attribute, `mesh` has a `nodeSets` attribute that is a dictionary mapping names of node sets to the indices of the nodes in the set. We'll make the node sets we need by performing range queries over the nodal coordinates, storing the indices of the nodes on the left edge under the key \"left\", etc." ] }, { "cell_type": "code", "execution_count": 10, - "id": "4d67689c", + "id": "50df5e33", "metadata": {}, "outputs": [], "source": [ @@ -301,7 +310,7 @@ }, { "cell_type": "markdown", - "id": "413f82f9", + "id": "244059e7", "metadata": {}, "source": [ "Now we're going to register the essential boundary conditions so that the optimism solvers will know how to deal with them. This is done with an `EssentialBC` object. Each `EssentialBC` represents a boundary condition on one field component of a node set. As en example, let's create one to represent the $x$-component displacement boundary condition on the nodes of the left edge:" @@ -310,7 +319,7 @@ { "cell_type": "code", "execution_count": 11, - "id": "ebedd68d", + "id": "8df8bbf2", "metadata": {}, "outputs": [], "source": [ @@ -319,16 +328,16 @@ }, { "cell_type": "markdown", - "id": "3ffe98a4", + "id": "1b5f25c8", "metadata": {}, "source": [ - "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` we just created and create all of the essential boundary conditions in one statement:" + "This is one boundary condition; we have three essential boundary conditions in total to apply. The thing to do is to collect them all in a list. So the code looks like it would in a real application, we'll ignore the `ebcLeft` object above and create all of the essential boundary conditions in one statement:" ] }, { "cell_type": "code", "execution_count": 12, - "id": "70344154", + "id": "a8f61e3f", "metadata": {}, "outputs": [], "source": [ @@ -339,16 +348,16 @@ }, { "cell_type": "markdown", - "id": "9f671f9e", + "id": "91ac4490", "metadata": {}, "source": [ - "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions. It's purpose will become clearer when we use it later." + "Next, we create a `DofManager` object. What is this for? It's a class to help us index into our nodal arrays, keeping track of which degrees of freedom have essential boundary conditions, and which are unrestrained (hence are unknowns to be solved for). Its usage will become clearer later." ] }, { "cell_type": "code", "execution_count": 13, - "id": "a34991d4", + "id": "01f306eb", "metadata": {}, "outputs": [], "source": [ @@ -358,18 +367,18 @@ }, { "cell_type": "markdown", - "id": "0c4b89c1", + "id": "d0291001", "metadata": {}, "source": [ - "The variable `fieldShape` tells `dofManager` what the array of the unknowns will look like. In solid mechanics, the unknown field is the displacement, which happen to have the same shape as the nodal coordinates, so the `mesh.coords` array is a convenient place to grab this information. You could also manually enter `(nNodes, 2)`, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", + "The variable `fieldShape` tells `dofManager` what the shape of the array of DOFs will be. In solid mechanics, the degrees of freedom are nodal displacements, which have the same shape as the nodal coordinates `mesh.coords`. Thus, line 1 above is a convenient way to obtain this information. You could also manually enter `(nNodes, 2)` for the second argument of the `DofManager` constructor, where `nNodes` is the total number of nodes in the mesh, and 2 is the number of components of displacement in plane strain.\n", "\n", - "We use the vis tools again to check that we've done the essential boundary condition specification correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`. This is just an array of size `fieldShape` that indicates whether each dof has an essential boundary condition. Numpy can cast this to an array of integers (more precisely, the `int64` type in numpy) with values 0 or 1 which can be plotted in Paraview. The `dataType` argument is different now; for block ID, it was a scalar field, but for boundary conditions, we want it to be a vector field (one component for each displacement component)." + "We use the vis tools again to check that we've specified the essential boundary conditions correctly. First, we'll use a little trick to turn the boundary conditions into a viewable nodal field. We'll use an attribute of the `DofManager` class called `isBc`, which is a boolean array of shape `fieldShape` that indicates whether each DOF has an essential boundary condition. We cast this to an array of integers (more precisely, the `int64` type in numpy) which gives them values of 0 (for no essential boundary condition) or 1 (for an essential boundary condition) which can be plotted in ParaView. The `dataType` argument is different from the `add_nodal_field` call for the element block ID, since that was a *scalar* field. For boundary conditions, we want it to be a *vector* field (one component for each displacement component)." ] }, { "cell_type": "code", "execution_count": 14, - "id": "3d92324c", + "id": "6ea61275", "metadata": {}, "outputs": [], "source": [ @@ -381,7 +390,7 @@ }, { "cell_type": "markdown", - "id": "1103ae8f", + "id": "8a3cc28d", "metadata": {}, "source": [ "Note that `writer` still refers to the same `VTKWriter` object as before. A VTKWriter object is always associated with the same filename, so when we add a new field and then call the `write()` method, it will overwrite the previous VTK file. Indeed, if you open `check_problem_steup.vtk`, you'll see that it now contains two output fields, \"block_id\" and \"bcs\".\n", @@ -390,9 +399,8 @@ ] }, { - "attachments": {}, "cell_type": "markdown", - "id": "ceeb6e96", + "id": "93df5456", "metadata": {}, "source": [ "![contour plots of boundary condition fields](boundary_conditions.jpeg)" @@ -400,7 +408,7 @@ }, { "cell_type": "markdown", - "id": "f4ee2543", + "id": "fb692616", "metadata": {}, "source": [ "The first plot shows all nodes with boundary conditions on the $x$ component of the displacement. We see that the left edge has a value of 1 (which means that the boundary condition flag is \"True\" there), and every other node has a value of 0, which means they are unrestrained. This is exactly what we want. The $y$ component plot also confirms that the top and bottom edges correctly have their boundary conditions. Of course, the top edge has an *inhomogeneous* boundary condition. We'll enforce the actual value of this boundary condition later." @@ -408,17 +416,17 @@ }, { "cell_type": "markdown", - "id": "f8cef770", + "id": "c02948e6", "metadata": {}, "source": [ "## Build the function space\n", - "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation." + "The next step is to create a representation of the discrete function space to help us do things like interpolate our fields and do calculus on them. The first ingredient we need is a quadrature rule. In optimism, quadrature rules are specified by the highest degree polynomial that they can exactly integrate. A smart way to do this is to choose it in relation to $p$, the polynomial degree of our interpolation, which we can obtain like this:" ] }, { "cell_type": "code", "execution_count": 15, - "id": "68a56fe0", + "id": "e6cc7a6c", "metadata": {}, "outputs": [], "source": [ @@ -427,16 +435,16 @@ }, { "cell_type": "markdown", - "id": "d0b6b351", + "id": "61bbbb95", "metadata": {}, "source": [ - "In the linearized theory of solid mechanics, the natural trial/test function space is $H^1$, because the operator contains products of gradients of the displacement. Since our interpolation is of degree $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" + "In the linearized theory of quasi-static solid mechanics, the natural trial/test function space is $H^1$, because the operator contains inner products of gradients of the displacement. Since our interpolation uses polynomials of degrees up to $p$, the displacement gradient is of degree $p-1$, and the inner product of gradients is of degree $2(p-1)$. Thus, we choose this as our quadrature rule precision to avoid under-integrating our operator:" ] }, { "cell_type": "code", "execution_count": 16, - "id": "cf915ada", + "id": "eb4e827a", "metadata": {}, "outputs": [], "source": [ @@ -446,18 +454,18 @@ }, { "cell_type": "markdown", - "id": "fbb9d8e4", + "id": "6b967ae4", "metadata": {}, "source": [ - "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is a good rule of thumb.\n", + "The benefit of specifying our quadrature rule in this way is that if we later decide to modify the mesh to use higher-order elements, the quadrature rule will be updated automatically. This helps keep us out of trouble with hourglassing problems. Note that our operator is *nonlinear*, so the quadrature rule won't integrate it exactly, but the accuracy requirement of the linearized theory is still a good rule of thumb.\n", "\n", - "With the quadrature rule (and the mesh), we can now construct the function space:" + "With the quadrature rule and the mesh, we can now construct the function space:" ] }, { "cell_type": "code", - "execution_count": 17, - "id": "6c36cd8b", + "execution_count": null, + "id": "fcc72db5", "metadata": {}, "outputs": [], "source": [ @@ -466,15 +474,15 @@ }, { "cell_type": "markdown", - "id": "f79a4bc8", + "id": "c570688e", "metadata": {}, "source": [ - "We'll see this object in operation later when we set up our energy function." + "The function space holds data like the values of shape functions and their gradients at all of the quadrature points in the mesh. We'll see this object again later when we set up our energy function." ] }, { "cell_type": "markdown", - "id": "d44fd1a0", + "id": "dffeee71", "metadata": {}, "source": [ "## Material models\n", @@ -483,8 +491,8 @@ }, { "cell_type": "code", - "execution_count": 18, - "id": "2c1f89f0", + "execution_count": null, + "id": "4f7d7828", "metadata": {}, "outputs": [], "source": [ @@ -496,15 +504,15 @@ }, { "cell_type": "markdown", - "id": "66eab9e6", + "id": "81d32097", "metadata": {}, "source": [ - "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in optimism/materials." + "TODO: The property names are not documented yet. For now, you can find them by inspecting the code in `optimism/materials`." ] }, { "cell_type": "markdown", - "id": "176065d4", + "id": "c9a4acb7", "metadata": {}, "source": [ "Now we'll instantiate the other material model for the right-hand side. We'll pick a $J_2$ plasticity model, which is a bit more interesting (and thus has more parameters)." @@ -512,8 +520,8 @@ }, { "cell_type": "code", - "execution_count": 19, - "id": "24778eaf", + "execution_count": null, + "id": "d5d64680", "metadata": {}, "outputs": [], "source": [ @@ -531,7 +539,7 @@ }, { "cell_type": "markdown", - "id": "19011b92", + "id": "34971c24", "metadata": {}, "source": [ "The meaning of the parameters is clear from the keys. There are several hardening models currently available, such as linear hardening, a version of power law hardening, and the Voce exponential saturation law. We'll keep it simple and just use linear hardening.\n", @@ -541,8 +549,8 @@ }, { "cell_type": "code", - "execution_count": 20, - "id": "aa2b0a4e", + "execution_count": null, + "id": "e6040db7", "metadata": {}, "outputs": [], "source": [ @@ -551,7 +559,7 @@ }, { "cell_type": "markdown", - "id": "3f1723e8", + "id": "0eb8be55", "metadata": {}, "source": [ "## Write the energy function to minimize" @@ -559,7 +567,7 @@ }, { "cell_type": "markdown", - "id": "94039b0f", + "id": "47989532", "metadata": {}, "source": [ "Numerical solutions to PDEs are obtained in optimism by minimizing an objective function, which may be thought of intuitively as en energy. In fact, for hyperelastic materials, the objective function *is* the total energy. For history-dependent materials, one can often write an incremental functional that is minimized by the displacement field that carries the body from one discrete time step to the next. We're going to write the function for our problem now.\n", @@ -569,8 +577,8 @@ }, { "cell_type": "code", - "execution_count": 29, - "id": "e20f5769", + "execution_count": 34, + "id": "637b03bf", "metadata": {}, "outputs": [], "source": [ @@ -582,21 +590,23 @@ }, { "cell_type": "markdown", - "id": "35578c4a", + "id": "39c60c2d", "metadata": {}, "source": [ - "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action).\n", - "To use this new function, we'll invoke it like this:\n", + "The `create_multi_block_mechanics_functions` writes some functions for us to help write the problem in the right form. The most important part of the output is `mechanicsFunctions.compute_strain_energy(...)`, which writes all the loops over the blocks, elements, and quadrature points that you would need to compute the energy from the nodal displacements. (Now we finally see the `FunctionSpace` object `fs` in action). To be able to write these loops, `create_multi_block_mechanics_functions` needs to know we are working in plane strain, and it needs to use the material models as well. \n", + "\n", + "The `mechanicsFunctions.compute_strain_energy(...)` function is invoked like this:\n", "```python\n", "mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", "```\n", - "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point. It has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (NOTE: For multi-material simulations, this is currently sized such that `nint` is the *maximum* number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency).\n", + "where `U` is the array of the degrees of freedom (the nodal displacements), and `internalVariables` is an array of the internal variables at every quadrature point; that is, it has shape `(ne, neqp, nint)`, with `ne` being the number of elements, `neqp` the number of quadrature points per element, and `nint` the number of internal variables for the material model. (*NOTE: For multi-material simulations, this is currently sized such that `nint` is the maximum number of internal variables over all materials, so that the array is padded for materials using fewer internal variables. We may remove this restriction in the future to improve memory efficiency*).\n", + "Under the hood, this function tells optimism everything it needs to solve the problem through automatic differentiation. The derivative of `mechanicsFunctions.compute_strain_energy(U, internalVariables)` with respect to the nodal displacements `U` gives the residual function, and the Hessian yields the stiffness matrix (more precisely, the Hessian-vector product is evaluated to give the action of the stiffness matrix on a vector). All of these derivatives are taken automatically by optimism in the solver, so we don't need to worry about them at the application level.\n", "\n", "All of the minimization solvers in optimism require the objective function to have a signature like this:\n", "```python\n", "energy_function(Uu, p)\n", "```\n", - "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", + "where `Uu` is the array of unknowns, and `p` is called the parameter set. The parameter set essentially holds all of the information that is needed to specify the problem but *isn't* the set of unknowns. These are things like values of the boundary conditions and the internal variables, as well as some other categories. We need to cast our energy function in this form. Let's write it like this and work out what the intervening steps must be:\n", "\n", "```python\n", "def energy_function(Uu, p):\n", @@ -604,15 +614,15 @@ " internalVariables = p.state_data\n", " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", "```\n", - "On the first line, we use the parameter set to pass from the array of unknowns to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", + "On the first line, we use the parameter set to extend the array of unknown displacements to the full array of nodal displacements. This means we need to fill in the values of the inhomogeneous boundary conditions, which is what we'll do when we implement `create_field(...)`. Next, we pull out the internal variables from the parameter set. Finally, we use the canned `compute_strain_energy(...)` function with these variables in order to compute the total energy.\n", "\n", "The inhomogeneous boundary condition part is handled like so:" ] }, { "cell_type": "code", - "execution_count": 34, - "id": "217248f1", + "execution_count": null, + "id": "379d4116", "metadata": {}, "outputs": [], "source": [ @@ -627,73 +637,80 @@ }, { "cell_type": "markdown", - "id": "c578d78c", + "id": "559b6254", "metadata": {}, "source": [ - "We will store the applied displacement in the first slot of the parameter set. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", + "We will store the applied displacement in the first slot of the parameter set, `p[0]`. In line 4 above we extract it. Then we make an array of the same size as the nodal displacements, set it to zero, and replace the values in the DOFs on the top edge with the value of the applied displacement.\n", "\n", - "Now we can write the `create_field` function shown above in the proposed `energy_function`:" + "Now we can write the `create_field(...)` function shown above in the proposed `energy_function(...)`:" ] }, { "cell_type": "code", - "execution_count": 36, - "id": "2755c6db", + "execution_count": 35, + "id": "3a23d5ce", "metadata": {}, "outputs": [], "source": [ "# helper function to go from unknowns to full DoF array\n", "def create_field(Uu, p):\n", - " return dofManager.create_field(Uu, get_ubcs(p))\n", + " Uebc = get_ubcs(p)\n", + " return dofManager.create_field(Uu, Uebc)\n", " \n" ] }, { "cell_type": "markdown", - "id": "34511c5d", + "id": "fab84ba9", "metadata": {}, "source": [ - "The pieces are finally in place to define the energy function for our problem:" + "The `create_field(...)` method in the `dofManager` takes in the unknown displacements `Uu` and an array of the same size with the values of the essential boundary conditions `Uebc` and packages them up together to create an array of all DOF values. The pieces are finally in place to define the energy function for our problem:" ] }, { "cell_type": "code", "execution_count": 37, - "id": "be482478", + "id": "6ebd34c5", "metadata": {}, "outputs": [], "source": [ "# write the energy to minimize\n", "def energy_function(Uu, p):\n", " U = create_field(Uu, p)\n", - " internalVariables = p.state_data\n", - " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n", - "\n", - "\n" + " internalVariables = p[1]\n", + " return mechanicsFunctions.compute_strain_energy(U, internalVariables)\n" ] }, { "cell_type": "markdown", - "id": "5fab41f4", + "id": "7d271df0", + "metadata": {}, + "source": [ + "The slot `p[1]` is always reserved for the internal variable field." + ] + }, + { + "cell_type": "markdown", + "id": "2815942c", "metadata": {}, "source": [ "## Set up the optimization solver\n", - "We have an objective function - `energy_function`, which we will hand to a routine that will find the unknowns displacements that minimize it. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", + "We have an objective function - `energy_function(...)`, which we will hand to an optimization routine that will find the unknowns displacements. In this section, we specficy which optimization solver we want to use, and tell it how to work. \n", "\n", - "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization on the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" + "We will use the Steighaug trust region method. This method uses linear conjugate gradient iterations as part of its algorithm, which in turn need to be preconditioned in order to be effective. Currently, the only available preconditioner in optimism is a Cholesky factorization of the stiffness matrix. We need to intruct the solver how to assemble the stiffness matrix like this:" ] }, { "cell_type": "code", "execution_count": 38, - "id": "cd2bb5dd", + "id": "396fcac4", "metadata": {}, "outputs": [], "source": [ "# Tell objective how to assemble preconditioner matrix\n", "def assemble_sparse_preconditioner_matrix(Uu, p):\n", " U = create_field(Uu, p)\n", - " internalVariables = p.state_data\n", + " internalVariables = p[1]\n", " elementStiffnesses = mechanicsFunctions.compute_element_stiffnesses(U, internalVariables)\n", " return SparseMatrixAssembler.assemble_sparse_stiffness_matrix(\n", " elementStiffnesses, mesh.conns, dofManager)\n" @@ -701,65 +718,63 @@ }, { "cell_type": "markdown", - "id": "2b5da4b6", + "id": "50fbde23", + "metadata": {}, + "source": [ + "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses(...)` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "7aeca02f", + "metadata": {}, + "outputs": [], + "source": [ + "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n" + ] + }, + { + "cell_type": "markdown", + "id": "6e4b621c", "metadata": {}, "source": [ - "We see once again that `mechanicsFunctions` provides a helper function. In this case, the function `compute_element_stiffnesses` takes the same inputs as the energy function, but instead of returning the total energy, it returns an array containing the stiffness matrix for each element. The elemental stiffness matrices are the Hessians of the total energy in each element, and automatic differentiation is again used to perform this calculation. The `assemble_sparse_stiffness_matrix(...)` function takes these elemental stiffness matrices and contructs the system's stiffness matrix using a sparse matrix data structure (from SciPy). We tell the solver how to use this capability by building something called a `PrecondStrategy`: " + "Finally, we can specify custom settings for the solver is we wish." ] }, { "cell_type": "code", "execution_count": 40, - "id": "b1bc05ed", + "id": "01556d5f", "metadata": {}, "outputs": [], "source": [ - "precondStrategy = Objective.PrecondStrategy(assemble_sparse_preconditioner_matrix)\n", - "\n", "# solver settings\n", "solverSettings = EquationSolver.get_settings(use_preconditioned_inner_product_for_cg=True)\n" ] }, { "cell_type": "markdown", - "id": "d4b55505", + "id": "dbf08b3f", "metadata": {}, "source": [ - "Finally, we can specify custom settings for the solver is we wish." + "## Solve!" ] }, { "cell_type": "markdown", - "id": "3e10766a", + "id": "031e405b", "metadata": {}, "source": [ - "## Solve!" + "Parameter sets are constructed by calling the `Params` function in the `Objective` module. This is to help organize them in certain categories that the solver needs to be aware of, such as which ones have derivatives taken with respect to inside the solver." ] }, { "cell_type": "code", - "execution_count": 30, - "id": "28b893dd", + "execution_count": null, + "id": "c18895d2", "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Assembling preconditioner 0\n", - "min, max diagonal stiffness = 0.9999999999999998 1.0000000000000002\n", - "Factorizing preconditioner\n", - "num warm start cg iters = 1\n", - "Assembling preconditioner 0\n", - "min, max diagonal stiffness = 0.9875365254464693 1.0110599625088337\n", - "Factorizing preconditioner\n", - "\n", - "Initial objective, residual = 0.04170630570636782 1.5580526760827706e-05\n", - "obj= -1.47133479e-10 , model obj= -1.47130159e-10 , res= 5.10352716e-10 , model res= 7.43900588e-21 , cg iters= 1 , tr size= 2.0 , interior , accepted= True\n", - "\n" - ] - } - ], + "outputs": [], "source": [ "# initialize unknown displacements to zero\n", "Uu = dofManager.get_unknown_values(np.zeros(fieldShape))\n", @@ -783,7 +798,7 @@ }, { "cell_type": "markdown", - "id": "2b1968f3", + "id": "782f1ac5", "metadata": {}, "source": [ "## Post-process\n", @@ -792,8 +807,8 @@ }, { "cell_type": "code", - "execution_count": 31, - "id": "e14a722a", + "execution_count": 41, + "id": "f5ac29dc", "metadata": {}, "outputs": [], "source": [ @@ -801,22 +816,23 @@ "writer = VTKWriter.VTKWriter(mesh, baseFileName='uniaxial_two_material')\n", "\n", "U = create_field(Uu, p)\n", - "writer.add_nodal_field(name='displacement', nodalData=U, fieldType=VTKWriter.VTKFieldType.VECTORS)\n", + "writer.add_nodal_field(name='displacement', nodalData=U, \n", + " fieldType=VTKWriter.VTKFieldType.VECTORS)\n", "\n" ] }, { "cell_type": "markdown", - "id": "da5415dd", + "id": "be97ee7c", "metadata": {}, "source": [ - "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial conditions of the internal variables. That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new displacement field and the old internal variables as inputs:" + "Let's also show an example of plotting quadrature field data. A commonly needed output is the stress field. We make use of another function in `mechanicsFunctions` to help us. However, before we write out any quadrature field data, we should update the internal variables. Currently, `state` still holds the initial values of the internal variables (before the load step was taken). That is, the solver finds the equilibrium displacement field, but doesn't change `state` in place. This is again due to Jax's functional programming paradigm. The following function call returns the updated internal variables using the new equilibrium displacement field and the old internal variables as inputs:" ] }, { "cell_type": "code", "execution_count": null, - "id": "e46298c0", + "id": "d2926f79", "metadata": {}, "outputs": [], "source": [ @@ -827,7 +843,7 @@ }, { "cell_type": "markdown", - "id": "d23fc4a2", + "id": "57ea016b", "metadata": {}, "source": [ "We make use of another of the `mechanicsFunctions` member functions to get the stress field, using the updated internal variables. Note that we don't pay the cost of iteratively solving the history-dependent material model response again. This function expects that the internal variables are already updated when it is called, and simply uses theses values in the evaluation of the material model energy density functions." @@ -836,7 +852,7 @@ { "cell_type": "code", "execution_count": null, - "id": "d3b3e743", + "id": "dc41d135", "metadata": {}, "outputs": [], "source": [ @@ -846,18 +862,18 @@ }, { "cell_type": "markdown", - "id": "b4221fcc", + "id": "b2cae332", "metadata": {}, "source": [ - "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function coputes them both.\n", + "As the left-hand side of this assignment implies, the scalar energy density field (at every quadrature point) is returned in addition to the stresses. In fact, only the scalar-valued energy density function for each material is implemented in each material model. The stress field is the derivative of the energy density function with respect to the displacement gradient, and optimism uses Jax's automatic differentiation to compute this derivative. When computing this derivative, evaluating the energy density at the same time is essentially free, so the above function computes them both.\n", "\n", - "The `stresses` array contains the stress tensor for every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, this must be remedied, either on the physics solver side or in the visualization software. One simple way to accomplish this is to compute element averages of the field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this:" + "The `stresses` array contains the stress tensor at every quadrature point of every element. There is no interpolation associated with a quadrature field, so in order to visualize the stresses, one must be created, either on the application side or in the visualization software, if it has this capability. One simple way to accomplish this is to compute element-wise averages of the quadrature field and plot it as a cell-based field (instead of a nodal field). Optimism provides a method to do this called `FunctionSpace.project_quadrature_field_to_element_field(...)`:" ] }, { "cell_type": "code", "execution_count": null, - "id": "e0d3a7bb", + "id": "55774221", "metadata": {}, "outputs": [], "source": [ @@ -872,22 +888,30 @@ }, { "cell_type": "markdown", - "id": "2558d9ca", + "id": "1f3b4a7e", "metadata": {}, "source": [ - "The above code can be wrapped in a loop that will compute the response at multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" + "This concludes the tutorial. As an enhancement, the load stepping part above can be wrapped in a loop that will compute the response over multiple time steps. Note that before moving to the next time step, you should update the internal variables in the parameter set:" ] }, { "cell_type": "code", - "execution_count": 33, - "id": "7e0acf41", + "execution_count": 42, + "id": "6fcb6668", "metadata": {}, "outputs": [], "source": [ "# update the state variables in the new equilibrium configuration\n", "p = Objective.param_index_update(p, 1, state)" ] + }, + { + "cell_type": "markdown", + "id": "782062d0", + "metadata": {}, + "source": [ + "You can look at other example simulations to see how to do this." + ] } ], "metadata": {