From 6e682e1e719e0f1c855bb6aa1cd0915e66842eba Mon Sep 17 00:00:00 2001 From: iasonmanolas Date: Thu, 6 May 2021 16:45:52 +0300 Subject: [PATCH 1/2] Update ground of truths to encompesate for changing the rotational components of the displacements of the simulation results to euler angles. --- .../LongSpanGridshell_displacements.eigenBin | Bin 6736 -> 6736 bytes .../ShortSpanGridshell_displacements.eigenBin | Bin 1024 -> 1024 bytes .../SimpleBeam_displacements.eigenBin | Bin 256 -> 256 bytes 3 files changed, 0 insertions(+), 0 deletions(-) diff --git a/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin b/formFinder_unitTestFiles/LongSpanGridshell_displacements.eigenBin index cd8534c3a97b8327ea2848a12da66fca9b764047..45332e0fb3be7c09ee3de9f11ab5439f68192ba8 100644 GIT binary patch literal 6736 zcmb`Jc{r5q`~L|uH%Vp9E+j2Pk|mGpMu=7+WzSYmBwNT@5|wq3HTzE4%2t`{Mr27@ zBI|^ZvPDWOziWQKU!UW9JUxGXkI($`KIWYBe9d*g&-1)y!WsVGi}e5Xzb%=qlleNd zEYQ^M+pTGv4D{V%z0Cw(+xLOd8ctJz-XKNAhdS5u$ecT9%uj<=Ir$FpEWN47&{LM`|p-7b01vIH1i zyE5f^J{|ZMKWvk74hKo&M6Lw0wqsUtZTaJ)n6xCPz)kmpi@|Vy3F-( zc#++}io-n}Qlps?_u7^~=DO3qr;=-7b7NG-2IgkasaDkQE?Wa`UEs3n*3JX^eS?Da z2Pwb;Gib}KbO8@Im%`87V}i9KL5fZG&!}IthdcIhZjO>8q=wi`Z{=FU9^W> z=L@Lte_uVoe(6l4^Q}5?RDY+|@}?qixcK7xZ7>$Nsfj!|K$IY zQ~g2rhupeTYajBbQTmy%1@dQ@F)>~i`SZ{uAdc?O(xp&RJ@O~T`13Ky->h!debRC( z@+aS2v%3uWV;6mF07L#P#7OSCi~Law+Sg0>hv$-pA_wX(i-bvw2kI}!Lw4^T)Zf=E zI#La&zk&J~O9^^SQgOzgQd1=Zh6pP7gGnr9LaY(r7+sTYl9|qWOH$zNe@Z z&F8zd%HYhFdi9Fgs9Ym7pRpHK9>}5jwEqDaH>3F!=i{~5hvu`2E#jUDnos$%H=ep^ zzpBlz9Xp8jYcPaa_ZZr*P=jz>744V&#sPU&v|o)C&b3qAWDwJN>eK$K+)&u}$N?_2 zU#@qKu1ul*YS=FJQvmH(R_d#QQnX*ImE*Z>Xutm3_W>sHe&OqA7Vu<`c(*;N7RV_- zuK0lb55=0c>2$sTzvuP0-qA3%*GBPIcH&%*P#(q0kpF68E#R-9g)t+}@AF-5vebWm zJ{jcooaUfdmjKS2>?gpU8t_e(7q^9V;7PT!dq!9ts4AFuDnv@EzOhjR#aH~^v{a$K zPo$=XPoVyM+9g^2%Rc!(+~0iDJC&?|pNa-QT0L@woYR50GN&KIp*%1{%u3CN76PMp z3L;V5MZgO<_}oFtaBk2$0mapyHY&>M<^kVKcUv3BbfE5NbYa*q`VW5gtNS?yBEdj} z=lho+-*6DIFnHHdp919fsQKhC$AASvktx-tSU_@70nR{{)6w3cpxpn#T~S>>nwlvbpeh$6~J=43rJw|9x#b<1sWR8ntSHl zfNrSJ!fB*W+9N&HQT!8AkCH0N`+P!;h1Cbusf7087v2BhdrTPq7-8E4UbahAePmgo ziaMQ$m=;^2x~iP{Vqm{QJ&-vTFH%*9i?C%BADmac6L4%v{lm2{gQmMy)CCn?81Mbu2!HqWw+7&9A{H+0 zP`Y6Kkv3l!ILF5a>)Tx*{MEW;Ym6)0?!1fR>#Q4`v+&~*pL2tIu@)&-6z4n0|DF}) z_3IWseRIVH3fM$yZ@uXLCqLfD02c{{+l5*?9eu-Lo4K%*r9K5(4EyCIF2}&J<=WU? zO|j6+z2M@0q<5LF6rV=%;ogJA6g>+3nLHbE*gG7$iruP+6$t)=A3-^)WPK_c_WDgv zZ*oqDe1XvCn{FP=xqi=ODzp%)Ye#(!;Vy#yx~HrnkhXS2evCr#)2FiBx1&6sNCs&W z$8^ZqY0Lb|F!~R^YyZ~9#`DP#K3TpLW?ceZ4n6%#ZBGq+`8;XS3f957D+x1q!|Gt( zgSVo>bs2eLI#3~Q8sV_5>zza*xAwom>Kzqc-zwmg~zO?SSzwJoU@ zT7UyN@0pw7Rl8|_wD-S%U;d+kv0GaXqWD0@#VxfcZ$4c4oqu92l*-_F!hkrx`}yBi zL6U1aXLSQ;Xz8;yW$gw>MLmDIKkNq=L(R1p=4imhO)QcYJPzVC*S#1)2PWpt`*hxK%FV{PTS5kJO~O(dVyC*gii=83o&&w)oI4O#xBwFD9qv zK7e&~L0^|T=Yad6@)5e;9^aK?j^dJqozBB3Pj29VQbhW6)2xIq<7xEc@8MsTh27QECRWam7J)_C2&e7@g`l*1$TQiqPVPPgrYIZQ~ePf za|PA8SRD#WK%9T_dp8+ob|W87jtL(>fzFjWJtNzS>epFy%#g&F!N-rykz?tffkV!O z8eNU$1msq?E`!&zcqxq=<3NE|bwHtO4(L5- z(N9C?EHkp5^C{8j6dDPMdk5ytkLU?@*9$O{WIaZlM1B z^ZXJkx_oqgtB7wOozHjK&XCThyxk;8=jT7I6QXn1Rv9PhI?N*KMdu3~Ic!Jglk7fI z==^=W2P6ND@93YCO!tS(ooYk(XNpBtpYD%C(5)?GdmK zH2N@%?$5MV^Aoy1TM`*=|C>K+9_?@G{Vh4F>PqkLWuf_F9;k2ZwB^0eP`wt#_$R%u zR_7Yt&^1f+>=Sx_n{TKm(EICj6n4=2yK~yL^56Qa-IQ^GKA&EXO83>FIlIt8ou$vG zOuOed`h4Qfvdy2+^B?lHEu!nu*Y4UDD1LNh_$_@t_3zVa==15d^W24ho6mgC%d_GDmK_jg>8 zr4OppwlaV8Z~OJne*bUZYYdoXo3k{6x`}knpXNQl?Onl&3>*Z`sgv5TE{}mk+b@pa zca4I$%rJ!?NYz!=2hX9nwe@z&X_VJ##mS^%(F3vuHfQZ6HG;b@TQ;9iX#ipyHgC?6 zYXrZ%tjk@?JAvb)?}@z|UxR%cFPdG2Lm(tXb6n1O2(T46Ka@qv!Z~3VfZ`XMo>{*@ zd2zbK9v4uZ4};(HUc_nP*0MH0T!RFyt!#lcV4%EBJ!!ZOxccY*z0leU)^Ye={3X@} zzO#|@>p6Qs)5WO1xm`VAL%c-!4AS}Q-WzgIyzH#)R1?ZuS?6J?1(AnV7GWPCan)GIDlA`q(4FI|jv% zO&rY@L3zrIoPudGb)es`KRNVm4Y0}QQWEwk15Yz#G(_0*K*-{Q#DmKB0K53ifM;6) z=qeV|QL!ihD$M@#w`U4~bdRV4drcw0Wu&46k)GhIm-ayMPQSzb5#|LTCtK=Lc1r;W z5&xxgcmF-$?T|6Y!IlS<_R4j%wPJV&HBe2U9K&-zKk@FOV))*EcCQcu!&}}(Djr4e zb>cSf6MT%}Ic_0e1dxhv_KT*+?@Y!?-a>i#f>O=>=>4+{_YpC~p}P6A^&qaSrBV`} zhT&2zG4{=-7|znCXzQnp;ek)9dKR~0_~z3DXEZm4KVo%CkK2IZjM4kc=&F=x^H&Ut zzZPiq_=fTXYPp>vQJv?xNx?UW)5X$dV~V(A(fKyD*%-cXb624K1q?rPeq1kdo{SIj zJd=&gB;$`H`z6+2CF2p-*M9BxCgV0y7)1u@yo&iR4itBkcIZezdBJCcV!a& z{}U2!u8}fA*MmDeJM&Q7Q(x!KILc#4nyfpA>RzeVJQhWq2p7if8xi*mAIA}gRWhz^ zmAzpReZOi@Hc!WTGCo{a5zT;;aMRv2`SnFCcrF$wtoV=xFF1Y7maftdF1XO+LaJxI z7*QTQ+fTxWG-{1LBAmk}C zPMgsP|H;dOPhIaarR(ZSTYDgiyI$M9YZJ-~-EAphgz8d<7$g%BrzLe%age0}9vFQ? zCD~U){|Tp@-9z`Ft&L_`ujDLZjlt(st`hL8?llE#r~LV4x~b=b&*7Z~_X2v13%gg>`&?f4&jyJh^4C zxVHh^hKaHvx;?OcH^X=TmQKigYnHjxs12^0H+$zd*9bF%&%dq@p+d?Fk=2T<HzM zC*#qlLhxu%G>$)%3TBrQll0~rK~daVR;zIva8>*9T&uMcXjYAXvexYZsx`8Xq7MdP z5p%Jbn8<7Ral5kn7Q-&MU0uX&uV)LC3E{r;w+R6&*|rQOOcp@52gZ{IQw5-ut(^U} z2?2IbFZy-xsBe)5IE3?pE!SU3^>#mpElkx3T>`e-KxDZ1Vi=K^o!(s;PNaN*V&vl z=x`H!_B>nGjX@pNKLnJ5H1=1%_PJyH+kJ@t!na@#;_GHb+{y*)sybaKx0 z$`Hsm&^I-c9tEFs*_reBScs&NtKv^9SctEF%`DB2Nkm>XGhXRMCRQr3w^|!8Vt^)L z@Z%9iJk9BK{{9Hq|P zX?$2M!-?TmrMbs*{m3}aH;bSvh_3~!B~RwD;5u%&c4jj(4ogniCA6~;ONvdnk~)b{ zBhDR2Tp94;BlgEtUz0=rlsfg<>_+~)HBR?!kiqcgzY|M%-;(i~cv@*1dX8m7uI_G( zjBEO`@7`ET!gr*HJuF-%6GG1CExl&Q#NQvsUN7y$h({bAoQ75y@t&6WB1RV@_^Y8= zUoJ-M84`HXk&EFH$uVjvx)?5Et-aI03d2hcz02lzVt8G%*nOuNGF~rJLdjT0_n>l5 zVJag=81CLE(Bps+U)596B1cpx+1L zmeI+C668->kzlU_h6_rl>%=jl?+u#j>*y?ko?bgVCAJsCv-{MK_OsatWbJMZ$Lsd~mxlX*ULs9a0r>WhcDw9$xTZ51%p?8+q>+C!i;PwAiV3j?ub zZzNRxa{$A1gmp`je?GTO4Edwm4(zZW<-QJ)3G%~edV zc4dMFl|rFbrDVWu#`bxmd@%6YUUFcwVlY(Mq>#;~k_>e$Hg_a-XF|1t3D&O-S@7Bw zj$Ll+BH-t%DRqVQ<m!lCK*QAS2PqhAGLZ}2_RI$e&qzbvd^-TCwszD) zJ1G6UIjrcJ4@|?}r6`kwVYs07*zF^caJ%K5Ctj2Ba3bnJv{zp?tgca;P3^|;2Ut0$ z^fHFe1x9}=VPeH48OKbQH?ZRB9uAEO+^l%gnBes;Ayz!q%5$<)h?V&M@!eGvC$?C} z3!^-O@xF8`sv}N|h1$q1V+1aN1!Q((gsS9xpG6pkhg{6{(x72DN6(%<-cbyX++e@@ zX&(LWP?>#4UB&R@B<)(%}HEzM&DTfk5T*C<9z$yzk( z(l7#a(kg7iFrvw1zH}-U{hxaCFl2$56@OvL)A*?m!>u!a7*-BqxYO!teb^j^n=hIi zY(vlEULm5@jLyqr#tUxEp*_2JWWH<=-8(^-JQw=V{oq|N^_7{G2vlMboR7tbEU%=W ztiBlDD>v^RKw`x|a$a&Gbz=D10ISmZb_}2Apbe@FVR%@C2~Q>qE6)Bxw)8FwE0H*u ze?oZ(Bi?)+53*^;h;<_S%osZ{!ej1k;dK%#;qRuV$nJ~oZO_G&d&r-9CI3}T# zEq~)Fh9}szk61Kd_^(9Bu3V3vZ)oH5-bDM`P#yw@-D28>`| z7+1Ub6eDQ(HTmPn{~z(Y77LI+dapVKYN{}z%QUmWvJxZYvsC1^Kf?&eXE_@1B}Vk} xoLVsZg%P~*J2CBF7>?b4&QSCcJ=Z!{9gb%h{=xESk4`0qr$(IHpIe3D{|8a>nu!1a literal 6736 zcmb`Lc{G;k|Hlhajby8gqO?s-`=-M9Xd|Xg+0sG_6)8n3X+*XP>9K{9r7T6I$n)gk zV;O{$q!N-PqGFIjq~GnjuI7(A^Uv>m{c{~x_x(Eed2gM=dl3B*XXpR=IMo^9w75D7 z)LthTj@=asb+c3rr*7~Ao2P4ZI;7nI@ARK=HL(M`_s?!_TTBIU8rwoE6? z$qu7KR^Oa6$qi5cs_LMz*$*8*T9-bw4#h1=(*2cnQK;-(?vKAkd=-Tw^y8; zU-AeSke{dUb!9_NdWrv|tzECMo`$a~gsq?#< zPBA<*hG`4Cy2i`rM1N0K7ZTm2Rn$McMG@phnpVoRbpUQh_J~gsvK9}!ICHWkuNAvk$AP$3PqUo#>%i}{U)Nc`l!Ix%!dbtvJ88dW-KPEeU+;l!+9;L& znkevGy*zt%hydR5-5$1Z{W?>r7NJw}2)^g_`+4U4`1e}e{g{a%K6kn2!prPuAH665 zpF<^SUtdST_wU2(Uqc`54TZH!YNiC6--O=iy0R5>_;5*X%o*uMfD!U3kpx$FifLeN7gj|8l?cjq0z=vo`XBi@MUSlMh@0 z-t810@zpUfL698%(K8s=L-Iwoujk){De;Lj7(!BsspnpT_m^WJF63DF0&9}|7BzKH09_e_XB$h%4OL7Oz94<7hP^g+p8M4zAdF250fanwTM zFXrVFe=&^ui=D5Dzj$aR@fVj?5PwnKnfQzA9uR+VqB-&RCw>ymC)(+ee4<4X$tUik z`NWy&B%gSP<`esmlYHWRnooQmO7e*TG@n25=jD)oVa8(8FMLn?g(GRj~J7Fp&acOCOVRS{mVZ}_W(cAJ;1F3vIiJR_u%{gFaMrr=^o%8bPurUDcJ*@ za**u7zufP>{;y7_7qCn*H5+a#3$nKjD}|IJvn7utm&2#~idz4)R0D5RhhMc@9q25% zFEem{Z7%=yTjyL6L=ZsqRpS_~1ZM7E=kp2b{5j5#e8-Kw(_eBr#2m?XG~m8#oVR3mA@}~Q zz#_4s6FWgaSmSN&l}{ik-{n!mokjZ`r2AAlA$|PpGqK!r4LVQF%09dUm;Bt!ryr_+ z;@7{B^I7Fy2j(%C{Fgpyf&H=LJR-R7T)TB{Wi!V!PYKX`J+%uS{8^k-&7FPjy%Vaq z`@1~RyLQ}jL9Jgmoa1zxyr}7k3+K;|-~XBMH!{90d3e1+aERTp5LYLkB@FJxg$gRj|*moie4ku|ls#%qy+_&v^+rmmP zLEZndaaaMaOI?~HEt!K?d3#rCUw@UuUq0EOmh9XZ*LJms?b-mE+4uhUv~)FOW?f| z$nA}K4iz68H{bPm3Dt8n+;(rQ!70<*+WQ+TFgL7fgi~QDw##49U)xrIit;Uuf`UvO zF2AokhWq}yg->PoG{?cQg+fo`!c4fZ_0iR%?FCTq_*C(C^Ys zuYoM}dks=gt8wrnbNAw4jz2ANQFLu7{&p}9cdV5J0zI(D;wJ13r)4ZQO zHcp9$LbL6rJC|m`rmDau-~EMfQ6gSXr>+!aGDm%}4Xyyi_LExco>YT`2IY4aQ+^ub zZ)f~J7=I<>JNZz43FB8)QvOlKpT_u$8DDBH)n^3Lr;O>deiYTGK#uC;D@XOYr%LsS zVfy^a^tsFQ8GMxLqs8=*s-*t_Oz%-`jF>hCS)?;_@JEAuy>`J2rA z^~xw>QPzhYUx zOj*A+vVOI(evS60{d$#6`!$O7OUj=1D~$DP4%>rY*dDm1&^>Tud+?0y!Q~#h2b06; z9vozQP|Ef|$o4=ZxbzyI+XDrQ3}+Lrw=H>d^*(a_^1-$cORfjLDbDv5N~_W8bdkvj z?mRO8p5jsN{ec5#A;0N(b zCHU9%oI|cr&#`{uy5q;*R-o6B4-<2^{LfOi7i4p1lx(%yQSQEP!J>=Dx##k#YGuA~ z`w_FfsOJ&KiGJ@=HdHPL@%S@dVwMmS<15}(M+))u%8ipFLJJTwgT44|CFpQAW0ob? z*WrR~j(uFOei3z+@8M2grwMOXarfmHTVq_f=N9BH9orUCfIGbA7F$IKv4489-{edo zI%|1beB<`DNs_N?$n9an;3=Db;dX1C_T6n&HTn4E+46bYxE;8)Biimd*H?#Adt<}7 z(}FL%X$*HC+Are?+;cfS67PD4=irBJQzSCyr=ZK07k&B$DR?gL~Wt$MgV^S0PxqTUH(>S+_+kx_EJ(W)G^uPG=$W!hFfyvy1o0e72>uIYFehZisG{NlPb z4kyJX{T+vqpNnuoE_ZxD)(=k!QHF5SKdZ$H>E{*({O(u5C_HWNNH`w>34#1K9VJxKWAvx@Lx>I1@uJ)VRQ6_gLg;|L%A zFe7~UO7($RR3BItO!R@zR3Ff%`oIvX51gd>fGO1nmQj7cp^E4O8>v3e z-87%Dh2|3$(R@O^2+1d8(0oD%%_rQT`GioKPdM{C$tSeXenA867u4I5enB_w7ev#3 zK?&^_d{HC)f(qI%5YT>sXf^2!7yPG$$-lW_?R z2JX&ycw>rkXs9zx=nmCU%nyL$ZLMB=o&`g1_l=Rp9U)LAl4w}!$piaBrRGJad9XVM z)kj5M!zb3sd-D!nMqyu<3~mWT@$SV7mZi9%Lsf}{N&OkHw3(uwBo+plXLpP04v&CG zBHG3|W>GMD`^FhwDp625xc}lCk0=B6k5K!4(vwBlv4CXr;z81ZF z8E>iG8FTek0G?aE@{>oED+-1SoT6Gzz+a!dUfkOq2%kR72O4HvhM|GA;;VF{L89z} zW1HtScz1i*8mqQ=h;YI8d)h)UeY2o^WJxgGTQ89DNQy*#lX%ngA%U2GZrN}jYj3<> zmaC#a=s2o-KOEEJX%B`X-OAWah#`NjNMAQhv%diz)Z7wOB?iU zVE?2nqt^TOupuqiS4H6@c<{v;G=dE;p0?y!pC-P!p8-ngpaQ42_JjPl?M2b`HJX+%?(5!kc%by zz^N-lACP%W^a1Hmq7PopCi@5Nbg32~LAaK5<+J$tP}VB>BXi3nZUtSx)kagA+(Tac&sNC)`2OFO*6o z{eq|aNWWlqJ?R&``+dnkzuknRj zM^I0KlbVai&JN0eb_cob_8cFq+OB)j;Qj#o#fQQ#D*XcnWiT_+KN%Y&Ve#bW)>K0% zm|L`e6y8k&Z+@2KM#y0a0pf^cToUn;DpKE0>$vtR->Y+nJ93VBNT}4w*dIPXx$mpC z-MUA5dxAQp{4TqI05#4&|65ah5Tx&;MEL~ZZ{JfJ)bosbtXj!?Q5v?^^LJrF#KYyy z;6kga=;`BGs;x#~KZ|y6jEVENxcP88CQ!tye$^w*nit#JK!W;==Hia~GPv{7XJWE( zZ3p!3mrjnF?ig+0%@VyLui62+)ax^({kFAS)mwmb{hke23kFi01}6w8>&=oPtiL)m zAVT_z=J79SbZ=rT=^!=%Bw;+^>AizkczT z@ZXT1*wi#vf%Fq6qLNO(#{YdrktALi1&6bVUHg28K$K5hPTh^WpjtlgXAHurnQzSb zB7Vhn+PEiR1PriUKdF1$!Mqeck4IS4)Db$fg@s#f!pZNqvT*P2F_&8WBf>K?RK^10 z_qHb{4;Es0*SW|{$9X1xB-dGPM_8+=|2lUkhVwsRb*}g@@Ejg{lfHq0GsR}v3xpM? ztmy*88$2esQ8NN7pY1O>aC8{DqKXq9aTci{^Kq%|=rFjf7uLx%BjBHf4Sq1h4FAXy z_h9TX>iTT5J3qt>%7$8?wt+#n6)tL@oM+PBp0(7?LpIcB$!_d4q1+>v4(;m+X3!w= ze!LR(4^+gqybNaG%U?OR3q5|cvDnyMMPhXKhC1Ro+EuIgrf`hJ@ODFU{Dj93FTgmJ z*Y2o5S6XV|xUT_-9tNE&E!qS+BeXj*Favy-v%5UeqJs0^cDTaijnE>G&J?no2}!$y zi4!G_AbOfUS-OQqyCqYbe@6LxmpoK@e<(CuX^t)1%A#A%(P2F(*Wyg~>4Qu74j<|cO_c5O#1Ea?5K}Az}+^=-|zh$;eFwSL+BrAOpg+7tbomk6{vNAFIZ(Ed{?MVzu6Lh_L}qJ0;d zS*eu5P5Mw&(@9}Xzlc6Ua`Oq;5ku)hGhUEtHRZ51ywYFWW6-SpphYPvo||6IHj1Gj?OOjQ z2I-nlh*zt*Sicw@Z`Cfh-EieV zYx&FVW3fH=hG9Y#E6bDHU@FA!ai(|JrF(8#?^GIwVS9-xXS!z$yoHl*h%uJ1cy1!R zpJ9xa_TqkJq}FNJI1VE_nJ<3B*v9tA%f<7QZJ!LnF;Kl$S4X}Pb>g%2}>l(Q5n|DGL1_cWuwN~=Q7-`0`VSQJ0wVTP5!=$eRG z3LH0|*A<#gsMheuwyrHBq**I3qr+_7>|=ZPx0B}tefcF z^K(O&287`afyEei?eFBS;5+y*-{D!X?~-%1R~6st0pb(reuE|$NP5C>e@cNhhLm)? zt%CfPQx(#_D^L;F`7(K79hfDHcinuji%~5Bj)#YMzDP_Ed9$5S-%pG&ezdudeq#QwAK`xSi?TgA^yv vRu6i~P?At@Jk3zS%g09@hj!H={}8{>hF_YHr^2eU|6K(a!k2|ydSdiCYZO(mrl?VAEv_ zU~(v5SbKB!ekOBZ$ckf9g^)s;#C&2J8b@TNN4vyb_X+Q>&w!@oDS@lwgRGr<;n7c_O0f}V-Ea1YH!#R<{6W|(f&%Zj8jbd G!u Date: Mon, 24 May 2021 14:43:32 +0300 Subject: [PATCH 2/2] Refactoring --- csvfile.hpp | 44 ++- drmsimulationmodel.cpp | 486 +++++++++++++++++++----------- drmsimulationmodel.hpp | 43 +-- edgemesh.cpp | 3 +- reducedmodeloptimizer_structs.hpp | 8 +- simulation_structs.hpp | 42 ++- simulationhistoryplotter.hpp | 105 ++++--- simulationmesh.hpp | 5 +- topologyenumerator.cpp | 2 +- trianglepatterngeometry.cpp | 4 +- trianglepatterngeometry.hpp | 2 +- utilities.hpp | 34 ++- 12 files changed, 511 insertions(+), 267 deletions(-) diff --git a/csvfile.hpp b/csvfile.hpp index de28e8c..bf17814 100755 --- a/csvfile.hpp +++ b/csvfile.hpp @@ -1,10 +1,13 @@ #ifndef CSVFILE_HPP #define CSVFILE_HPP +#include +#include #include #include #include #include #include +#include class csvFile; @@ -60,7 +63,46 @@ public: csvFile &operator<<(const std::string &val) { return write(escape(val)); } - template csvFile &operator<<(const T &val) { return write(val); } + template + csvFile &operator<<(const T &val) + { + return write(val); + } + + template + static std::vector> parse(const std::filesystem::path &csvFilepath) + { + std::vector> resultCSV; + if (!std::filesystem::exists(csvFilepath)) { + std::cerr << "The file does not exist:" << csvFilepath.string() << std::endl; + return resultCSV; + } + + std::ifstream inputfile(csvFilepath.string().c_str()); + if (!inputfile.is_open()) { + std::cerr << "Can't open file:" << csvFilepath.string() << std::endl; + return resultCSV; + } + std::vector row; + std::string line; + using Tokenizer = boost::tokenizer>; + while (std::getline(inputfile, line)) { + Tokenizer tokenizer(line); + row.resize(std::distance(tokenizer.begin(), tokenizer.end())); + std::transform(tokenizer.begin(), tokenizer.end(), row.begin(), [](const std::string &el) { + return boost::lexical_cast(el); + }); + // std::cout << std::endl; + // row.assign(tokenizer.begin(), tokenizer.end()); + // for (const auto &el : row) { + // std::cout << el << " "; + // } + // std::cout << std::endl; + resultCSV.push_back(row); + } + + return resultCSV; + } private: template csvFile &write(const T &val) { diff --git a/drmsimulationmodel.cpp b/drmsimulationmodel.cpp index 33b40c5..ddeadb6 100755 --- a/drmsimulationmodel.cpp +++ b/drmsimulationmodel.cpp @@ -22,7 +22,7 @@ void DRMSimulationModel::runUnitTests() // mesh.createSpanGrid(spanGridSize); beam.load(std::filesystem::path(groundOfTruthFolder).append("simpleBeam.ply").string()); std::unordered_map> fixedVertices; - fixedVertices[0] = std::unordered_set{0, 1, 2, 3}; + fixedVertices[0] = std::unordered_set{0, 1, 2}; fixedVertices[beam.VN() - 1] = std::unordered_set{1, 2}; std::unordered_map nodalForces{ {beam.VN() / 2, Vector6d({0, 1000, 1000, 0, 0, 0})}}; @@ -45,9 +45,10 @@ void DRMSimulationModel::runUnitTests() settings.Dtini = 0.1; settings.xi = 0.9969; settings.maxDRMIterations = 0; - totalResidualForcesNormThreshold = 1; - // settings.shouldDraw = true; - settings.beVerbose = true; + formFinder.totalResidualForcesNormThreshold = 1; + settings.shouldDraw = false; + settings.beVerbose = false; + // settings.shouldCreatePlots = true; SimulationResults simpleBeam_simulationResults = formFinder.executeSimulation(std::make_shared(beamSimulationJob), settings); simpleBeam_simulationResults.save(); @@ -58,10 +59,19 @@ void DRMSimulationModel::runUnitTests() assert(std::filesystem::exists(std::filesystem::path(simpleBeamGroundOfTruthBinaryFilename))); Eigen::MatrixXd simpleBeam_groundOfTruthDisplacements; Eigen::read_binary(simpleBeamGroundOfTruthBinaryFilename, simpleBeam_groundOfTruthDisplacements); - if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements)) { - std::cerr << "Error:Beam simulation produces wrong results." << std::endl; + + double error; + if (!simpleBeam_simulationResults.isEqual(simpleBeam_groundOfTruthDisplacements, error)) { + std::cerr << "Error:Beam simulation produces wrong results. Error:" << error << std::endl; // return; } +#ifdef POLYSCOPE_DEFINED + beam.registerForDrawing(); + simpleBeam_simulationResults.registerForDrawing(); + polyscope::show(); + beam.unregister(); + simpleBeam_simulationResults.unregister(); +#endif // Second example of the paper VCGEdgeMesh shortSpanGrid; @@ -109,13 +119,21 @@ void DRMSimulationModel::runUnitTests() // shortSpanGridshellSimulationResults.registerForDrawing( // shortSpanGridshellSimulationJob); // polyscope::show(); - assert( - shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements)); - if (!shortSpanGridshellSimulationResults.isEqual( - shortSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Short span simulation produces wrong results." << std::endl; - return; + assert(shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements, + error)); + if (!shortSpanGridshellSimulationResults.isEqual(shortSpanGridshell_groundOfTruthDisplacements, + error)) { + std::cerr << "Error:Short span simulation produces wrong results. Error:" << error + << std::endl; + // return; } +#ifdef POLYSCOPE_DEFINED + shortSpanGrid.registerForDrawing(); + shortSpanGridshellSimulationResults.registerForDrawing(); + polyscope::show(); + shortSpanGrid.unregister(); + shortSpanGridshellSimulationResults.unregister(); +#endif // Third example of the paper VCGEdgeMesh longSpanGrid; @@ -163,6 +181,7 @@ void DRMSimulationModel::runUnitTests() << std::endl; } longSpanGridshell_simulationJob.pMesh->setBeamCrossSection(CrossSectionType{0.03, 0.026}); + SimulationResults longSpanGridshell_simulationResults = formFinder.executeSimulation(std::make_shared( longSpanGridshell_simulationJob), @@ -177,12 +196,21 @@ void DRMSimulationModel::runUnitTests() Eigen::MatrixXd longSpanGridshell_groundOfTruthDisplacements; Eigen::read_binary(longSpan_groundOfTruthBinaryFilename, longSpanGridshell_groundOfTruthDisplacements); - assert( - longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements)); - if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements)) { - std::cerr << "Error:Long span simulation produces wrong results." << std::endl; - return; + assert(longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements, + error)); + if (!longSpanGridshell_simulationResults.isEqual(longSpanGridshell_groundOfTruthDisplacements, + error)) { + std::cerr << "Error:Long span simulation produces wrong results. Error:" << error + << std::endl; + // return; } +#ifdef POLYSCOPE_DEFINED + longSpanGrid.registerForDrawing(); + longSpanGridshell_simulationResults.registerForDrawing(); + polyscope::show(); + longSpanGrid.unregister(); + longSpanGridshell_simulationResults.unregister(); +#endif std::cout << "Form finder unit tests succesufully passed." << std::endl; @@ -1130,14 +1158,11 @@ void DRMSimulationModel::updateNodalMasses() const double SkTranslational = elem.material.youngsModulus * elem.A / elem.length; translationalSumSk += SkTranslational; const double lengthToThe3 = std::pow(elem.length, 3); - const long double SkRotational_I2 - = elem.material.youngsModulus * elem.I2 - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const long double SkRotational_I3 - = elem.material.youngsModulus * elem.I3 - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia - const long double SkRotational_J = elem.material.youngsModulus * elem.J - / lengthToThe3; // TODO: I2->t2,I3->t3,t1->polar inertia + const long double SkRotational_I2 = elem.material.youngsModulus * elem.I2 + / lengthToThe3; + const long double SkRotational_I3 = elem.material.youngsModulus * elem.I3 + / lengthToThe3; + const long double SkRotational_J = elem.material.youngsModulus * elem.J / lengthToThe3; rotationalSumSk_I2 += SkRotational_I2; rotationalSumSk_I3 += SkRotational_I3; rotationalSumSk_J += SkRotational_J; @@ -1336,9 +1361,8 @@ void DRMSimulationModel::updateElementalFrames() } void DRMSimulationModel::applyForcedDisplacements( - const std::unordered_map nodalForcedDisplacements) + const std::unordered_map &nodalForcedDisplacements) { - const int gradualDisplacementSteps = 500; for (const std::pair vertexIndexDisplacementPair : nodalForcedDisplacements) { const VertexIndex vi = vertexIndexDisplacementPair.first; @@ -1347,9 +1371,9 @@ void DRMSimulationModel::applyForcedDisplacements( VectorType displacementVector(vertexDisplacement(0), vertexDisplacement(1), vertexDisplacement(2)); - if (mCurrentSimulationStep < gradualDisplacementSteps) { + if (mCurrentSimulationStep < mSettings.gradualForcedDisplacementSteps) { displacementVector *= mCurrentSimulationStep - / static_cast(gradualDisplacementSteps); + / static_cast(mSettings.gradualForcedDisplacementSteps); } pMesh->vert[vi].P() = node.initialLocation + displacementVector; node.displacements = Vector6d({displacementVector[0], @@ -1553,6 +1577,13 @@ void DRMSimulationModel::printCurrentState() const std::cout << "Simulation steps executed:" << mCurrentSimulationStep << std::endl; std::cout << "Residual forces norm: " << pMesh->totalResidualForcesNorm << std::endl; std::cout << "Kinetic energy:" << pMesh->currentTotalKineticEnergy << std::endl; + static std::chrono::steady_clock::time_point begin = std::chrono::steady_clock::now(); + const auto timePerNodePerIteration = std::chrono::duration_cast( + std::chrono::steady_clock::now() - begin) + .count() + * 1e-6 / (mCurrentSimulationStep * pMesh->VN()); + std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" << std::endl; + std::cout << "time(s)/(iterations*node) = " << timePerNodePerIteration << std::endl; } void DRMSimulationModel::printDebugInfo() const @@ -1694,22 +1725,26 @@ void DRMSimulationModel::draw(const std::string &screenshotsFolder) ->addEdgeScalarQuantity("I3", I3); // Specify the callback - PolyscopeInterface::addUserCallback([&]() { - // Since options::openImGuiWindowForUserCallback == true by default, - // we can immediately start using ImGui commands to build a UI + static bool calledOnce = false; + if (!calledOnce) { + PolyscopeInterface::addUserCallback([&]() { + // Since options::openImGuiWindowForUserCallback == true by default, + // we can immediately start using ImGui commands to build a UI - ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, - // instead of full width. Must have - // matching PopItemWidth() below. + ImGui::PushItemWidth(100); // Make ui elements 100 pixels wide, + // instead of full width. Must have + // matching PopItemWidth() below. - ImGui::InputInt("Simulation drawing step", - &mSettings.drawingStep); // set a int variable - ImGui::Checkbox("Enable drawing", - &mSettings.shouldDraw); // set a int variable - ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); + ImGui::InputInt("Simulation drawing step", + &mSettings.drawingStep); // set a int variable + ImGui::Checkbox("Enable drawing", + &mSettings.shouldDraw); // set a int variable + ImGui::Text("Number of simulation steps: %zu", mCurrentSimulationStep); - ImGui::PopItemWidth(); - }); + ImGui::PopItemWidth(); + }); + calledOnce = true; + } if (!screenshotsFolder.empty()) { static bool firstDraw = true; @@ -1742,16 +1777,27 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptr &externalForce : pJob->nodalExternalForces) { totalExternalForcesNorm += externalForce.second.norm(); } - totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination - * totalExternalForcesNorm; - // if (!pJob->nodalExternalForces.empty()) { - // double externalForcesNorm = 0; - // for (const auto &externalForce : pJob->nodalExternalForces) { - // externalForcesNorm += externalForce.second.norm(); - // } - // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; - // } + if (!pJob->nodalExternalForces.empty()) { + totalResidualForcesNormThreshold = settings.totalExternalForcesNormPercentageTermination + * totalExternalForcesNorm; + } else { + totalResidualForcesNormThreshold = 1; + std::cout << "No forces setted default residual forces norm threshold" << std::endl; + } + + if (mSettings.beVerbose) { + std::cout << "totalResidualForcesNormThreshold:" << totalResidualForcesNormThreshold + << std::endl; + } + + // if (!pJob->nodalExternalForces.empty()) { + // double externalForcesNorm = 0; + // for (const auto &externalForce : pJob->nodalExternalForces) { + // externalForcesNorm += externalForce.second.norm(); + // } + // mSettings.totalResidualForcesNormThreshold = externalForcesNorm * 1e-2; + // } constrainedVertices = pJob->constrainedVertices; if (!pJob->nodalForcedDisplacements.empty()) { @@ -1773,18 +1819,20 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrVN() << " nodes and " << pMesh->EN() << " elements." << std::endl; } + vcg::tri::UpdateBounding::Box(*pMesh); computeRigidSupports(); for (auto fixedVertex : pJob->constrainedVertices) { assert(fixedVertex.first < pMesh->VN()); } #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw ) { - PolyscopeInterface::init(); - polyscope::registerCurveNetwork( - meshPolyscopeLabel, pMesh->getEigenVertices(), pMesh->getEigenEdges()); - // registerWorldAxes(); - } + if (mSettings.shouldDraw || mSettings.isDebugMode) { + PolyscopeInterface::init(); + polyscope::registerCurveNetwork(meshPolyscopeLabel, + pMesh->getEigenVertices(), + pMesh->getEigenEdges()); + // registerWorldAxes(); + } #endif updateElementalFrames(); @@ -1793,130 +1841,203 @@ SimulationResults DRMSimulationModel::executeSimulation(const std::shared_ptrnodalExternalForces, constrainedVertices); - while (settings.maxDRMIterations == 0 || - mCurrentSimulationStep < settings.maxDRMIterations) { - // while (true) { - updateNormalDerivatives(); - updateT1Derivatives(); - updateRDerivatives(); - updateT2Derivatives(); - updateT3Derivatives(); - updateResidualForcesOnTheFly(constrainedVertices); + const size_t movingAverageSampleSize = 50; + std::queue residualForcesMovingAverageHistorySample; - // TODO: write parallel function for updating positions - // TODO: Make a single function out of updateResidualForcesOnTheFly - // updatePositionsOnTheFly - // updatePositionsOnTheFly(constrainedVertices); - updateNodalMasses(); - updateNodalAccelerations(); - updateNodalVelocities(); - updateNodalDisplacements(); - applyDisplacements(constrainedVertices); - if (!pJob->nodalForcedDisplacements.empty()) { - applyForcedDisplacements( + double residualForcesMovingAverageDerivativeMax = 0; + while (settings.maxDRMIterations == 0 || mCurrentSimulationStep < settings.maxDRMIterations) { + // while (true) { + updateNormalDerivatives(); + updateT1Derivatives(); + updateRDerivatives(); + updateT2Derivatives(); + updateT3Derivatives(); + updateResidualForcesOnTheFly(constrainedVertices); - pJob->nodalForcedDisplacements); - } - // if (!pJob->nodalForcedNormals.empty()) { - // applyForcedNormals(pJob->nodalForcedNormals); - // } - updateElementalLengths(); - // pMesh->previousTotalPotentialEnergykN = - // pMesh->currentTotalPotentialEnergykN; - // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + // TODO: write parallel function for updating positions + // TODO: Make a single function out of updateResidualForcesOnTheFly + // updatePositionsOnTheFly + // updatePositionsOnTheFly(constrainedVertices); + updateNodalMasses(); + updateNodalAccelerations(); + updateNodalVelocities(); + updateNodalDisplacements(); + applyDisplacements(constrainedVertices); + if (!pJob->nodalForcedDisplacements.empty()) { + applyForcedDisplacements( - // if (mCurrentSimulationStep != 0) { - // history.stepPulse(*pMesh); - // } - - if (std::isnan(pMesh->currentTotalKineticEnergy)) { - if (mSettings.beVerbose) { - std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; + pJob->nodalForcedDisplacements); + } + // if (!pJob->nodalForcedNormals.empty()) { + // applyForcedNormals(pJob->nodalForcedNormals); + // } + updateElementalLengths(); + mCurrentSimulationStep++; + + double sumOfDisplacementsNorm = 0; + for (size_t vi = 0; vi < pMesh->VN(); vi++) { + sumOfDisplacementsNorm += pMesh->nodes[vi].displacements.getTranslation().norm(); + } + sumOfDisplacementsNorm /= pMesh->bbox.Diag(); + pMesh->sumOfNormalizedDisplacementNorms = sumOfDisplacementsNorm; + // pMesh->residualForcesMovingAverage = (pMesh->totalResidualForcesNorm + // + mCurrentSimulationStep + // * pMesh->residualForcesMovingAverage) + // / (1 + mCurrentSimulationStep); + pMesh->residualForcesMovingAverage += pMesh->totalResidualForcesNorm + / movingAverageSampleSize; + residualForcesMovingAverageHistorySample.push(pMesh->residualForcesMovingAverage); + if (movingAverageSampleSize < residualForcesMovingAverageHistorySample.size()) { + const double firstElementValue = residualForcesMovingAverageHistorySample.front(); + pMesh->residualForcesMovingAverage -= firstElementValue / movingAverageSampleSize; + residualForcesMovingAverageHistorySample.pop(); + + pMesh->residualForcesMovingAverageDerivativeNorm + = std::abs((residualForcesMovingAverageHistorySample.back() + - residualForcesMovingAverageHistorySample.front())) + / (movingAverageSampleSize); + residualForcesMovingAverageDerivativeMax + = std::max(pMesh->residualForcesMovingAverageDerivativeNorm, + residualForcesMovingAverageDerivativeMax); + pMesh->residualForcesMovingAverageDerivativeNorm + /= residualForcesMovingAverageDerivativeMax; + // std::cout << "Normalized derivative:" + // << residualForcesMovingAverageDerivativeNorm + // << std::endl; + } + + // pMesh->previousTotalPotentialEnergykN = + // pMesh->currentTotalPotentialEnergykN; + // pMesh->currentTotalPotentialEnergykN = computeTotalPotentialEnergy() / 1000; + // timePerNodePerIterationHistor.push_back(timePerNodePerIteration); + if (mSettings.beVerbose && mSettings.isDebugMode + && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + printCurrentState(); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver", + // movingAverages); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // movingAveragesDerivatives); + // draw(); + std::cout << "Mov aver deriv:" << pMesh->residualForcesMovingAverageDerivativeNorm + << std::endl; + // SimulationResulnodalForcedDisplacementstsReporter::createPlot("Number of Steps", + // "Time/(#nodes*#iterations)", + // timePerNodePerIterationHistory); + } + + if ((mSettings.shouldCreatePlots || mSettings.isDebugMode) && mCurrentSimulationStep != 0) { + history.stepPulse(*pMesh); + } + + if (std::isnan(pMesh->currentTotalKineticEnergy)) { + if (mSettings.beVerbose) { + std::cout << "Infinite kinetic energy detected.Exiting.." << std::endl; + } + break; } - break; - } #ifdef POLYSCOPE_DEFINED - if (mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) /* && + if ((mSettings.shouldDraw && mCurrentSimulationStep % mSettings.drawingStep == 0) + || (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0)) /* && + currentSimulationStep > maxDRMIterations*/ - { - // std::string saveTo = std::filesystem::current_path() - // .append("Debugging_files") - // .append("Screenshots") - // .string(); - draw(); - // yValues = history.kineticEnergy; - // plotHandle = matplot::scatter(xPlot, yValues); - // label = "Log of Kinetic energy"; - // plotHandle->legend_string(label); + { + // std::string saveTo = std::filesystem::current_path() + // .append("Debugging_files") + // .append("Screenshots") + // .string(); + draw(); + // yValues = history.kineticEnergy; + // plotHandle = matplot::scatter(xPlot, yValues); + // label = "Log of Kinetic energy"; + // plotHandle->legend_string(label); - // shouldUseKineticEnergyThreshold = true; - } + // shouldUseKineticEnergyThreshold = true; + } #endif - if (mSettings.shouldCreatePlots && mCurrentSimulationStep % 10 == 0) { - printCurrentState(); - SimulationResultsReporter::createPlot( - "Number of Steps", "Log of Kinetic energy", history.kineticEnergy); - } - // t = t + Dt; - mCurrentSimulationStep++; - // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy - // << std::endl; - // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm - // << std::endl; - // Kinetic energy convergence - if ((mSettings.shouldUseTranslationalKineticEnergyThreshold || mCurrentSimulationStep > 500000) - && pMesh->currentTotalTranslationalKineticEnergy - < mSettings.totalTranslationalKineticEnergyThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" - << std::endl; - } -#ifdef POLYSCOPE_DEFINED - std::cout << "Warning: The kinetic energy of the system was " - " used as a convergence criterion" - << std::endl; -#endif - break; - } - // Residual forces norm convergence - if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy - /*|| + if (mSettings.isDebugMode && mCurrentSimulationStep % mSettings.debugModeStep == 0) { + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces mov aver deriv", + // movingAveragesDerivatives_norm); + SimulationResultsReporter::createPlot("Number of Steps", + "Residual Forces mov aver", + history.residualForcesMovingAverage); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Residual Forces", + // history.residualForces); + // SimulationResultsReporter::createPlot("Number of Steps", + // "Log of Kinetic energy", + // history.kineticEnergy); + } + // t = t + Dt; + // std::cout << "Kinetic energy:" << mesh.currentTotalKineticEnergy + // << std::endl; + // std::cout << "Residual forces norm:" << mesh.totalResidualForcesNorm + // << std::endl; + // Residual forces norm convergence + if (pMesh->previousTotalKineticEnergy > pMesh->currentTotalKineticEnergy + && mCurrentSimulationStep > movingAverageSampleSize + && (pJob->nodalForcedDisplacements.empty() + || mCurrentSimulationStep > mSettings.gradualForcedDisplacementSteps) + /*|| mesh->previousTotalPotentialEnergykN > mesh->currentTotalPotentialEnergykN*/ - /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { - if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { - if (mSettings.beVerbose) { - std::cout << "Simulation converged." << std::endl; - printCurrentState(); - std::cout << "Total potential:" << pMesh->currentTotalPotentialEnergykN << " kNm" - << std::endl; - } - break; - // } - } - // printCurrentState(); - // for (VertexType &v : mesh->vert) { - // Node &node = mesh->nodes[v]; - // node.displacements = node.displacements - node.velocity * Dt; - // } - // applyDisplacements(constrainedVertices); - // if (!pJob->nodalForcedDisplacements.empty()) { - // applyForcedDisplacements( + /*|| mesh->currentTotalPotentialEnergykN < minPotentialEnergy*/) { + // if (pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold) { + const bool fullfillsKineticEnergyTerminationCriterion + = mSettings.shouldUseTranslationalKineticEnergyThreshold + && pMesh->currentTotalTranslationalKineticEnergy + < mSettings.totalTranslationalKineticEnergyThreshold + && mCurrentSimulationStep > 20 && numOfDampings > 2; + const bool fullfillsResidualForcesNormTerminationCriterion + = pMesh->totalResidualForcesNorm < totalResidualForcesNormThreshold; + const bool fullfillsMovingAverageNormTerminationCriterion + = pMesh->residualForcesMovingAverage + < mSettings.residualForcesMovingAverageNormThreshold; + /*pMesh->residualForcesMovingAverage < totalResidualForcesNormThreshold*/ + const bool shouldTerminate + = fullfillsKineticEnergyTerminationCriterion + // || fullfillsMovingAverageNormTerminationCriterion + || fullfillsResidualForcesNormTerminationCriterion; + if (shouldTerminate) { + if (mSettings.beVerbose) { + std::cout << "Simulation converged." << std::endl; + printCurrentState(); + if (fullfillsResidualForcesNormTerminationCriterion) { + std::cout << "Converged using residual forces norm threshold criterion" + << std::endl; + } else if (fullfillsKineticEnergyTerminationCriterion) { + std::cout << "Warning: The kinetic energy of the system was " + " used as a convergence criterion" + << std::endl; + } else if (fullfillsMovingAverageNormTerminationCriterion) { + std::cout << "Converged using residual forces moving average derivative norm " + "threshold criterion" + << std::endl; + } + } + break; + // } + } + // printCurrentState(); + // for (VertexType &v : mesh->vert) { + // Node &node = mesh->nodes[v]; + // node.displacements = node.displacements - node.velocity * Dt; + // } + // applyDisplacements(constrainedVertices); + // if (!pJob->nodalForcedDisplacements.empty()) { + // applyForcedDisplacements( - // pJob->nodalForcedDisplacements); - // } - // updateElementalLengths(); - resetVelocities(); - Dt = Dt * mSettings.xi; - ++numOfDampings; - } - - if (mSettings.debugMessages) { - printDebugInfo(); - } + // pJob->nodalForcedDisplacements); + // } + // updateElementalLengths(); + resetVelocities(); + Dt = Dt * mSettings.xi; + ++numOfDampings; + // std::cout << "Number of dampings:" << numOfDampings << endl; + } } SimulationResults results = computeResults(pJob); @@ -1929,20 +2050,21 @@ mesh->currentTotalPotentialEnergykN*/ << std::endl; results.converged = false; } else if (std::isnan(pMesh->currentTotalKineticEnergy)) { - results.converged = false; + std::cerr << "Simulation did not converge due to infinite kinetic energy." << std::endl; + results.converged = false; } else if (mSettings.beVerbose) { - auto t2 = std::chrono::high_resolution_clock::now(); - double simulationDuration = - std::chrono::duration_cast(t2 - t1).count(); - simulationDuration /= 1000; - std::cout << "Simulation converged after " << simulationDuration << "s" - << std::endl; - std::cout << "Time/(node*iteration) " - << simulationDuration / - (static_cast(mCurrentSimulationStep) * pMesh->VN()) - << "s" << std::endl; - std::cout << "Number of dampings:" << numOfDampings << endl; + auto t2 = std::chrono::high_resolution_clock::now(); + double simulationDuration = std::chrono::duration_cast(t2 - t1) + .count(); + simulationDuration /= 1000; + std::cout << "Simulation converged after " << simulationDuration << "s" << std::endl; + std::cout << "Time/(node*iteration) " + << simulationDuration / (static_cast(mCurrentSimulationStep) * pMesh->VN()) + << "s" << std::endl; + std::cout << "Number of dampings:" << numOfDampings << endl; + std::cout << "Number of steps:" << mCurrentSimulationStep << endl; + results.converged = true; } // mesh.printVertexCoordinates(mesh.VN() / 2); #ifdef POLYSCOPE_DEFINED @@ -1994,7 +2116,7 @@ mesh->currentTotalPotentialEnergykN*/ } } - if (mSettings.shouldCreatePlots) { + if (mSettings.shouldCreatePlots && results.converged) { SimulationResultsReporter reporter; reporter.reportResults({results}, "Results", pJob->pMesh->getLabel()); } diff --git a/drmsimulationmodel.hpp b/drmsimulationmodel.hpp index 928ef08..4e790fc 100755 --- a/drmsimulationmodel.hpp +++ b/drmsimulationmodel.hpp @@ -23,28 +23,33 @@ struct DifferentiateWithRespectTo { class DRMSimulationModel { public: - struct Settings { - bool debugMessages{false}; - bool shouldDraw{false}; - bool beVerbose{false}; - bool shouldCreatePlots{false}; - int drawingStep{1}; - double totalTranslationalKineticEnergyThreshold{1e-10}; - double totalExternalForcesNormPercentageTermination{1e-3}; - double Dtini{0.1}; - double xi{0.9969}; - int maxDRMIterations{0}; - bool shouldUseTranslationalKineticEnergyThreshold{false}; - Settings() {} - }; + struct Settings + { + bool isDebugMode{false}; + int debugModeStep{10000}; + bool shouldDraw{false}; + bool beVerbose{false}; + bool shouldCreatePlots{false}; + int drawingStep{1}; + double totalTranslationalKineticEnergyThreshold{1e-10}; + double totalExternalForcesNormPercentageTermination{1e-3}; + double residualForcesMovingAverageDerivativeNormThreshold{1e-8}; + double residualForcesMovingAverageNormThreshold{1e-8}; + double Dtini{0.1}; + double xi{0.9969}; + int maxDRMIterations{0}; + bool shouldUseTranslationalKineticEnergyThreshold{false}; + int gradualForcedDisplacementSteps{100}; + Settings() {} + }; + double totalResidualForcesNormThreshold{1}; -private: + private: Settings mSettings; double Dt{mSettings.Dtini}; bool checkedForMaximumMoment{false}; double externalMomentsNorm{0}; size_t mCurrentSimulationStep{0}; - double totalResidualForcesNormThreshold{1e-3}; matplot::line_handle plotHandle; std::vector plotYValues; size_t numOfDampings{0}; @@ -80,8 +85,7 @@ private: void updateNodalDisplacements(); void applyForcedDisplacements( - const std::unordered_map - nodalForcedDisplacements); + const std::unordered_map &nodalForcedDisplacements); Vector6d computeElementTorsionalForce( const Element &element, const Vector6d &displacementDifference, @@ -208,8 +212,7 @@ private: executeSimulation(const std::shared_ptr &pJob, const Settings &settings = Settings()); - void runUnitTests(); - void setTotalResidualForcesNormThreshold(double value); + static void runUnitTests(); }; template PointType Cross(PointType p1, PointType p2) { diff --git a/edgemesh.cpp b/edgemesh.cpp index d1f4f2d..0fe03e1 100755 --- a/edgemesh.cpp +++ b/edgemesh.cpp @@ -200,7 +200,6 @@ bool VCGEdgeMesh::load(const string &plyFilename) { getEdges(eigenEdges); getVertices(eigenVertices); vcg::tri::UpdateTopology::VertexEdge(*this); - label = std::filesystem::path(plyFilename).stem().string(); const bool printDebugInfo = false; if (printDebugInfo) { @@ -208,7 +207,7 @@ bool VCGEdgeMesh::load(const string &plyFilename) { std::cout << "Mesh has " << EN() << " edges." << std::endl; } - label=std::filesystem::path(plyFilename).stem(); + label = std::filesystem::path(plyFilename).stem().string(); return true; } diff --git a/reducedmodeloptimizer_structs.hpp b/reducedmodeloptimizer_structs.hpp index fa881f5..7c3250a 100644 --- a/reducedmodeloptimizer_structs.hpp +++ b/reducedmodeloptimizer_structs.hpp @@ -126,12 +126,14 @@ struct Colors } xRange x; x.fromString(json.at(jsonXRangeKey)); + xRanges.push_back(x); xRangeIndex++; } numberOfFunctionCalls = json.at(JsonKeys::NumberOfFunctionCalls); solverAccuracy = json.at(JsonKeys::SolverAccuracy); objectiveWeights.translational = json.at(JsonKeys::TranslationalObjectiveWeight); + objectiveWeights.rotational = 2 - objectiveWeights.translational; return true; } @@ -197,6 +199,7 @@ struct Colors inline void updateMeshWithOptimalVariables(const std::vector &x, SimulationMesh &m) { + assert(CrossSectionType::name == RectangularBeamDimensions::name); const double E = x[0]; const double A = x[1]; const double beamWidth = std::sqrt(A); @@ -207,7 +210,7 @@ struct Colors const double I3 = x[4]; for (EdgeIndex ei = 0; ei < m.EN(); ei++) { Element &e = m.elements[ei]; - e.setDimensions(RectangularBeamDimensions(beamWidth, beamHeight)); + e.setDimensions(CrossSectionType(beamWidth, beamHeight)); e.setMaterial(ElementMaterial(e.material.poissonsRatio, E)); e.J = J; e.I2 = I2; @@ -494,6 +497,7 @@ struct Colors const ReducedPatternOptimization::Results &reducedPattern_optimizationResults, const shared_ptr &pTiledReducedPattern_simulationMesh) { + assert(CrossSectionType::name == RectangularBeamDimensions::name); // Set reduced parameters fron the optimization results std::unordered_map optimalXVariables(reducedPattern_optimizationResults.optimalXNameValuePairs.begin(), @@ -504,7 +508,7 @@ struct Colors const double A = optimalXVariables.at(ALabel); const double beamWidth = std::sqrt(A); const double beamHeight = beamWidth; - RectangularBeamDimensions elementDimensions(beamWidth, beamHeight); + CrossSectionType elementDimensions(beamWidth, beamHeight); const double poissonsRatio = 0.3; const std::string ymLabel = "E"; diff --git a/simulation_structs.hpp b/simulation_structs.hpp index fd75566..9d60766 100755 --- a/simulation_structs.hpp +++ b/simulation_structs.hpp @@ -9,11 +9,14 @@ struct SimulationHistory { size_t numberOfSteps{0}; std::string label; - std::vector residualForces; + std::vector residualForcesLog; std::vector kineticEnergy; std::vector potentialEnergies; std::vector redMarks; std::vector greenMarks; + std::vector residualForcesMovingAverage; + std::vector sumOfNormalizedDisplacementNorms; + // std::vector residualForcesMovingAverageDerivativesLog; void markRed(const size_t &stepNumber) { redMarks.push_back(stepNumber); } @@ -22,13 +25,21 @@ struct SimulationHistory { void stepPulse(const SimulationMesh &mesh) { kineticEnergy.push_back(log(mesh.currentTotalKineticEnergy)); // potentialEnergy.push_back(mesh.totalPotentialEnergykN); - residualForces.push_back(mesh.totalResidualForcesNorm); + residualForcesLog.push_back(std::log(mesh.totalResidualForcesNorm)); + residualForcesMovingAverage.push_back(mesh.residualForcesMovingAverage); + sumOfNormalizedDisplacementNorms.push_back(mesh.sumOfNormalizedDisplacementNorms); + // residualForcesMovingAverageDerivativesLog.push_back( + // std::log(mesh.residualForcesMovingAverageDerivativeNorm)); } - void clear() { - residualForces.clear(); - kineticEnergy.clear(); - potentialEnergies.clear(); + void clear() + { + residualForcesLog.clear(); + kineticEnergy.clear(); + potentialEnergies.clear(); + residualForcesMovingAverage.clear(); + sumOfNormalizedDisplacementNorms.clear(); + // residualForcesMovingAverageDerivativesLog.clear(); } }; @@ -399,6 +410,9 @@ public: } void unregister(const std::string &meshLabel) { + if (polyscope::getCurveNetwork(meshLabel) == nullptr) { + return; + } polyscope::getCurveNetwork(meshLabel)->removeQuantity("External force_" + label); polyscope::getCurveNetwork(meshLabel)->removeQuantity("Boundary conditions_" + label); } @@ -478,23 +492,29 @@ struct SimulationResults void save(const std::string &outputFolder = std::string()) { - const std::filesystem::path outputPath(outputFolder); - job->save(std::filesystem::path(outputPath).append("SimulationJob").string()); + const std::filesystem::path outputFolderPath = outputFolder.empty() + ? std::filesystem::current_path() + : std::filesystem::path(outputFolder); + std::filesystem::path simulationJobOutputFolderPath + = std::filesystem::path(outputFolderPath).append("SimulationJob"); + std::filesystem::create_directories(simulationJobOutputFolderPath); + job->save(simulationJobOutputFolderPath); const std::string filename(getLabel() + "_displacements.eigenBin"); Eigen::MatrixXd m = Utilities::toEigenMatrix(displacements); - Eigen::write_binary(std::filesystem::path(outputPath).append(filename).string(), m); + Eigen::write_binary(std::filesystem::path(outputFolderPath).append(filename).string(), m); - saveDeformedModel(outputFolder); + saveDeformedModel(outputFolderPath); } // The comparison of the results happens comparing the 6-dof nodal // displacements - bool isEqual(const Eigen::MatrixXd &nodalDisplacements) + bool isEqual(const Eigen::MatrixXd &nodalDisplacements, double &error) { assert(nodalDisplacements.cols() == 6); Eigen::MatrixXd eigenDisplacements = Utilities::toEigenMatrix(this->displacements); const double errorNorm = (eigenDisplacements - nodalDisplacements).norm(); + error = errorNorm; return errorNorm < 1e-10; // return eigenDisplacements.isApprox(nodalDisplacements); } diff --git a/simulationhistoryplotter.hpp b/simulationhistoryplotter.hpp index 9beb5fd..422a4a8 100755 --- a/simulationhistoryplotter.hpp +++ b/simulationhistoryplotter.hpp @@ -38,11 +38,11 @@ struct SimulationResultsReporter { file << "\n"; } - if (!history.residualForces.empty()) { + if (!history.residualForcesLog.empty()) { file << "Residual forces" << "\n"; for (size_t step = 0; step < numberOfSteps; step++) { - file << history.residualForces[step] << "\n"; + file << history.residualForcesLog[step] << "\n"; } file << "\n"; } @@ -92,27 +92,31 @@ struct SimulationResultsReporter { } } - static void createPlot(const std::string &xLabel, const std::string &yLabel, + static void createPlot(const std::string &xLabel, + const std::string &yLabel, const std::vector &YvaluesToPlot, - const std::string &saveTo = {}) { - auto x = - matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); - createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo); - // matplot::figure(true); - // matplot::hold(matplot::on); + const std::string &saveTo = {}) + { + if (YvaluesToPlot.size() < 2) { + return; + } + auto x = matplot::linspace(0, YvaluesToPlot.size() - 1, YvaluesToPlot.size()); + createPlot(xLabel, yLabel, x, YvaluesToPlot, saveTo); + // matplot::figure(true); + // matplot::hold(matplot::on); - // ->marker_indices(history.redMarks) - // ->marker_indices(truncatedRedMarks) - // .marker_color("r") - // ->marker_size(1) - ; - // auto greenMarksY = matplot::transform( - // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; - // }); - // matplot::scatter(history.greenMarks, greenMarksY) - // ->color("green") - // .marker_size(10); - // matplot::hold(matplot::off); + // ->marker_indices(history.redMarks) + // ->marker_indices(truncatedRedMarks) + // .marker_color("r") + // ->marker_size(1) + ; + // auto greenMarksY = matplot::transform( + // history.greenMarks, [&](auto x) { return history.kineticEnergy[x]; + // }); + // matplot::scatter(history.greenMarks, greenMarksY) + // ->color("green") + // .marker_size(10); + // matplot::hold(matplot::off); } void createPlots(const SimulationHistory &history, @@ -124,27 +128,54 @@ struct SimulationResultsReporter { std::filesystem::create_directory(graphsFolder.string()); if (!history.kineticEnergy.empty()) { - createPlot("Number of Iterations", "Log of Kinetic Energy", - history.kineticEnergy, - std::filesystem::path(graphsFolder) - .append("KineticEnergy" + graphSuffix + ".png") - .string()); + createPlot("Number of Iterations", + "Log of Kinetic Energy log", + history.kineticEnergy, + std::filesystem::path(graphsFolder) + .append("KineticEnergyLog_" + graphSuffix + ".png") + .string()); } - if (!history.residualForces.empty()) { - createPlot("Number of Iterations", "Residual Forces norm", - history.residualForces, - std::filesystem::path(graphsFolder) - .append("ResidualForces" + graphSuffix + ".png") - .string()); + if (!history.residualForcesLog.empty()) { + createPlot("Number of Iterations", + "Residual Forces norm log", + history.residualForcesLog, + std::filesystem::path(graphsFolder) + .append("ResidualForcesLog_" + graphSuffix + ".png") + .string()); } if (!history.potentialEnergies.empty()) { - createPlot("Number of Iterations", "Potential energy", - history.potentialEnergies, - std::filesystem::path(graphsFolder) - .append("PotentialEnergy" + graphSuffix + ".png") - .string()); + createPlot("Number of Iterations", + "Potential energy", + history.potentialEnergies, + std::filesystem::path(graphsFolder) + .append("PotentialEnergy_" + graphSuffix + ".png") + .string()); + } + if (!history.residualForcesMovingAverage.empty()) { + createPlot("Number of Iterations", + "Residual forces moving average", + history.residualForcesMovingAverage, + std::filesystem::path(graphsFolder) + .append("ResidualForcesMovingAverage_" + graphSuffix + ".png") + .string()); + } + // if (!history.residualForcesMovingAverageDerivativesLog.empty()) { + // createPlot("Number of Iterations", + // "Residual forces moving average derivative log", + // history.residualForcesMovingAverageDerivativesLog, + // std::filesystem::path(graphsFolder) + // .append("ResidualForcesMovingAverageDerivativesLog_" + graphSuffix + ".png") + // .string()); + // } + if (!history.sumOfNormalizedDisplacementNorms.empty()) { + createPlot("Number of Iterations", + "Sum of normalized displacement norms", + history.sumOfNormalizedDisplacementNorms, + std::filesystem::path(graphsFolder) + .append("SumOfNormalizedDisplacementNorms_" + graphSuffix + ".png") + .string()); } } }; diff --git a/simulationmesh.hpp b/simulationmesh.hpp index 95e02b6..5dea7af 100755 --- a/simulationmesh.hpp +++ b/simulationmesh.hpp @@ -9,7 +9,7 @@ struct Element; struct Node; using CrossSectionType = RectangularBeamDimensions; -// using CrossSectionType = CylindricalBeamDimensions; +//using CrossSectionType = CylindricalBeamDimensions; class SimulationMesh : public VCGEdgeMesh { private: @@ -44,6 +44,9 @@ public: double totalResidualForcesNorm{0}; double currentTotalPotentialEnergykN{0}; double previousTotalPotentialEnergykN{0}; + double residualForcesMovingAverageDerivativeNorm{0}; + double residualForcesMovingAverage{0}; + double sumOfNormalizedDisplacementNorms{0}; bool save(const std::string &plyFilename = std::string()); void setBeamCrossSection(const CrossSectionType &beamDimensions); void setBeamMaterial(const double &pr, const double &ym); diff --git a/topologyenumerator.cpp b/topologyenumerator.cpp index 42f242f..9aa76d9 100755 --- a/topologyenumerator.cpp +++ b/topologyenumerator.cpp @@ -464,7 +464,7 @@ void TopologyEnumerator::computeValidPatterns( const bool tiledPatternHasEdgesWithAngleSmallerThanThreshold = patternGeometry.hasAngleSmallerThanThreshold(numberOfNodesPerSlot, 5); if (tiledPatternHasEdgesWithAngleSmallerThanThreshold) { - if (debugIsOn || savePlyFiles) { + if (debugIsOn /*|| savePlyFiles*/) { if (savePlyFiles) { auto danglingEdgesPath = std::filesystem::path(resultsPath) .append("ExceedingAngleThreshold"); diff --git a/trianglepatterngeometry.cpp b/trianglepatterngeometry.cpp index 9eb7e61..e1cbe47 100755 --- a/trianglepatterngeometry.cpp +++ b/trianglepatterngeometry.cpp @@ -435,7 +435,7 @@ bool PatternGeometry::hasAngleSmallerThanThreshold(const std::vector &nu thetaAnglesOfIncidentVectors[thetaAngleIndex] - thetaAnglesOfIncidentVectors[thetaAngleIndex - 1]); if (absAngleDifference < angleThreshold_rad - && absAngleDifference > vcg::math::ToRad(0.01)) { + /*&& absAngleDifference > vcg::math::ToRad(0.01)*/) { // std::cout << "Found angDiff:" << absAngleDifference << std::endl; // vert[vi].C() = vcg::Color4b::Magenta; // hasAngleSmallerThan = true; @@ -817,7 +817,7 @@ std::shared_ptr PatternGeometry::tilePattern( const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, const std::vector &perSurfaceFacePatternIndices, - std::vector &tileIntoEdgesToTiledVi, + std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexToTiledPatternEdgeIndex) { perPatternIndexToTiledPatternEdgeIndex.resize(patterns.size()); diff --git a/trianglepatterngeometry.hpp b/trianglepatterngeometry.hpp index 4b19f44..5bbc857 100755 --- a/trianglepatterngeometry.hpp +++ b/trianglepatterngeometry.hpp @@ -121,7 +121,7 @@ private: const std::vector &connectToNeighborsVi, const VCGPolyMesh &tileInto, const std::vector &perSurfaceFacePatternIndices, - std::vector &tileIntoEdgesToTiledVi, + std::vector &tileIntoEdgesToTiledVi, std::vector> &perPatternIndexTiledPatternEdgeIndex); }; diff --git a/utilities.hpp b/utilities.hpp index 8a1be83..0289725 100755 --- a/utilities.hpp +++ b/utilities.hpp @@ -130,13 +130,13 @@ inline void parseIntegers(const std::string &str, std::vector &result) { } inline Eigen::MatrixXd toEigenMatrix(const std::vector &v) { - Eigen::MatrixXd m(v.size(), 6); + Eigen::MatrixXd m(v.size(), 6); - for (size_t vi = 0; vi < v.size(); vi++) { - const Vector6d &vec = v[vi]; - for (size_t i = 0; i < 6; i++) { - m(vi, i) = vec[i]; - } + for (size_t vi = 0; vi < v.size(); vi++) { + const Vector6d &vec = v[vi]; + for (size_t i = 0; i < 6; i++) { + m(vi, i) = vec[i]; + } } return m; @@ -220,7 +220,7 @@ inline void deinitPolyscope() return; } - polyscope::shutdown(); + polyscope::render::engine->shutdownImGui(); } inline void init() @@ -302,4 +302,24 @@ std::string to_string_with_precision(const T a_value, const int n = 2) return out.str(); } +template +size_t computeHashUnordered(const std::vector &v) +{ + size_t hash = 0; + for (const auto &el : v) { + hash += std::hash{}(el); + } + return hash; +} + +inline size_t computeHashOrdered(const std::vector &v) +{ + std::string elementsString; + for (const auto &el : v) { + elementsString += std::to_string(el); + } + + return std::hash{}(elementsString); +} + #endif // UTILITIES_H