From dcf2fe5cadade3e9581d68b44e6b9d5db5040247 Mon Sep 17 00:00:00 2001 From: eap Date: Mon, 15 Jun 2015 21:31:23 +0300 Subject: [PATCH] [HYDRO module - Feature #523] river, channel, embankment meshing --- .../SMESH/images/quad_from_ma_medial_axis.png | Bin 0 -> 9580 bytes .../gui/SMESH/images/quad_from_ma_mesh.png | Bin 0 -> 10961 bytes .../gui/SMESH/input/basic_meshing_algos.doc | 3 +- .../gui/SMESH/input/quad_from_ma_algo.doc | 30 + idl/SMESH_BasicHypothesis.idl | 7 + resources/StdMeshers.xml.in | 13 + src/SMESH/SMESH_MesherHelper.cxx | 23 + src/SMESH/SMESH_MesherHelper.hxx | 8 +- src/SMESHUtils/CMakeLists.txt | 2 + src/SMESHUtils/SMESH_MAT2d.cxx | 1571 +++++++++++++++++ src/SMESHUtils/SMESH_MAT2d.hxx | 224 +++ src/SMESH_SWIG/StdMeshersBuilder.py | 48 +- src/StdMeshers/CMakeLists.txt | 2 + src/StdMeshers/StdMeshers_MEFISTO_2D.cxx | 16 +- .../StdMeshers_QuadFromMedialAxis_1D2D.cxx | 1251 +++++++++++++ .../StdMeshers_QuadFromMedialAxis_1D2D.hxx | 65 + src/StdMeshers/StdMeshers_Quadrangle_2D.cxx | 32 +- src/StdMeshers/StdMeshers_Quadrangle_2D.hxx | 2 + src/StdMeshersGUI/StdMeshers_images.ts | 4 + .../StdMeshers_Quadrangle_2D_i.cxx | 38 + .../StdMeshers_Quadrangle_2D_i.hxx | 23 + src/StdMeshers_I/StdMeshers_i.cxx | 2 + 22 files changed, 3335 insertions(+), 29 deletions(-) create mode 100644 doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png create mode 100644 doc/salome/gui/SMESH/images/quad_from_ma_mesh.png create mode 100644 doc/salome/gui/SMESH/input/quad_from_ma_algo.doc create mode 100644 src/SMESHUtils/SMESH_MAT2d.cxx create mode 100644 src/SMESHUtils/SMESH_MAT2d.hxx create mode 100644 src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx create mode 100644 src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx diff --git a/doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png b/doc/salome/gui/SMESH/images/quad_from_ma_medial_axis.png new file mode 100644 index 0000000000000000000000000000000000000000..b02ceabcd1ec527a60281d981bafcb0a4ab79d56 GIT binary patch literal 9580 zcmV-yC6n5TP)Px#24YJ`L;(K){{a7>y{D4^000SaNLh0L01m_e01m_fl`9S#00007bV*G`2j2z{ z5(g>j&X^eh03ZNKL_t(|+U;F!kXzSTel|_0wy8z&OloR;xhRWoMb@3W+q`LZJ1uMCf*RpSYc|` zIE~`$N^yuC+x0r1PHG~4oU5aw`=u-059vtnJN{wsu64Cn&wbDP@x0GF)@U?9My7~e zMz?^B@JdETcqJpe(k&n(yplQX?Afz)69})$*z)r7%*@OiZ@fXbB8rVsJu4mi;0Hgr zX3d(TM~~ig(@k_Y2(JTwqsPVo(A}nZKmjY4%K$9PqFX_DB_q6&#cGvG1wb~NrCUMq zD+RCQU&-RLdc6*SF-Et7@T!1St5pEHuG6g`yebpCl2QCh!7KS!G74S^uN1t}tsuOT zQT$3q!7KS!GBU#Jpc~O#2(M&>SF&V`XfBj*DP#3|J(I~$vPJ$?jmVXQ70s3b8R2zs z!J5zK#Xn1>65SKRt75SM-2}p`8fKa%@5+)9UL_xP8uily@lFiG5P#=PmwY}yJ3C8W zmg)llP_Ce|1alDp*j5zogtl%VLRiLFp-`w+t0J7WEK5c+*~sCwTtR-W?JI)-SvPG! zBw-xHnx+x>MlK%MyLWH?6Jt10#))z(UR*@!aG#wkT%1?KFnD9uG|fyV!+(aXYip=% z9>34Dg=vc&=;yoQxU?)Q;JX1oKo($X#{=aG@(t_lE;3=-Bu6`1@hdlWpn*O1A9kQoJLM5FT9)M=`r}`Z>uuC(H4&RK#=4Dn zDWVOPu?)lS&;+Y>BosJ23DnWCD?bli;Z=q#p>~+_Ek_93nCYePBy&^bn4zoFS%j^^ z0`e`2Xh%05u#R-$T!W8tx8G%qh4j1b|J$3JBu6{s@jxtibz|2Z>~ieZYBkUQ;;D$y zGvFwq9a_AK;57;J!!CEa0=vNzh2&@nuT4CiN*1pniwqruE^@Sl*GBDBtCAI0Fwhlt z{m|*BmPL+sbn$A3t-EX3UE^ifzP%F{aCw~IS2ZI?J6d>s@4feqy>J%0YnZl>W5_W~ zTj{%4zhui%d5saZO?JDi_? zj$zue7Ha*ui6JW4JU@e%{VpeN+z;%Y05tfN&j-}>fb4eYM9 zAlr;${fhQ}b^cg_lWswiw4<$^YFI}a))5ch&k&fB{DJyuln%pOXfm}Yvg{5`%xu-` zbv{4sCI-Yv;;z^0rBX?d6}x2%I5ah)a0XhFVfS0O-r#eW_G|8hFts`6cC4?=t} z)=p$MdJh;1?7E4t5xKAq!1aco*|!u*T6t_!wx&+}nXB3e17 z$BD8S7+zRd*uH%`#jpKqah*ypI1^z*$DMCcykjhmR>btJ;PrlM^XAPDKKLNT1L`GP zafJ=Rl)`?YRt76O&|nAb66~>K$0*+#UA)RY0tbv_c=3!){6Np$3diFS0NJjFE(x#c zAyuB7{2rg~_`SD^CaQc*m`U z?zAboD278^1LKBJ5p7X#zt(>(QOdU(jYi|_*|Ycl`4di>OF!y0;R!B&W!kvHg8^Xx zQ#{~dufF=~z4zYhO(RT<&Jq!?sWlgdvec?JUf+KE?R&oQ{Qw%F+?Pre!M&YJF{0Kk z=Q(agjtyBW%hN7bC2qO$$}4~AC#D{&oN<0g8+7*LHkZZCXc3is#H8`o^GaC5FdB`9 zBd+Dja50L{UH*X~mv>Pq-+J~qKK-baZ%q|2uE`mO=u7z#vK*|!u7JW{Ji1$Bx~4>#c;>K`~i9Gt8cXxDf1MyGGZyiileKKFFmK0FUIKdFGkV zeC9KR*WqCf^oC%^IUd)}3^!`|bvvTAjVb`Jtf+PG!xI<|yb3i!wiU%)@e^DE z5FBGAwIgbSUtB(9YF`5|>^0XRhzF8x6NE%?yr>Ns9*7fW_{LobuZbYZ!75A+wIV#I zpGNg0DoX%bkw6k$$^o~_LDs$T(AKgnj;$C)h8<-u{z%w^gbU4ep5q#SVn@-`Di;f~ zOPYmoKY3Tto^TxE@ix+0CK?1I(|uK*pOu<=rXO^;;*cAo4hn3HsAcz$ZueE;{?;Iv zIPo=2EE8%#+&$gz0=}M*I8EY1cMd5B8 zIaxlR&*$^yayjCv*z86-U~RKU^4b;8QT(dLSV8(nVVHXrKq`nBA+Q7O`b*Q+=FOX* zdFGiBj0e=HXI*#;fBRjSi@-{qJDW6#$vj(Uj1`N;9*Y0bi?wzYZ>3cqSo$q&UO$TN zRgL)d>EkdL>Nr}m|Ghyk9X2|W2a>2^g(*mJPwcMR_FVR|yIJEJ@4J31gPX1!)c{tB zQ(*J&EjTMXyKAWxboe!UBh&_HaC935rC$}Hc%Tocj!8L2%kB`fLE5OEhE@{ zRT8gny@RV+o4^io40jc+yX)Cakja7Y8ZiiiW=|-)Tlrz~t&pfiAo-QE1zM3SIah`7 z3YnhFZHp1iC#=`&KvXFAc6>*4w{3&k&0+{SSHrC-%D-R|#!^XWu?5amF_Fc%X6?jXQ$M=miYo}O z%J5%B0GHjER>WMba^%US03P7-hupz`bp#JlbFy)b-J$Ua!O#H^SScQeoFNs9#Q>5( zZ$W}6+ZA`AxQp5Zgx8QBy;AW^yn=!2$$|}^la-si@DUv+9Gdtw(SGWGVa#n zT9#FL0p@AL7H`|2JJvUospnrsRF}e}p@V}`tyW{ftMIsvikN2u)}T$LKIO2z}vy#W31Q{4#@u-XN&5?sX5 zk$<@EmKR|6Tfg?DZ4(m{gwZw|;w_kn-GS7^8fP_c8(i zz&F0}4I0=|X@2XQ-~8s@y?a+yRsw!T#*<@&Rg_M<5`E`ASV!%N2b~U){9ssYV!Vd>BEvYR@8W%U-DfSbEhUyeb?IwCY*hrV5F2my%#+`Tkz_&$n$(xn%cE|?2%IeP2KahDv-fl{dymtdCLT65dtSIaQ~+cL-B zfB${LrF#B#d3kweW@i8X{lP!x8EE>Q=6=I1+>d&zStY=6fupX>X|G$1gA16$z ziPzU&d+j&>`sv3OsspkYrK=V|NXUxHD;LhD;(*%mt4$_+YPRjrysIeJeLTsjj>eTR z&owXbY1oPt-2wE%DDX{35sdJfdcL(NDQQk*W%5AG*P5M`f$LZrAKa4V$nGj~pf>Cf z;uWzfX?BHIKTrx+t?Rlcd^@^cvAkF`naRH@Ump-==)j!x_Y?zGIaCNw@^Uk_Wca!n zJb+d|O>Kf<@vlg0I94XRxXA;u%5iwXZuc`_{yaZzKKaQ{_VEoQ%FhxRB+bAKNj$K1>sGeEW`3{QHqXtM`+O1EO`T#1pO$QR|H`wa z?i<_yrMR-jj9!St~?1%sK>mT~ihX{@x zjrMH}fX$n=T8(Z4;k65MjLO8Qm?=wVvIC8_-Gubp$pW-$Ybfyg;)^eOD-X<@ z0L*zMZ+Gc1Y)!)EFpP-;K>Y&s^R2hudh2ig>0Dbk2x0jyJ#_Z{`mynv`6N6a!UrCB z;HzK#D#1}<|9aI`R~0^cmHCp^s{X*%do#Aje>u05M(kS;0_^$cpC>G;jaLBLI#@eV zSu)H;11oimYn<UO4FKCpU>%nFK-;^VvHy~sUPhQrMr!>=Lutg%g7Bb) zLV=(7NrP4HX_YI~;uiH+;mv0Gdx$~aNZKU^;z#%t#ZibD~< z+H9t#{%fN8Ebq>vt>4NDzDtZpwj<3dl}h68Q4lFc zy7>c8J)_wXIkmWk39qYGt@`vQZe>5b)qLF6cmr&{OZLNQ5Jo-sY6Sm)8=&sTNz-5f z`8fky?0k;8218OGn79V|B-j_6e3x|C&3r$M+vED-IQ(`lLq~06Yc4Gz1`px24`%4d zZiKlAjB7x$dsN&KoWTI)Xznc_ksrD0Nt!Ph7QCW=3Vv?6!YBZBTj5UIzgCtg9vC9L z@(7QAv?=fX;j$|UNorP4HXn>;R6t&~qp8u475W3E!wGiV9R)1Ln5Jn7fwg`bm8Di! z)?BD(_Tj=?1jr8dfdTKSwED4&2XX15b@uWbv3Ao0)?OX|wbrg(yL$C%zMT-j-FM%; zY15{)Yu6_6d;RT`%~AD#{tAvhgw;2q_Vj-ofBC-8ZY59#*g*GVFTRH>|Fpdy%jXTh ztnq-)uh;4+uHDtLV-XJ|;a^)i%&x|)67){OtIy|JuV?Nn^>Lqk>=Xh{RPcRBrmYM6 z?!4yeA0u2UnQy)H(o4Z$_3O2IE{AdVSMBQ5d06RQYdz9 zhFKWP7j?@Gu$EBaN}YE@Xf{T9yuu|wnL?eC&YLw z`j$DCUny7@k=^7xzCAmx5r?Z1a;1&~1h8k;%@e>x7SMj_C*R0ySVI6+$-gcyFVD=( z?BBmX-k0_4inUf|wwD4@(0Xxl333d;_<#M``SYZrp?JO}K2bi;?e4p}o37ouM*a#M zvW&4r!3x2GeMLGM=1l-Az?m~=uD||z!f2Z{2HT!^wt?NX*1RErX^S0b)J`RPUG9Tx zwOVqQsMTtSE}4KlY4x_|xjN5JmTOcb!8=}Na4`efqi2>0CGN@j%2_WZ`e zS~ig-PpV`*FmVkZAJxo78#n0){PL(EuB}GAA|dN_wguDc3+8jcxoWEhGOp3ktg7*= zwhqNzt&tdo7)*k$_~RDOew9d!sYWfO{sbEB3le)**cyvOp$t zUqDvZALjJ3AYu{u+w-LQiVgGvj?YAF2dbEfC z=&r*%nqk46F!wrH`ISBM{kUib!o||WfqL=lJ-`1qM<4ok=ag)YvHJ>z9lmIbZtiG- zAMq&=3x!YR@efQlr;&@eKzhU|T93PC4-Kw8FOd@v9 zrT>*bcjFXT=-R26ZEdH<|a2QJ5OSMj>*Zy>zVes2Fb#4vaO zV47yq&b33IY=h27W|Q#Rck_g)lVR9gctlk#CGpx9wf{^e6Jx&R+0m6N5B1Z)YM6^H zhGBLiVHLo@&%@ZGBc~edKm&X1Kh}|yt}C>*MTsE#42-++#u)P)ypexv8jsgLG zym0k&M3tTl=D|4HTJDp&SS-dEo_4I}%Jdha6Q~gnR8Ka+DyhKcAGRe3(201JC=Vtf zn97Abu^WdNh8X<*)~QpcZv2&pVw8jr=D|4X01mW;3}J((wYqJO5yMa=ydFAqD3{B5 ze4{BrfKJ4#h;>t$s1-)R*bk<97b|$B@_TVQ+JVJ-&e%%DPmD2ihRp9O+IOtb(g0Z|wF%YC;4q z5fY3ziCdX@ZuDcM+s2Bd5Feb$B}X-+!WOs7K-S5*hMC!_)oR6JF*!~QFKmed-^iLu zcoiR5ax0I!EsgP9Dp8~g6ADH7Ig~3+)mXwSnWG%$*FjAOrXIFZ&p%wx&tdi;F1|Nf z9|!TXdU310b!8)3^(3$&_4#tldE>FDKGAL)bSNJ zA;HldNG%?4w$OM}y7|y1x~iP~YcF=}*s=Zh@5P}*5dr5&`MZ4^g+f7yyD-L*2CKW% z%NWzQ7b>r6<{~gv0KivzkDLtmfb$z(fBp6EfB*Yze>>j={n<#-k@ynJB^-0)j-dwY zNZY7IYNQ|9w-L(LtXVTPHFf4kz+FXfRU$x);0=>Ge5zk}*U(|;DBNjROGb^vtxP>j ztmx*JGUcUpzj_sk@9O9e!$TNN#+EYuEklR3v-yc`?VJ7?Kde1B(Ds1T+1c4lCL{aQ z@NF~Eoh}`wd}|t!y z@En-2EUSk-9G+$?Z*m4qj$woT8-ID`f#3ebhd=z`p}p@4)!wWleUNwwR($cBC3H_9 z#xi*5F@%_N_Dc5b*+cc$fvIPmDD!{j=M1(KTT~E!`H{m}x)ZO2j((4v?__HyS<7Pd3Gf@=6{@v27VrRYhDGeqx(N(Y4u~NmiTpMgQ=&xUzAiVa| z&YwSj=zITyr|NB#mw6NOW`1uXR}3H|I79_2m5Q4U6Hy(b>jxvpS`3(p?DbbuzSUQ& z)oPDE`Y7ff>lS1Pv2zR21Yk_lOyTJFM5-+l`;vy0IyUG4_%f0qSKb`_5c{2X-g#o- ztW)A8$FRGWe%8e)FP=|Jyy;Z1tf#mZ*W8pDFvNH1v#naS>hs%Iu~#Qd6k1=MDzG~S z!j~i(q;_P$XzN-@=SVE}R^`ERNEWjLeWTJj0*~k-#$i~69i8`xXly~Pw<{pqN8Tr6 zOn@nNF05r)-6FO1S&jfHmik*=Td`Q=s#* zR4PSBs>~{k!hjLcytz2E^cTWwB#Y8_h!LqG_Y#(WxSrXE(!WvtwXdd60)rvuFHjYR z_!9fk@BPZ$p^8Pxz%Q>9ez5%Vo_n`beV{*lm%eN@ib8yenDuwf2ije8GhpWD=PBPx z+nzn%q!2Vm4)b${g%}27JcnC?cg`RL@g*If+V&=wfOA!+P(zMgbKR{J52S4`z2;Jn z%CXAL1vIxB!D|V83QU*29Q~fxGVd?B!HTKE`FDqBQL-X>@b&uc-Mb%q?6IUT6o!N# zo`e9!m(-BjVNvq5g)PP~4bZCLUstbQy>{(dk8em@jRLS0JIzJBrT~g~wSFO$U#Rq~ zW-lh;w#Be^4p)0tl|cK6XB+OVM*ZzQ%G?*Z=>;!cmKZM@+yfhK0T^ZQdj9nW_B61& zhG`2q##e9+REZ#=tiQMN!m!c_l*zZo*P!+_ln!fR%X+^gUILQ2ByzeU|V~ zG`_3u$}TnV+GsTX`}tRI|D8uX_#A-Tl1+im;@}uKuIN{SP1T|DD{}S1r)J%&b$;-J zAKY@wEffzVWtT2p`p)AIVgBpRZ_Jwj%=s?h#m!AK@T*33gp{Zcz^er7=}_a+C+rYd zzgWjp-i)oAHfd1Ql<2};x^(G3o?Es%vGT2-wycGkZ2ck<@TWrN3?$9Y4tO(eb;OqS zn{*M;l!;%@zljH)g!y=L%VCc3Eir>YmXjCEa-gK|PPUC2?|tl%$hUZX zpt1yWp_xC8YXfrYXJ=>KadQe0svk_nrR%3rIt+8MCBWTZY4XvkXBj#lS=M+pabW#K zj0)P}D`K8c2v;l8Cjn^kHYmbtN+uuugDk|C99VnMLpJo57J<1){&nQdL3~M~o*GfM zTc&F%EeJ<=jmMi7i$ww1o{0pp!?dmLHZal+BZ}97qIHTc8pl?|Oe={c2(ns{x_DKo zJ^-0Q(7*bCXM@)bFYa+;foWU;@$A7C+dV&pk1N>WQjqTwO)?|=Di~pq#(_`(00iJk zL_t(Z_jr?rseP5xt-_YM04!&p+kp$`o2F=39V2{pYYq?X_w5%G3>4l6njG`Xb9qqn|;ES zHInAOi0k?iUjR&7|8R)(2NJVg+i~G7%|t=0n__opg*ykS#MVsR^7*`j_=0peuP!g% z^5yFb=J%>L_k}<0TZ~QDQavltqwO@RCm~wl1Epjpv$PyoNj5h1$EcjW4F~r{Vg3Ny zI_Q1H;Mtg{A!++at-z=seNk{_MNAp&|(a}f?XrYDYYK2UjVtL@zB1afQcMa}0p>hb{jg<1aNFct#zZOCT4y`hwB6DClGptSim1 ziffJlpVH=c?ASs6H8wl*!*_1@t1H_iv~mo06$ffiwp=a?xtHWp`<@@9QmJG<2~Sz4 zcfRwxt5&Td{~DWJf9;jlPUPndOg9%@Ok@X87~Y*Z^#!9Gpq}}hjnNG{CbAgU0I;nX zc&sW|A24)SKz`0{aw21C6;2eMloy~-D3lJte98Hi^hwIMlFV>9 zGkBGXM+HgW-T;95)g0U@mqH3f#2<*?G)?2hUB+09v)zv8ftT(4oM=Fx%-8G) zhnn4}Qp-R||2p6ua7Pt3I!mhuCpxS~{xzNxnHzVjj{K%lsq}zXzt-1)R2vAdT{u@k zAut!APsV?ud3J7NY@3l29U;6ra)ta{*r>qm2xCE^+&1GqGlx?QnC&JfI^=jQS3G)& z<^u3G3OkZTVfeP0VHmPVd`Krs039~G9zA+A{|8Un|7h>fn1I`nkWnoJp3Fy0d!RG+GsSs{`IflF!j6Ua{$1*h%J98_E`EV^0EEW+V#> zH!F)Z924qhNgk#+r@qC@SIV6EB)CcgmIX^oKjf*6G>-s-jV8NGASR(sk`%7%#c=uhG)N{Vff~I_Ar%AD1Vz6A zOEx5#k$;eci_N*f25XkR8Mf68gav_@UBO<;*^_(xB&8mSXd(^p0>k+MJJJt9^}q#@ z)=uh)&{&L+VE9ohf9O20`AI7LTBn0im%k4}GY}97AtWT{>0Iv{rGP5?E{kevI{W?2( zQ9&R-NH>6{&c3GhZ=kPi6k`{(6#vHBnSDezf(P z#y_jOL?LUoIZPIQK6hP*?^(DJ#=w1T80dty$gG4#H#ED`AY!9-fx$HsUuajfFfo=? zB+1OIEXQo|l#~kgBH1-AJwc_cf8JV`jnlUKINHJla0eg(ddM|S&Oa1|yf;$mA`&{X zVLu1#io6*Huq>X#im)0t&Vm<5mq!<*vUg4Sw`N>P|CG!9kg~c$jdrpx8jM5>OLE|^ zAlc-ZwANx`pL%T8OFR@S0fMiZhHQ&1Jni)t^}fVJ7kDxw_lGDaG#IuP)>UD;YaEQa zz)r1nUipK5&Vs5ab1s%3l~lSX#M+*zu#@3jzwBEd$)y};0ah-GA)3q7(cEz=-oa4j zTdheq;R1R!lM>bhYeMo&vrIgWUx`+pk>7Fsuy#?&lQeID9@SQm?zHqyr$E=3;XM~E z4M3!X)zamM4gX#+bc=5R+SM}R98w|Z!puWwBvWm=#f1 z4Z9^A3FGB;IbAF6vy0psujU6{sEW8ooBPSMFP1h+t#&j&=^sS6Ao!*v^c6HDIOkO3 zd_OoLJhxs$fJ&?r{5*6FVf9EM`ZUM46#~Q|6LdDiECL95QYzYz!L7Ind)1$Ibz~<| z4;FPQkS>0Kef$+1C9GF3#y9-L7adkNE$R*j@shDlFxwek^$%~1&WJzW+RO!t7hqKq zX<@zUPn#ugA}eGJr_HsP(vZiU7H$@Ib!bLEfO;ERvzOb9Mn&qXT=^2b3jhV`27qFd zLj@+<%f>Ha%HR$&5`8*y{FKR9mAHB63#hH$el@DVI!A4^rr>UA5T;DdcbiYll zS8hojudjO58%y0xC`HxGveRpnD7mrH+FZ#g3LcQhMM=PXpcT~GHL~7Ywo*6B?oWuX z^Fhu3uOI{(f)C*|!ribx-a~QTfBE!LsYwvP)o5+enm6*SdUJ*cD5# zi^Mhu=yLks`tUp)Z|=j;0)CQ~*v$p!$l$|8T|xz!dh)StBH9Hi;2f0a8jwea8Okmx z<0OST5ewH;v^1IcD~U@D&~y|rpYUN|pHQYmbD>m^szE zc@q_}4U{5iicnz9g`{j$*NFFa#dZ;UXhy&ls}gNP$p5Jwm${>%`n0yteeA-@t1mI_ z=OTUG7Cg;(=u;(faOcPGj&H!d_a%vAk)IX67NbL9)$uUzbf0glP+nK8SN@s!3Rr%e zA`Op=>C-wAoV(9r%>xu!iqT1})S6U8c+d044IXL|)TRl-I+u-ak``LG4RSvs#!JTV;U$%cLbTzp(OI;6D_HO3lP12<874G=PO$?o6Oi-u3#XnjkS);6cZ|) zRankFLG9I@-vvX&U?3Rk%*p=NTy1h0B(huR=h7qifI*QPDtkJ$sVwVpHc$t^fW1n# z{@TV=^=Z)w2SvDHz|60{zFcF(IKUmd`>R;sl9uMcixv%K1D>_#CpVpkvcmRR5rt61 zM+^k9A?a9*)3XCF)MBZ<4h1HPO)8ru9wxl5r|`n?yPl1a3BQYS8ppqoORUkS?K(r) zV<714F!}FvO_OvSqjUp^*H5r>3Do|v4t!@t^TLLpNl_G!l&OL}#)Mtb^Ux7osU*#` z{5~WAkP%>rE%jxb8BgPyKr8e&O;=j|IHFumF1{W=lL!?-E<7non9TsH*LO0Ves5;j zObDgLyWag)ploKku<@r*d1HP+=!uE_DwK(EiptikmV5V9;!d-nLgSs2yz( zb!>fkiHdD4taVIkTisHUFia^LX|}GfiA$8KY^G=X`Zn;Po0r^ESZL7$X>&-q+Kv^h zT;tSYbOb|!=9w;g-Cmb_R~QQ)fD?z+F`|rHJhmKE-vxWkgO0#F%;UFX0$1Fno6V-6 z=U)&^gk-Sz2ffx*y?QS=*sOA9a=SV~Y4PlHC*ZEESvYn}mZKui__=`Kf-sBR!h8g( zcC)>1XfR7^*w<=dl1mZRg?>CF?kNXjHaUz~o>}IRDIk#cNBPA}3Q{HNEb+ly_bfP1 zfn}8#g~mO z7Pm?*ee!qMIKj4{=L*uPJmMT6|NF;y&h2zAuqNhX4et~YO3V3kZ@D0f-MZKmY3)Kv zl5AW&jBX<0i=lhkvUuN^Kg{9R!9filu16_S8 zKU?&sl8+;8lJq&jCOxryr;eloB=;pLl~kYx-eyuQ=fsMFUg##9ILTDD_RMDvBJHGn zzne<8-O}yEP&P*r*&t+BbRx+czYp>V>!|H6^D6t8=q7BR%b;59xQ}x(o$2PD50RPy zj|xzWPcze=poi2yKMvhhnG|ZTw>DdgxpgaSfS;W3Btr%ly%Ny^17xTX58}@Db)~%6;Zr z9)86vg-5R=u?7|gOpTmyax>mCI+41^_nNN8@vmil8`_uo?PC-J0=h~R=<5`n*dP@S zER35Kh^fZwKj3otGJTDRQwPRw$7IGQ;D7)HXfE1dZW|q~{TC$H%bu6NL(DB@yS;>M ze_sD1aYpY|ObI5P8{E_%oae8n;+@x#-*ipibkw1IRrX=#%MczR5hKmv^^?Sd{@U#T zF`g9iT1zc>@s-Zo;@gTjMyzK6>Qns{%okC&sj&;jXSIdN>03qqGY7~W+q0H!66T=H zuIs7rZVm#^XGvOTqMLTEZv))kK3e^-oZt#ip=0Mj2q5@!bwk?9M>s*l{FOoBPwTGK zwE2BL#vqcvTP(fH=S1j{QoXK@=h0^CpNKiSf7kEq}s361EI&CP$Q zB@tdLwK$hPmvZF~6QZ(G!CgrbtW?+!UCNhburZ!;&iA!_Oy8VV&NQB~9~%L8^SY%xM?=4gO)0O;z( z8X{z`_QtK1wx+y(=b|S~^Ty6U79b=rIIV#E;nBAkSz_Y77x+6~_wW|4>pE8Vn~*+J zA=O+utiR982J0DGGk(~JW}2*;eQKFoL!a8=e<*5`nUB_GGHS`M!?m*qa{&{o%EKQu zdQ62YOf>Y{nCqpak1aHFduV>*lGNsDFk$xd;jMYJN@9cSFR2C*lyAOyj9OwIh4>h= z(^C4=cY%G`buu>h)H6BRvsC%koD+udx@6aNItdP~;>dlyTN-3JBzf?=K3JqB8oq| zFj$O}Gu|BK-^uOpi_dL9SE%BKG;(q^f`77sMUD%FQ3!I{U^|XeKBLR$0dl+lKr4qg zWg_hC&wugY1GK`@CTMQEZFTV@rJTmEDS>KBf7Zw za1_l;oe3BI{x@DVHfsMFQbGDi7bBjuiP>=Blac&h5d!3rz4~%iZBvLersiJ~ZC8K* z7HEVFK^tQ(lpo~UVwKVSF)WX1!=1C=eD}V-e>-3P38h$a&u4HR)&X<`y`t(bf55J~aDuX{jao1G7r^L|S4R2l$A6`(BxBVqiw!SGL(x_|`|1(pBg^v;9qlI{;Etjf>21sEm z71ESL)s9-sbT}-kc>x6N%A|Wli$_YKv>kd#aRJbEwhnEv7@555H@y+eDSe|yg^CXZ zN0^!r=AS=y7`kq#UZkj%cD&Qi|BTl11HfQc*72~EI@z7NZ@5Ux&V!&N7C;zCC?%lY$@4a=~2v87(B2oo#l^Y(0mpw7oSrlr+@|I;^6wZ1>DAZiX_)g(p#X& zQAfv6QTf2fOIyq5ItzL6Kj={2_Chaq7C*JMeF7<;Z9v|5j?M!>Cb8Cz)b9I;mcaHM zqf87d6z|jKAF>Hy_LimWsDyyl4of-jh2BCBR2iVf9VjvTZK^IuL!N;SnD#8XVkT&& zN#vUxj2d36zl%amiK=yIG`m@Y-gpB5T;Sz7qmwP<@sXmuCyfJH%5pZP?rs3Oz=}}f z{gyFM&dmCBFRHHIRYx9L z0<;t;3DBIVyhe@gBA^Sl4t4*kF|B=bbBGbN_eD927wAeB_Fz;p5th6Ek3R$*vcz;W z%`3K+-$cF;XpZK;!E)a~FmIiBykwG;j8dAG=p^`Kx7zfb02=TLr?^(cUgsk2<_Z90 zd#@XwYf6~%s}KS&gOKlsxp=(*e<#-a29adO8l4Lb4@`ehfbJ6QXe3ggn1%xatcwa~ zei&f5BreEPGA1_==kMyDw?Iz4_bA;azgwW^oQ8!gS0)lUCgzCc7z!27>{tWfUQ%6 zoIG)jkW`G0O8F%V+v9WJ=TsZ)6yaTj=y|n5a9=rJAN1o;xLUzHSa*qsQ)xWcdOZ0i zJBGfwDnZ#)8WUm&Mlf#1AVuo-H)Pg&ojlN5z9rb+ypX`FV(R`^j&&J(ssbF&Fo%u2 z7+Fl|E69tAa+2*@lj=}}KI3 zz$@A>Y<_ztSWQtXF&y$(+ zAQ$0lbfKqW&|qY*$|=CvJ!nzJK@6sZ-3912V@$cm8epE2MykIj6| z^eq`Bnl80>{!g;qMw?eKL5&bh(I(~3NBh5z(LTH+Y{p2H&l+GhW^1@^ZHdLiG(zB$ zKTU&P~tbrYS?wC|?)8=(jaw)N@Q)r3= z&m_fYEh`}cQcl^voRmfRQ1W;vB~r>c!v3$y+eN2Cd)IXUaf4P09Mo|DtJi7@QstM z5&Kmasc*zuI97sUK>syWubxP?U6YX*XmEx5hVD)OPT)ZQXwlmqRq^*kVQLuq5pplT z!2R0PwS(z)e8vV_Ptf9I`MFC{>54~_ucuXO0(1sUnIAra&NPdu_E*k zwv-q~lP07RMbMYlnwfKZrM!4E=`8$((x;@m{%>np=zddR<#xsto6ezzPTSt`Hz17 zUuxKS=?+-8@IvKHpvO+0>UlTa;Gx3k3)6H+uRTQWs5I_<@g)>n<1d7Ez>fGvOsr)TX zo+%*@hVv>yU1!QmX4PI7G*fUGkm!{U^BBV9Yui_@A1S#xUBaiQo(kGYf**dWf2?Ud z8u?P`ecrpoJ@q^$$5-(5W9HOnLylKXL6~l$5wv($hgbYf_qVuBDP~nK<(OhpK2B$5 zx_M4Gmf)q<))a#1y#vbpbs>nJl}+j#SemSVB9rq~+N zEYzP}Qv2=@gmg+1h_koa7BqE7xfzt9u0K=y;Z+LG8I`GO>JjNQp1+B-dT9z-j`AmtFlrEk`YYj9T(ObD`NpD15il=|v$yK(TfkQ{*}OS((iBGI;gZ0}q8%(oUeY7V6Zc zM&H__wl&z7WpN~O4m$gVY{TB)OYNE={A2Cz7@8${5J!z z54KM(MZZjB-Q(w3>1f7!E42(wzU=Y!O!P$RG$tq+_GQA;nqVT8DzlOI_@=M>1aH&- zCH^*1QG`0gEN2;&edEv%=b%meRsF(pX~1GuJ;r-&-8N5)SBZ<>4rAjyB>9GC{_1A+ zI&@}7a^uV+(KqitZDg9|!yDdlAC~*e4p|83!wjjC#Depr9Hy977<;ymlUI5NBY24z zMA<0Xc`i3W6lf7P64@xKnB0sLl*=S_6MrtX=6?Hwc&uf!xM5T6k5gvZ5x(p>yICBk z3jc!WKVHMiV9@=~0q`^WOoYq$4wKeg=NCV=Bhs!1!(Ok)O7z_(E903rMV(i-M{6;xhm(&VC$G?bV~^oPf_uI}|8(BrDYHcjWY zquOQ4(8A<#q8nU*6Oy9tP7p?#8DFc)6zt$GI48Yjq3Yzj2%-*JeOgSeJzdtTv_$9- z+*xlh)e)9znIKain3p-8X%5*@WF*ottvTVVOsO9%*)PD_J-9m7Qo=?2jK7{SdmM1R zletq|(qcXeI#%T0aZGUw{~I0R69-z$VB=$J6cY$C1>%$o+KlYLif7l0_Rk|`p6VZ@ z+!~h=2~uOW9-yYb^v;&R(?z!%4*F}=d=n8$sSI z=TTFG@boHgo_`hGtK(;Nlz3d!Iz7y=lkmKW)D|koc20U8o)OXs-IXI7&NW$Y3YD6B zhrc<{&CzB?5FsiU7dEf=(v0^bp&xkukLoGc(9F*?)Mgc>F2!2q@rDvkIgZM$+A-@Z{A-)%j^*@mM2_&f-QI14exJj@h7{3} z2%D~lz-Et}6Q43~wj_+&9y_M!*67=+j!J{T1HF%nq?T%V;(B<4PY-__Kn%X7`!%c# zCpDWo>PLo;%FE*Pk|IpgUzp{szSmCx-MF`WzrH4!j(wq?nF?F>o$3(TR;}tLPs;5g zycK)-`I1$r=A2XiP0M%A-=K+#Xv@`-q4lJlCIUb3>;ByycdFUgl*Iu49bADj2aT=% z+?Uvi)&yIVGT-Og(c@|c7fKWFGIteMoUmU}gXGUNa&bnm*hcr~Tb(LQq(1QKui`j6 zW8B|GiMC{a0moRT^f8}C8Ek&;%xn+yH;+FTsj2$akm|9_Z)Ppys)12zJvu8*=7y(j z`J|=^k@Gk(++gvxGE}+-vhpJ`>>e3$?jW;Hi{>c65HrRn!KtqD|wyZjRokBUdI*Un7FG>v2^eNqeXdYK+g zci<6_w~}Tj;xr-=$kyu@&@biWvL_KW8t9r%x6IcR$%sR77;>s+o&CoB&(}s~L3H=l zjTuPf;XTzX{B`E%wzr+FIrc>E7?e991!UD8V=|GgC!T@I!uK?9`_q1T%slCi9&^^w zjOh3%;#_eMRjA`z&_<)Eb)wdnZmoJ(+uc#z%>9>9HAvo!+L(DFVvQO^FVsjzCLd7< zH<{f4z~96UULF_j+f3eAZf{%}QHuYrmpqgw#1#Pm|LuRZ0B)_#b0+hG zvqUdrW9W4&IJVW;3G>beRBH<82V3l$sU|2AX8GSX^r=TpCBJ=x=8XPCtBRm>D9rNa zY~s3Kh)v$EioS*|V@JGDd^4cZ6)gQMx0Gndc!qAdHX{6h{CBV(2Kh=2ztY+c|6j`j zmxuL1C3Pc8PAlbdM)_A1u7wzC_VG`vUlPH2i`__%-kxI2P_VG7mtDj7L>&~~!ysgF zo#Z!QrDdf;JLMVqx{u!(r;}VD6xI^tR7OYvL;$repSIC}+^ENOlHJIX z88Tg;hhr)X+!j_TT9fHfwuS zGHPUZ(9D$7*Q{=)si+da??3>Ceed~SJwH%|C~Y?xGmN$OCM7S;p2iWsGgM6Tje#Yo z} z9a-%bf5y6=J6!rT+I&_S<)j2>2gm(Bo|id+w7G|{78FLvzY%d{sNi63HrL`?++P?vj0b;7~_Pup}M)GadTI@W6>d!r4yt^oqWFjnH&zPH8gF z`E%Yj^E(o8c&eOz(VD2O#K`yL(K7X$f|H?!#GRH*(^CDJCZP%5IK_Ypn4zAG;HsI7 zxm9fp(puhKkECxZV@ROUZ_M@Q*FjNk3W(J<*l z9jN$Rne89gn7SnoW$c=8CJn0Q{oH=6&&jgI!ea4I^7*4bs^b||yXVW4a^WazzUZQeI zVDs$H`;WFX|ANw&NsrN`LVT5WNWi`V=bdc-c+V-*sNg9CzW5*2n%iLVDsdYI=j}O7 z8<8Y}jY(_X8shRpjLl_phNykuXaB!bqZ~s;EkLpu)p^p;(>P`&U(U%v5lEFv*Z_X2 zhTrYeo+TXFPQS>#yfW_0@sK4_@Z_>~LzDv+sH5Pu)!q<1Bg_YtI=ae73E+jE>oAlB?b zmwO|Mxu9DzDzclqUYQXkJ%ia_(JfM@eNTwK`!;q_WP;6FsMX5~?gS~Aicy|JvDILn zZqKxi)pxV*HH*k*UDQNut{}T zEFimNAQOi?HCyXqKP{Rq#1Kz99_vzlNg4J$!kvaa)Rlj)F4^tOqu*fOSJQ^`Lr^Bb zpsJn6b}r>*7oqYMWJta?GO5nO)S`#o&qs|+-?5N+`tV)bL3zlUH}%^)d_iSV(4^Dh z?FMnVknfY=4kd7hD>ynS+0fvAOgYUbuBy9z*%bPN>{RbjWDfg|^RS<3us@meh-K+x zjW2pRUR1wKL0eO^n70u7v(@kj@8RLIrmYN*i~(Pp^R}(3o@7ufun%v#a_om~&JR1K z^4>=5zEaVsykx{gJ^|_X5zS>3nvEC0%KF$2_(sYz}A6k&$M#>5zq6 z7}rke=$uZ>AE}~Y;7fMn)T@1Oyz0aB_OH~38a2dbURaC2VjDovJ+<73X?0hN^1vNl3$YeQp#E(aC2LE2k#K*&DXZLhVknWF|r-u#0jV}Tbt`|?z|FEU|Ks0sQwJ0CH$?58rFcE_7dqn$mGwqs zkiWUjvz1HHQsn)s?Md(>WIIM%_RDOu_^MN4kTbVvopPSpK$|2$P4VrsZn3DzX5R+> z!p*6Y@4sW0RFz(?l zC*v!z+&;#V;n_sXx@cZSV?~MGmFw};`IGx}=LKj$3tI-xY71AknY~yMlJpTBfc3xE zEw@+*ER4d|qeK<-0F8O?md_i%&G;&yfaLrM&G5GadZUp-mv8`x<{ll^uF4Q+uM%LA zyy6VQ2fD0g-ba`}l_@qZjTxq0ZI1FmrMPY5usd{?uHXK7XL2#zzc5nrU;QEgL@B4K zrNMif=o0Q>A3`ec(ouyA06%MOeG0|F;lBg#QCeeiG$6V-Gu28~K-I`riP%9Gvltzi zQc6}MV{gb7Xsz{qwapW&iw1;Q$%=p@3XOC(!exq&M4hT5-?9T;YZ?(cnuStzj9J(V zVLo_yDAv#UL~|`0_Qm_$l`K)HRC&IjgOGczQkG55FHrcHNpBIVh6{KN)+|K2^~UvB zwG*L01kP^gmYyyz3KoePg)0ST;sgGwSKBv-)&t&u^Qy0xBiK-;`-7=3BlzETDtHKD z|3UeRvK*&y|4XM|#%_Lr`lk}~gmJ@fshVUcW(Ck=RKlUv_2Wst3vVzcu=v?CJHPqa zRjrqlkG%!19LJ8=jpH#j<}$x)#4XI{T;6p;D_I_?H<~lnSVVwUjDf HSOoqb3W$dn literal 0 HcmV?d00001 diff --git a/doc/salome/gui/SMESH/input/basic_meshing_algos.doc b/doc/salome/gui/SMESH/input/basic_meshing_algos.doc index f494926d1..b82b67f89 100644 --- a/doc/salome/gui/SMESH/input/basic_meshing_algos.doc +++ b/doc/salome/gui/SMESH/input/basic_meshing_algos.doc @@ -58,7 +58,8 @@ without geometrical objects. There is also a number of more specific algorithms:
    -
  • \subpage prism_3d_algo_page "for meshing prismatic shapes"
    • +
  • \subpage prism_3d_algo_page "for meshing prismatic 3D shapes"
    • +
  • \subpage quad_from_ma_algo_page "for meshing faces with sinuous borders"
  • \subpage projection_algos_page "for meshing by projection of another mesh"
  • \subpage import_algos_page "for meshing by importing elements from another mesh"
  • \subpage radial_prism_algo_page "for meshing geometrical objects with cavities"
    • diff --git a/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc b/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc new file mode 100644 index 000000000..976783d40 --- /dev/null +++ b/doc/salome/gui/SMESH/input/quad_from_ma_algo.doc @@ -0,0 +1,30 @@ +/*! + +\page quad_from_ma_algo_page Medial Axis Projection Quadrangle meshing algorithm + +Medial Axis Projection algorithm can be used for meshing faces with +sinuous borders and having channel-like shape, for which is it +difficult to define 1D hypotheses so that generated quadrangles to be +of good shape. + +\image html quad_from_ma_mesh.png "A mesh of a river model" + +The algorithm assures good shape of quadrangles by constructing Medial +Axis between sinuous borders of the face and using it to +discretize the borders. + +\image html quad_from_ma_medial_axis.png "Media Axis between two blue sinuous borders" + +The Medial Axis is used in two ways: +
      +
    1. If there is a sub-mesh on either sinuous border, then the nodes of + this border are mapped to the opposite border via the Medial + Axis.
      2. +
    2. If there is no sub-meshes on the sinuous borders, then a part of + the Medial Axis that can be mapped to both borders is discretized + using a hypothesis assigned to the face or its ancestor shapes, + and the division points are mapped from the Medial Axis to the both + borders.
      3. +
    + +*/ diff --git a/idl/SMESH_BasicHypothesis.idl b/idl/SMESH_BasicHypothesis.idl index 274365bf2..60f603caa 100644 --- a/idl/SMESH_BasicHypothesis.idl +++ b/idl/SMESH_BasicHypothesis.idl @@ -1090,6 +1090,13 @@ module StdMeshers { }; + /*! + * StdMeshers_QuadFromMedialAxis_1D2D: interface of "Quadrangle (Medial Axis Projection)" algorithm + */ + interface StdMeshers_QuadFromMedialAxis_1D2D : SMESH::SMESH_2D_Algo + { + }; + /*! * StdMeshers_Hexa_3D: interface of "Hexahedron (i,j,k)" algorithm */ diff --git a/resources/StdMeshers.xml.in b/resources/StdMeshers.xml.in index 350b9ba8c..4150b42dc 100644 --- a/resources/StdMeshers.xml.in +++ b/resources/StdMeshers.xml.in @@ -311,6 +311,19 @@ + + + QuadFromMedialAxis_1D2D=Quadrangle(algo=smeshBuilder.QUAD_MA_PROJ) + ViscousLayers2D=ViscousLayers2D(SetTotalThickness(),SetNumberLayers(),SetStretchFactor(),SetIgnoreEdges()) + + + ( hyp ); + SMESH_HypoFilter hypFilter( SMESH_HypoFilter::Is( h )); + + TopoDS_Shape shapeOfHyp; + mesh->GetHypothesis( shape, hypFilter, /*checkAncestors=*/true, &shapeOfHyp ); + return shapeOfHyp; +} + //======================================================================= //function : IsQuadraticMesh //purpose : Check mesh without geometry for: if all elements on this shape are quadratic, diff --git a/src/SMESH/SMESH_MesherHelper.hxx b/src/SMESH/SMESH_MesherHelper.hxx index 9b17d6bb9..05173ae06 100644 --- a/src/SMESH/SMESH_MesherHelper.hxx +++ b/src/SMESH/SMESH_MesherHelper.hxx @@ -236,6 +236,10 @@ class SMESH_EXPORT SMESH_MesherHelper static TopAbs_ShapeEnum GetGroupType(const TopoDS_Shape& group, const bool avoidCompound=false); + static TopoDS_Shape GetShapeOfHypothesis( const SMESHDS_Hypothesis * hyp, + const TopoDS_Shape& shape, + SMESH_Mesh* mesh); + public: // ---------- PUBLIC INSTANCE METHODS ---------- @@ -243,7 +247,9 @@ public: // constructor SMESH_MesherHelper(SMESH_Mesh& theMesh); - SMESH_Mesh* GetMesh() const { return myMesh; } + SMESH_Gen* GetGen() const { return GetMesh()->GetGen(); } + + SMESH_Mesh* GetMesh() const { return myMesh; } SMESHDS_Mesh* GetMeshDS() const { return GetMesh()->GetMeshDS(); } diff --git a/src/SMESHUtils/CMakeLists.txt b/src/SMESHUtils/CMakeLists.txt index 486454e74..8c2f3217d 100644 --- a/src/SMESHUtils/CMakeLists.txt +++ b/src/SMESHUtils/CMakeLists.txt @@ -63,6 +63,7 @@ SET(SMESHUtils_HEADERS SMESH_Utils.hxx SMESH_TryCatch.hxx SMESH_MeshAlgos.hxx + SMESH_MAT2d.hxx ) # --- sources --- @@ -76,6 +77,7 @@ SET(SMESHUtils_SOURCES SMESH_TryCatch.cxx SMESH_File.cxx SMESH_MeshAlgos.cxx + SMESH_MAT2d.cxx ) # --- rules --- diff --git a/src/SMESHUtils/SMESH_MAT2d.cxx b/src/SMESHUtils/SMESH_MAT2d.cxx new file mode 100644 index 000000000..390d337b8 --- /dev/null +++ b/src/SMESHUtils/SMESH_MAT2d.cxx @@ -0,0 +1,1571 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : SMESH_MAT2d.cxx +// Created : Thu May 28 17:49:53 2015 +// Author : Edward AGAPOV (eap) + +#include "SMESH_MAT2d.hxx" + +#include + +#include +#include +#include +#include +#include +#include +//#include +#include +// #include +// #include +#include +//#include +#include +#include +#include +#include +#include +#include + +#ifdef _DEBUG_ +#include "SMESH_File.hxx" +#include "SMESH_Comment.hxx" +#endif + +using namespace std; +using boost::polygon::x; +using boost::polygon::y; +using SMESH_MAT2d::TVD; +using SMESH_MAT2d::TVDEdge; +using SMESH_MAT2d::TVDCell; +using SMESH_MAT2d::TVDVertex; + +namespace +{ + // Input data for construct_voronoi() + // ------------------------------------------------------------------------------------- + + struct InPoint + { + int _a, _b; + double _param; + InPoint(int x, int y, double param) : _a(x), _b(y), _param(param) {} + InPoint() : _a(0), _b(0), _param(0) {} + + // working data + list< const TVDEdge* > _edges; // MA edges of a concave InPoint in CCW order + + size_t index( const vector< InPoint >& inPoints ) const { return this - &inPoints[0]; } + bool operator==( const InPoint& other ) const { return _a == other._a && _b == other._b; } + }; + // ------------------------------------------------------------------------------------- + + struct InSegment + { + InPoint * _p0; + InPoint * _p1; + + // working data + size_t _geomEdgeInd; // EDGE index within the FACE + const TVDCell* _cell; + list< const TVDEdge* > _edges; // MA edges in CCW order within _cell + + InSegment( InPoint * p0, InPoint * p1, size_t iE) + : _p0(p0), _p1(p1), _geomEdgeInd(iE) {} + InSegment() : _p0(0), _p1(0), _geomEdgeInd(0) {} + + inline bool isConnected( const TVDEdge* edge ); + + inline bool isExternal( const TVDEdge* edge ); + + static void setGeomEdgeToCell( const TVDCell* cell, size_t eID ) { cell->color( eID ); } + + static size_t getGeomEdge( const TVDCell* cell ) { return cell->color(); } + }; + + // check if a TVDEdge begins at my end or ends at my start + inline bool InSegment::isConnected( const TVDEdge* edge ) + { + return ((Abs( edge->vertex0()->x() - _p1->_a ) < 1.&& + Abs( edge->vertex0()->y() - _p1->_b ) < 1. ) || + (Abs( edge->vertex1()->x() - _p0->_a ) < 1.&& + Abs( edge->vertex1()->y() - _p0->_b ) < 1. )); + } + + // check if a MA TVDEdge is outside of a domain + inline bool InSegment::isExternal( const TVDEdge* edge ) + { + double dot = // x1*x2 + y1*y2; (x1,y1) - internal normal of InSegment + ( _p0->_b - _p1->_b ) * ( 0.5 * ( edge->vertex0()->x() + edge->vertex1()->x() ) - _p0->_a ) + + ( _p1->_a - _p0->_a ) * ( 0.5 * ( edge->vertex0()->y() + edge->vertex1()->y() ) - _p0->_b ); + return dot < 0.; + } + + // // ------------------------------------------------------------------------------------- + // const size_t theExternMA = 111; // to mark external MA edges + + // bool isExternal( const TVDEdge* edge ) + // { + // return ( SMESH_MAT2d::Branch::getBndSegment( edge ) == theExternMA ); + // } + + // // mark external MA edges + // void markExternalEdges( const TVDEdge* edge ) + // { + // if ( isExternal( edge )) + // return; + // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge ); + // SMESH_MAT2d::Branch::setBndSegment( theExternMA, edge->twin() ); + // if ( edge->is_primary() && edge->vertex1() ) + // { + // const TVDVertex * v = edge->vertex1(); + // edge = v->incident_edge(); + // do { + // markExternalEdges( edge ); + // edge = edge->rot_next(); + // } while ( edge != v->incident_edge() ); + // } + // } + + // ------------------------------------------------------------------------------------- +#ifdef _DEBUG_ + // writes segments into a txt file readable by voronoi_visualizer + void inSegmentsToFile( vector< InSegment>& inSegments) + { + if ( inSegments.size() > 1000 ) + return; + const char* fileName = "/misc/dn25/salome/eap/salome/misc/Code/C++/MAdebug.txt"; + SMESH_File file(fileName, false ); + file.openForWriting(); + SMESH_Comment text; + text << "0\n"; // nb points + text << inSegments.size() << "\n"; // nb segments + for ( size_t i = 0; i < inSegments.size(); ++i ) + { + text << inSegments[i]._p0->_a << " " + << inSegments[i]._p0->_b << " " + << inSegments[i]._p1->_a << " " + << inSegments[i]._p1->_b << "\n"; + } + text << "\n"; + file.write( text.c_str(), text.size() ); + cout << "Write " << fileName << endl; + } + void dumpEdge( const TVDEdge* edge ) + { + cout << "*Edge_" << edge; + if ( !edge->vertex0() ) + cout << " ( INF, INF"; + else + cout << " ( " << edge->vertex0()->x() << ", " << edge->vertex0()->y(); + if ( !edge->vertex1() ) + cout << ") -> ( INF, INF"; + else + cout << ") -> (" << edge->vertex1()->x() << ", " << edge->vertex1()->y(); + cout << ")\t cell=" << edge->cell() + << " iBnd=" << edge->color() + << " twin=" << edge->twin() + << " twin_cell=" << edge->twin()->cell() + << " prev=" << edge->prev() << " next=" << edge->next() + << ( edge->is_primary() ? " MA " : " SCND" ) + << ( edge->is_linear() ? " LIN " : " CURV" ) + << endl; + } + void dumpCell( const TVDCell* cell ) + { + cout << "**Cell_" << cell << " GEOM=" << cell->color() << " "; + cout << ( cell->contains_segment() ? " SEG " : " PNT " ); + if ( cell-> is_degenerate() ) + cout << " degen "; + else + { + cout << endl; + const TVDEdge* edge = cell->incident_edge(); + size_t i = 0; + do { + edge = edge->next(); + cout << " - " << ++i << " "; + dumpEdge( edge ); + } while (edge != cell->incident_edge()); + } + } +#else + void inSegmentsToFile( vector< InSegment>& inSegments) {} + void dumpEdge( const TVDedge* edge ) {} + void dumpCell( const TVDCell* cell ) {} +#endif +} +// ------------------------------------------------------------------------------------- + +namespace boost { + namespace polygon { + + template <> + struct geometry_concept { + typedef point_concept type; + }; + template <> + struct point_traits { + typedef int coordinate_type; + + static inline coordinate_type get(const InPoint& point, orientation_2d orient) { + return (orient == HORIZONTAL) ? point._a : point._b; + } + }; + + template <> + struct geometry_concept { + typedef segment_concept type; + }; + + template <> + struct segment_traits { + typedef int coordinate_type; + typedef InPoint point_type; + + static inline point_type get(const InSegment& segment, direction_1d dir) { + return *(dir.to_int() ? segment._p1 : segment._p0); + } + }; + } // namespace polygon +} // namespace boost + // ------------------------------------------------------------------------------------- + +namespace +{ + const int theNoBrachID = 0; // std::numeric_limits::max(); + + // ------------------------------------------------------------------------------------- + /*! + * \brief Intermediate DS to create InPoint's + */ + struct UVU + { + gp_Pnt2d _uv; + double _u; + UVU( gp_Pnt2d uv, double u ): _uv(uv), _u(u) {} + InPoint getInPoint( double scale[2] ) + { + return InPoint( int( _uv.X() * scale[0]), int( _uv.Y() * scale[1]), _u ); + } + }; + // ------------------------------------------------------------------------------------- + /*! + * \brief A segment on EDGE, used to create BndPoints + */ + struct BndSeg + { + InSegment* _inSeg; + const TVDEdge* _edge; + double _uLast; + int _branchID; // negative ID means reverse direction + + BndSeg( InSegment* seg, const TVDEdge* edge, double u ): + _inSeg(seg), _edge(edge), _uLast(u), _branchID( theNoBrachID ) {} + + void setIndexToEdge( size_t id ) + { + SMESH_MAT2d::Branch::setBndSegment( id, _edge ); + } + + int branchID() const { return Abs( _branchID ); } + + size_t geomEdge() const { return _inSeg->_geomEdgeInd; } + + void setBranch( int branchID, vector< BndSeg >& bndSegs ) + { + _branchID = branchID; + + if ( _edge ) // pass branch to an opposite BndSeg + { + size_t oppSegIndex = SMESH_MAT2d::Branch::getBndSegment( _edge->twin() ); + if ( oppSegIndex < bndSegs.size() /*&& bndSegs[ oppSegIndex ]._branchID == theNoBrachID*/ ) + bndSegs[ oppSegIndex ]._branchID = -branchID; + } + } + bool hasOppositeEdge( const size_t noEdgeID ) + { + if ( !_edge ) return false; + return ( _inSeg->getGeomEdge( _edge->twin()->cell() ) != noEdgeID ); + } + + // check a next segment in CW order + bool isSameBranch( const BndSeg& seg2 ) + { + if ( !_edge || !seg2._edge ) + return true; + + const TVDCell* cell1 = this->_edge->twin()->cell(); + const TVDCell* cell2 = seg2. _edge->twin()->cell(); + if ( cell1 == cell2 ) + return true; + + const TVDEdge* edgeMedium1 = this->_edge->twin()->next(); + const TVDEdge* edgeMedium2 = seg2. _edge->twin()->prev(); + + if ( edgeMedium1->is_secondary() && edgeMedium2->is_secondary() ) + { + if ( edgeMedium1->twin() == edgeMedium2 ) + return true; + // edgeMedium's are edges whose twin()->cell is built on an end point of inSegment + // and is located between cell1 and cell2 + if ( edgeMedium1->twin() == edgeMedium2->twin() ) // is this possible??? + return true; + if ( edgeMedium1->twin() == edgeMedium2->twin()->next() && + edgeMedium1->twin()->cell()->contains_point() ) + return true; + } + else if ( edgeMedium1->is_primary() && edgeMedium2->is_primary() ) + { + if ( edgeMedium1->twin() == edgeMedium2 && + SMESH_MAT2d::Branch::getBndSegment( edgeMedium1 ) == + SMESH_MAT2d::Branch::getBndSegment( edgeMedium2 )) + // this is an ignored MA edge between inSegment's on one EDGE forming a convex corner + return true; + } + + return false; + } + }; + + //================================================================================ + /*! + * \brief Computes length of of TVDEdge + */ + //================================================================================ + + double length( const TVDEdge* edge ) + { + gp_XY d( edge->vertex0()->x() - edge->vertex1()->x(), + edge->vertex0()->y() - edge->vertex1()->y() ); + return d.Modulus(); + } + + //================================================================================ + /*! + * \brief Compute scale to have the same 2d proportions as in 3d + */ + //================================================================================ + + void computeProportionScale( const TopoDS_Face& face, + const Bnd_B2d& uvBox, + double scale[2]) + { + scale[0] = scale[1] = 1.; + if ( uvBox.IsVoid() ) return; + + TopLoc_Location loc; + Handle(Geom_Surface) surface = BRep_Tool::Surface( face, loc ); + + const int nbDiv = 30; + gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax(); + gp_XY uvMid = 0.5 * ( uvMin + uvMax ); + double du = ( uvMax.X() - uvMin.X() ) / nbDiv; + double dv = ( uvMax.Y() - uvMin.Y() ) / nbDiv; + + double uLen3d = 0, vLen3d = 0; + gp_Pnt uPrevP = surface->Value( uvMin.X(), uvMid.Y() ); + gp_Pnt vPrevP = surface->Value( uvMid.X(), uvMin.Y() ); + for (int i = 1; i <= nbDiv; i++) + { + double u = uvMin.X() + du * i; + double v = uvMin.Y() + dv * i; + gp_Pnt uP = surface->Value( u, uvMid.Y() ); + gp_Pnt vP = surface->Value( uvMid.X(), v ); + uLen3d += uP.Distance( uPrevP ); + vLen3d += vP.Distance( vPrevP ); + uPrevP = uP; + vPrevP = vP; + } + scale[0] = uLen3d / ( uvMax.X() - uvMin.X() ); + scale[1] = vLen3d / ( uvMax.Y() - uvMin.Y() ); + } + + //================================================================================ + /*! + * \brief Fill input data for construct_voronoi() + */ + //================================================================================ + + bool makeInputData(const TopoDS_Face& face, + const std::vector< TopoDS_Edge >& edges, + const double minSegLen, + vector< InPoint >& inPoints, + vector< InSegment>& inSegments, + double scale[2]) + { + const double theDiscrCoef = 0.5; // to decrease minSegLen for discretization + TopLoc_Location loc; + + // discretize the EDGEs to get 2d points and segments + + vector< vector< UVU > > uvuVec( edges.size() ); + Bnd_B2d uvBox; + for ( size_t iE = 0; iE < edges.size(); ++iE ) + { + vector< UVU > & points = uvuVec[ iE ]; + + double f,l; + Handle(Geom_Curve) c3d = BRep_Tool::Curve ( edges[ iE ], loc, f, l ); + Handle(Geom2d_Curve) c2d = BRep_Tool::CurveOnSurface( edges[ iE ], face, f, l ); + if ( c2d.IsNull() ) return false; + + points.push_back( UVU( c2d->Value( f ), f )); + uvBox.Add( points.back()._uv ); + + Geom2dAdaptor_Curve c2dAdaptor (c2d, f,l ); + double curDeflect = 0.3; //0.01; //Curvature deflection + double angDeflect = 0.2; // 0.09; //Angular deflection + + GCPnts_TangentialDeflection discret(c2dAdaptor, angDeflect, curDeflect); + // if ( discret.NbPoints() > 2 ) + // { + // cout << endl; + // do + // { + // discret.Initialize( c2dAdaptor, 100, curDeflect ); + // cout << "C " << curDeflect << " " << discret.NbPoints() << endl; + // curDeflect *= 1.5; + // } + // while ( discret.NbPoints() > 5 ); + // cout << endl; + // do + // { + // discret.Initialize( c2dAdaptor, angDeflect, 100 ); + // cout << "A " << angDeflect << " " << discret.NbPoints() << endl; + // angDeflect *= 1.5; + // } + // while ( discret.NbPoints() > 5 ); + // } + gp_Pnt p, pPrev; + if ( !c3d.IsNull() ) + pPrev = c3d->Value( f ); + for ( int i = 2; i <= discret.NbPoints(); i++ ) // skip the 1st point + { + double u = discret.Parameter(i); + if ( !c3d.IsNull() ) + { + p = c3d->Value( u ); + int nbDiv = int( p.Distance( pPrev ) / minSegLen / theDiscrCoef ); + double dU = ( u - points.back()._u ) / nbDiv; + for ( int iD = 1; iD < nbDiv; ++iD ) + { + double uD = points.back()._u + dU; + points.push_back( UVU( c2d->Value( uD ), uD )); + } + pPrev = p; + } + points.push_back( UVU( c2d->Value( u ), u )); + uvBox.Add( points.back()._uv ); + } + // if ( !c3d.IsNull() ) + // { + // vector params; + // GeomAdaptor_Curve c3dAdaptor( c3d,f,l ); + // if ( useDefl ) + // { + // const double deflection = minSegLen * 0.1; + // GCPnts_UniformDeflection discret( c3dAdaptor, deflection, f, l, true ); + // if ( !discret.IsDone() ) + // return false; + // int nbP = discret.NbPoints(); + // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points + // params.push_back( discret.Parameter(i) ); + // } + // else + // { + // double eLen = GCPnts_AbscissaPoint::Length( c3dAdaptor ); + // int nbSeg = Max( 1, int( eLen / minSegLen / theDiscrCoef )); + // double segLen = eLen / nbSeg; + // GCPnts_UniformAbscissa discret( c3dAdaptor, segLen, f, l ); + // int nbP = Min( discret.NbPoints(), nbSeg + 1 ); + // for ( int i = 2; i < nbP; i++ ) // skip 1st and last points + // params.push_back( discret.Parameter(i) ); + // } + // for ( size_t i = 0; i < params.size(); ++i ) + // { + // points.push_back( UVU( c2d->Value( params[i] ), params[i] )); + // uvBox.Add( points.back()._uv ); + // } + // } + if ( points.size() < 2 ) + { + points.push_back( UVU( c2d->Value( l ), l )); + uvBox.Add( points.back()._uv ); + } + if ( edges[ iE ].Orientation() == TopAbs_REVERSED ) + std::reverse( points.begin(), points.end() ); + } + + // make connected EDGEs have same UV at shared VERTEX + TopoDS_Vertex vShared; + for ( size_t iE = 0; iE < edges.size(); ++iE ) + { + size_t iE2 = (iE+1) % edges.size(); + if ( !TopExp::CommonVertex( edges[iE], edges[iE2], vShared )) + continue; + if ( !vShared.IsSame( TopExp::LastVertex( edges[iE], true ))) + return false; + vector< UVU > & points1 = uvuVec[ iE ]; + vector< UVU > & points2 = uvuVec[ iE2 ]; + gp_Pnt2d & uv1 = points1.back() ._uv; + gp_Pnt2d & uv2 = points2.front()._uv; + uv1 = uv2 = 0.5 * ( uv1.XY() + uv2.XY() ); + } + + // get scale to have the same 2d proportions as in 3d + computeProportionScale( face, uvBox, scale ); + + // make scale to have coordinates precise enough when converted to int + + gp_XY uvMin = uvBox.CornerMin(), uvMax = uvBox.CornerMax(); + uvMin.ChangeCoord(1) = uvMin.X() * scale[0]; + uvMin.ChangeCoord(2) = uvMin.Y() * scale[1]; + uvMax.ChangeCoord(1) = uvMax.X() * scale[0]; + uvMax.ChangeCoord(2) = uvMax.Y() * scale[1]; + double vMax[2] = { Max( Abs( uvMin.X() ), Abs( uvMax.X() )), + Max( Abs( uvMin.Y() ), Abs( uvMax.Y() )) }; + int iMax = ( vMax[0] > vMax[1] ) ? 0 : 1; + const double precision = 1e-5; + double preciScale = Min( vMax[iMax] / precision, + std::numeric_limits::max() / vMax[iMax] ); + preciScale /= scale[iMax]; + double roundedScale = 10; // to ease debug + while ( roundedScale * 10 < preciScale ) + roundedScale *= 10.; + scale[0] *= roundedScale; + scale[1] *= roundedScale; + + // create input points and segments + + inPoints.clear(); + inSegments.clear(); + size_t nbPnt = 0; + for ( size_t iE = 0; iE < uvuVec.size(); ++iE ) + nbPnt += uvuVec[ iE ].size(); + inPoints.resize( nbPnt ); + inSegments.reserve( nbPnt ); + + size_t iP = 0; + if ( face.Orientation() == TopAbs_REVERSED ) + { + for ( int iE = uvuVec.size()-1; iE >= 0; --iE ) + { + vector< UVU > & points = uvuVec[ iE ]; + inPoints[ iP++ ] = points.back().getInPoint( scale ); + for ( size_t i = points.size()-1; i >= 1; --i ) + { + inPoints[ iP++ ] = points[i-1].getInPoint( scale ); + inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE )); + } + } + } + else + { + for ( size_t iE = 0; iE < uvuVec.size(); ++iE ) + { + vector< UVU > & points = uvuVec[ iE ]; + inPoints[ iP++ ] = points[0].getInPoint( scale ); + for ( size_t i = 1; i < points.size(); ++i ) + { + inPoints[ iP++ ] = points[i].getInPoint( scale ); + inSegments.push_back( InSegment( & inPoints[ iP-2 ], & inPoints[ iP-1 ], iE )); + } + } + } + return true; + } + + //================================================================================ + /*! + * \brief Create MA branches and FACE boundary data + * \param [in] vd - voronoi diagram of \a inSegments + * \param [in] inPoints - FACE boundary points + * \param [in,out] inSegments - FACE boundary segments + * \param [out] branch - MA branches to fill + * \param [out] branchEnd - ends of MA branches to fill + * \param [out] boundary - FACE boundary to fill + */ + //================================================================================ + + void makeMA( const TVD& vd, + vector< InPoint >& inPoints, + vector< InSegment > & inSegments, + vector< SMESH_MAT2d::Branch >& branch, + vector< const SMESH_MAT2d::BranchEnd* >& branchPnt, + SMESH_MAT2d::Boundary& boundary ) + { + const size_t noEdgeID = inSegments.size() + 1; // ID of non-existent geom EDGE + + // Associate MA cells with inSegments + for (TVD::const_cell_iterator it = vd.cells().begin(); it != vd.cells().end(); ++it) + { + const TVDCell* cell = &(*it); + if ( cell->contains_segment() ) + { + InSegment& seg = inSegments[ cell->source_index() ]; + seg._cell = cell; + seg.setGeomEdgeToCell( cell, seg._geomEdgeInd ); + } + else + { + InSegment::setGeomEdgeToCell( cell, noEdgeID ); + } + } + + vector< bool > inPntChecked( inPoints.size(), false ); + + // Find MA edges of each inSegment + + for ( size_t i = 0; i < inSegments.size(); ++i ) + { + InSegment& inSeg = inSegments[i]; + + // get edges around the cell lying on MA + bool hasSecondary = false; + const TVDEdge* edge = inSeg._cell->incident_edge(); + do { + edge = edge->next(); // Returns the CCW next edge within the cell. + if ( edge->is_primary() && !inSeg.isExternal( edge ) ) + inSeg._edges.push_back( edge ); // edge equidistant from two InSegments + else + hasSecondary = true; + } while (edge != inSeg._cell->incident_edge()); + + // there can be several continuous MA edges but maEdges can begin in the middle of + // a chain of continuous MA edges. Make the chain continuous. + list< const TVDEdge* >& maEdges = inSeg._edges; + if ( maEdges.empty() ) + continue; + if ( hasSecondary ) + while ( maEdges.back()->next() == maEdges.front() ) + maEdges.splice( maEdges.end(), maEdges, maEdges.begin() ); + + // remove maEdges equidistant from two neighbor InSegments of the same geom EDGE + list< const TVDEdge* >::iterator e = maEdges.begin(); + while ( e != maEdges.end() ) + { + const TVDCell* cell2 = (*e)->twin()->cell(); // cell on the other side of a MA edge + size_t geoE2 = InSegment::getGeomEdge( cell2 ); + bool toRemove = ( inSeg._geomEdgeInd == geoE2 && inSeg.isConnected( *e )); + if ( toRemove ) + e = maEdges.erase( e ); + else + ++e; + } + if ( maEdges.empty() ) + continue; + + // add MA edges corresponding to concave InPoints + for ( int is2nd = 0; is2nd < 2; ++is2nd ) // loop on two ends of inSeg + { + InPoint& inPnt = *( is2nd ? inSeg._p1 : inSeg._p0 ); + size_t pInd = inPnt.index( inPoints ); + if ( inPntChecked[ pInd ] ) + continue; + if ( pInd > 0 && + inPntChecked[ pInd-1 ] && + inPoints[ pInd-1 ] == inPnt ) + continue; + inPntChecked[ pInd ] = true; + + const TVDEdge* edge = // a TVDEdge passing through an end of inSeg + is2nd ? maEdges.front()->prev() : maEdges.back()->next(); + while ( true ) + { + if ( edge->is_primary() ) break; // this should not happen + const TVDEdge* edge2 = edge->twin(); // we are in a neighbor cell, add MA edges to inPnt + if ( inSeg.getGeomEdge( edge2->cell() ) != noEdgeID ) + break; // cell of an InSegment + bool hasInfinite = false; + list< const TVDEdge* > pointEdges; + edge = edge2; + do + { + edge = edge->next(); // Returns the CCW next edge within the cell. + if ( edge->is_infinite() ) + hasInfinite = true; + else if ( edge->is_primary() && !inSeg.isExternal( edge )) + pointEdges.push_back( edge ); + } + while ( edge != edge2 && !hasInfinite ); + + if ( hasInfinite || pointEdges.empty() ) + break; + inPnt._edges.splice( inPnt._edges.end(), pointEdges ); + inSeg.setGeomEdgeToCell( edge->cell(), inSeg._geomEdgeInd ); + + edge = is2nd ? inPnt._edges.front()->prev() : inPnt._edges.back()->next(); + } + } // add MA edges corresponding to concave InPoints + + } // loop on inSegments to find corresponding MA edges + + + // ------------------------------------------- + // Create Branches and BndPoints for each EDGE + // ------------------------------------------- + + if ( inPoints.front() == inPoints.back() /*&& !inPoints[0]._edges.empty()*/ ) + { + inPntChecked[0] = false; // do not use the 1st point twice + //InSegment::setGeomEdgeToCell( inPoints[0]._edges.back()->cell(), noEdgeID ); + inPoints[0]._edges.clear(); + } + + // Divide InSegment's into BndSeg's + + vector< BndSeg > bndSegs; + bndSegs.reserve( inSegments.size() * 3 ); + + list< const TVDEdge* >::reverse_iterator e; + for ( size_t i = 0; i < inSegments.size(); ++i ) + { + InSegment& inSeg = inSegments[i]; + + // segments around 1st concave point + size_t ip0 = inSeg._p0->index( inPoints ); + if ( inPntChecked[ ip0 ] ) + for ( e = inSeg._p0->_edges.rbegin(); e != inSeg._p0->_edges.rend(); ++e ) + bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p0->_param )); + inPntChecked[ ip0 ] = false; + + // segments of InSegment's + size_t nbMaEdges = inSeg._edges.size(); + switch ( nbMaEdges ) { + case 0: // "around" circle center + bndSegs.push_back( BndSeg( &inSeg, 0, inSeg._p1->_param )); break; + case 1: + bndSegs.push_back( BndSeg( &inSeg, inSeg._edges.back(), inSeg._p1->_param )); break; + default: + vector< double > len; + len.push_back(0); + for ( e = inSeg._edges.rbegin(); e != inSeg._edges.rend(); ++e ) + len.push_back( len.back() + length( *e )); + + e = inSeg._edges.rbegin(); + for ( size_t l = 1; l < len.size(); ++e, ++l ) + { + double dl = len[l] / len.back(); + double u = dl * inSeg._p1->_param + ( 1. - dl ) * inSeg._p0->_param; + bndSegs.push_back( BndSeg( &inSeg, *e, u )); + } + } + // segments around 2nd concave point + size_t ip1 = inSeg._p1->index( inPoints ); + if ( inPntChecked[ ip1 ] ) + for ( e = inSeg._p1->_edges.rbegin(); e != inSeg._p1->_edges.rend(); ++e ) + bndSegs.push_back( BndSeg( &inSeg, *e, inSeg._p1->_param )); + inPntChecked[ ip1 ] = false; + } + + // make TVDEdge's know it's BndSeg to enable passing branchID to + // an opposite BndSeg in BndSeg::setBranch() + for ( size_t i = 0; i < bndSegs.size(); ++i ) + bndSegs[i].setIndexToEdge( i ); + + + // Find TVDEdge's of Branches and associate them with bndSegs + + vector< vector > branchEdges; + branchEdges.reserve( boundary.nbEdges() * 4 ); + + map< const TVDVertex*, SMESH_MAT2d::BranchEndType > endType; + + int branchID = 1; // we code orientation as branchID sign + branchEdges.resize( branchID + 1 ); + + size_t i1st = 0; + while ( i1st < bndSegs.size() && !bndSegs[i1st].hasOppositeEdge( noEdgeID )) + ++i1st; + bndSegs[i1st].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg + branchEdges[ branchID ].push_back( bndSegs[i1st]._edge ); + + for ( size_t i = i1st+1; i < bndSegs.size(); ++i ) + { + if ( bndSegs[i].branchID() ) + { + branchID = bndSegs[i]._branchID; // with sign + + if ( bndSegs[i]._branchID == -bndSegs[i-1]._branchID && + bndSegs[i]._edge ) + { + SMESH_MAT2d::BranchEndType type = + ( bndSegs[i]._inSeg->isConnected( bndSegs[i]._edge ) ? + SMESH_MAT2d::BE_ON_VERTEX : + SMESH_MAT2d::BE_END ); + endType.insert( make_pair( bndSegs[i]._edge->vertex1(), type )); + } + continue; + } + if ( !bndSegs[i-1].isSameBranch( bndSegs[i] )) + { + branchEdges.resize(( branchID = branchEdges.size()) + 1 ); + if ( bndSegs[i]._edge ) + endType.insert( make_pair( bndSegs[i]._edge->vertex1(), + SMESH_MAT2d::BE_BRANCH_POINT )); + } + bndSegs[i].setBranch( branchID, bndSegs ); // set to i-th and the opposite bndSeg + if ( bndSegs[i].hasOppositeEdge( noEdgeID )) + branchEdges[ bndSegs[i].branchID() ].push_back( bndSegs[i]._edge ); + } + // define BranchEndType of the first TVDVertex + if ( bndSegs.front()._branchID == -bndSegs.back()._branchID ) + { + if ( bndSegs[0]._edge ) + { + SMESH_MAT2d::BranchEndType type = + ( bndSegs[0]._inSeg->isConnected( bndSegs[0]._edge ) ? + SMESH_MAT2d::BE_ON_VERTEX : + SMESH_MAT2d::BE_END ); + endType.insert( make_pair( bndSegs[0]._edge->vertex1(), type )); + } + else if ( bndSegs.back()._edge ) + { + SMESH_MAT2d::BranchEndType type = + ( bndSegs.back()._inSeg->isConnected( bndSegs.back()._edge ) ? + SMESH_MAT2d::BE_ON_VERTEX : + SMESH_MAT2d::BE_END ); + endType.insert( make_pair( bndSegs.back()._edge->vertex0(), type )); + } + } + // join the 1st and the last branch edges if it is the same branch + if ( bndSegs.back().branchID() != bndSegs.front().branchID() && + bndSegs.back().isSameBranch( bndSegs.front() )) + { + vector & br1 = branchEdges[ bndSegs.front().branchID() ]; + vector & br2 = branchEdges[ bndSegs.back().branchID() ]; + br1.insert( br1.begin(), br2.begin(), br2.end() ); + br2.clear(); + } + + // associate branchIDs and the input branch vector (arg) + vector< const SMESH_MAT2d::Branch* > branchByID( branchEdges.size() ); + int nbBranches = 0; + for ( size_t i = 0; i < branchEdges.size(); ++i ) + { + nbBranches += ( !branchEdges[i].empty() ); + } + branch.resize( nbBranches ); + for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID ) + { + if ( !branchEdges[ brID ].empty() ) + branchByID[ brID ] = & branch[ iBr++ ]; + } + + // Fill in BndPoints of each EDGE of the boundary + + size_t iSeg = 0; + int edgeInd = -1, dInd = 0; + while ( iSeg < bndSegs.size() ) + { + const size_t geomID = bndSegs[ iSeg ].geomEdge(); + SMESH_MAT2d::BndPoints & bndPoints = boundary.getPoints( geomID ); + + size_t nbSegs = 0; + for ( size_t i = iSeg; i < bndSegs.size() && geomID == bndSegs[ i ].geomEdge(); ++i ) + ++nbSegs; + size_t iSegEnd = iSeg + nbSegs; + + // make TVDEdge know an index of bndSegs within BndPoints + for ( size_t i = iSeg; i < iSegEnd; ++i ) + if ( bndSegs[i]._edge ) + SMESH_MAT2d::Branch::setBndSegment( i - iSeg, bndSegs[i]._edge ); + + // parameters on EDGE + + bndPoints._params.reserve( nbSegs + 1 ); + bndPoints._params.push_back( bndSegs[ iSeg ]._inSeg->_p0->_param ); + + for ( size_t i = iSeg; i < iSegEnd; ++i ) + bndPoints._params.push_back( bndSegs[ i ]._uLast ); + + // MA edges + + bndPoints._maEdges.reserve( nbSegs ); + + for ( size_t i = iSeg; i < iSegEnd; ++i ) + { + const size_t brID = bndSegs[ i ].branchID(); + const SMESH_MAT2d::Branch* br = branchByID[ brID ]; + + if ( bndSegs[ i ]._edge && !branchEdges[ brID ].empty() ) + { + edgeInd += dInd; + + if ( edgeInd < 0 || + edgeInd >= (int) branchEdges[ brID ].size() || + branchEdges[ brID ][ edgeInd ] != bndSegs[ i ]._edge ) + { + if ( bndSegs[ i ]._branchID < 0 && + branchEdges[ brID ].back() == bndSegs[ i ]._edge ) + { + edgeInd = branchEdges[ brID ].size() - 1; + dInd = -1; + } + else if ( bndSegs[ i ]._branchID > 0 && + branchEdges[ brID ].front() == bndSegs[ i ]._edge ) + { + edgeInd = 0; + dInd = +1; + } + else + { + for ( edgeInd = 0; edgeInd < branchEdges[ brID ].size(); ++edgeInd ) + if ( branchEdges[ brID ][ edgeInd ] == bndSegs[ i ]._edge ) + break; + dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1; + } + } + } + else + { + // no MA edge, bndSeg corresponds to an end point of a branch + if ( bndPoints._maEdges.empty() ) + { + // should not get here according to algo design + edgeInd = 0; + } + else + { + edgeInd = branchEdges[ brID ].size(); + dInd = bndSegs[ i ]._branchID > 0 ? +1 : -1; + } + } + + bndPoints._maEdges.push_back( make_pair( br, ( 1 + edgeInd ) * dInd )); + + } // loop on bndSegs of an EDGE + + iSeg = iSegEnd; + + } // loop on all bndSegs + + + // fill the branches with MA edges + for ( size_t iBr = 0, brID = 0; brID < branchEdges.size(); ++brID ) + if ( !branchEdges[brID].empty() ) + { + branch[ iBr ].init( branchEdges[brID], & boundary, endType ); + iBr++; + } + // set branches to branch ends + for ( size_t i = 0; i < branch.size(); ++i ) + branch[i].setBranchesToEnds( branch ); + + // fill branchPnt arg + map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* > v2end; + for ( size_t i = 0; i < branch.size(); ++i ) + { + if ( branch[i].getEnd(0)->_branches.size() > 2 ) + v2end.insert( make_pair( branch[i].getEnd(0)->_vertex, branch[i].getEnd(0) )); + if ( branch[i].getEnd(1)->_branches.size() > 2 ) + v2end.insert( make_pair( branch[i].getEnd(1)->_vertex, branch[i].getEnd(1) )); + } + branchPnt.resize( v2end.size() ); + map< const TVDVertex*, const SMESH_MAT2d::BranchEnd* >::iterator v2e = v2end.begin(); + for ( size_t i = 0; v2e != v2end.end(); ++v2e, ++i ) + branchPnt[ i ] = v2e->second; + + } // makeMA() + +} // namespace + +//================================================================================ +/*! + * \brief MedialAxis constructor + * \param [in] face - a face to create MA for + * \param [in] edges - edges of the face (possibly not all) on the order they + * encounter in the face boundary. + * \param [in] minSegLen - minimal length of a mesh segment used to discretize + * the edges. It is used to define precision of MA approximation + */ +//================================================================================ + +SMESH_MAT2d::MedialAxis::MedialAxis(const TopoDS_Face& face, + const std::vector< TopoDS_Edge >& edges, + const double minSegLen, + const bool ignoreCorners): + _face( face ), _boundary( edges.size() ) +{ + // input to construct_voronoi() + vector< InPoint > inPoints; + vector< InSegment> inSegments; + if ( !makeInputData( face, edges, minSegLen, inPoints, inSegments, _scale )) + return; + + //inSegmentsToFile( inSegments ); + + // build voronoi diagram + construct_voronoi( inSegments.begin(), inSegments.end(), &_vd ); + + // make MA data + makeMA( _vd, inPoints, inSegments, _branch, _branchPnt, _boundary ); +} + +//================================================================================ +/*! + * \brief Return UVs of ends of MA edges of a branch + */ +//================================================================================ + +void SMESH_MAT2d::MedialAxis::getPoints( const Branch& branch, + std::vector< gp_XY >& points) const +{ + branch.getPoints( points, _scale ); +} + +//================================================================================ +/*! + * \brief Returns a BranchPoint corresponding to a given point on a geom EDGE + * \param [in] iGeomEdge - index of geom EDGE within a vector passed at construction + * \param [in] u - parameter of the point on EDGE curve + * \param [out] p - the found BranchPoint + * \return bool - is OK + */ +//================================================================================ + +bool SMESH_MAT2d::Boundary::getBranchPoint( const std::size_t iEdge, + double u, + BranchPoint& p ) const +{ + if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() ) + return false; + + const BndPoints& points = _pointsPerEdge[ iEdge ]; + const bool edgeReverse = ( points._params[0] > points._params.back() ); + + if ( u < ( edgeReverse ? points._params.back() : points._params[0] )) + u = edgeReverse ? points._params.back() : points._params[0]; + else if ( u > ( edgeReverse ? points._params[0] : points._params.back()) ) + u = edgeReverse ? points._params[0] : points._params.back(); + + double r = ( u - points._params[0] ) / ( points._params.back() - points._params[0] ); + int i = int( r * double( points._maEdges.size()-1 )); + if ( edgeReverse ) + { + while ( points._params[i ] < u ) --i; + while ( points._params[i+1] > u ) ++i; + } + else + { + while ( points._params[i ] > u ) --i; + while ( points._params[i+1] < u ) ++i; + } + double edgeParam = ( u - points._params[i] ) / ( points._params[i+1] - points._params[i] ); + + const std::pair< const Branch*, int >& maE = points._maEdges[ i ]; + bool maReverse = ( maE.second < 0 ); + + p._branch = maE.first; + p._iEdge = maE.second - 1; // countered from 1 to store sign + p._edgeParam = maReverse ? ( 1. - edgeParam ) : edgeParam; + + return true; +} + +//================================================================================ +/*! + * \brief Check if a given boundary segment is a null-length segment on a concave + * boundary corner. + * \param [in] iEdge - index of a geom EDGE + * \param [in] iSeg - index of a boundary segment + * \return bool - true if the segment is on concave corner + */ +//================================================================================ + +bool SMESH_MAT2d::Boundary::IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const +{ + if ( iEdge >= _pointsPerEdge.size() || _pointsPerEdge[iEdge]._params.empty() ) + return false; + + const BndPoints& points = _pointsPerEdge[ iEdge ]; + if ( points._params.size() >= iSeg+1 ) + return false; + + return Abs( points._params[ iEdge ] - points._params[ iEdge+1 ]) < 1e-20; +} + +//================================================================================ +/*! + * \brief Creates a 3d curve corresponding to a Branch + * \param [in] branch - the Branch + * \return Adaptor3d_Curve* - the new curve the caller is to delete + */ +//================================================================================ + +Adaptor3d_Curve* SMESH_MAT2d::MedialAxis::make3DCurve(const Branch& branch) const +{ + Handle(Geom_Surface) surface = BRep_Tool::Surface( _face ); + if ( surface.IsNull() ) + return 0; + + vector< gp_XY > uv; + branch.getPoints( uv, _scale ); + if ( uv.size() < 2 ) + return 0; + + vector< TopoDS_Vertex > vertex( uv.size() ); + for ( size_t i = 0; i < uv.size(); ++i ) + vertex[i] = BRepBuilderAPI_MakeVertex( surface->Value( uv[i].X(), uv[i].Y() )); + + TopoDS_Wire aWire; + BRep_Builder aBuilder; + aBuilder.MakeWire(aWire); + for ( size_t i = 1; i < vertex.size(); ++i ) + { + TopoDS_Edge edge = BRepBuilderAPI_MakeEdge( vertex[i-1], vertex[i] ); + aBuilder.Add( aWire, edge ); + } + + // if ( myEdge.size() == 2 && FirstVertex().IsSame( LastVertex() )) + // aWire.Closed(true); // issue 0021141 + + return new BRepAdaptor_CompCurve( aWire ); +} + +//================================================================================ +/*! + * \brief Copy points of an EDGE + */ +//================================================================================ + +void SMESH_MAT2d::Branch::init( vector& maEdges, + const Boundary* boundary, + map< const TVDVertex*, BranchEndType > endType ) +{ + if ( maEdges.empty() ) return; + + _boundary = boundary; + _maEdges.swap( maEdges ); + + + _params.reserve( _maEdges.size() + 1 ); + _params.push_back( 0. ); + for ( size_t i = 0; i < _maEdges.size(); ++i ) + _params.push_back( _params.back() + length( _maEdges[i] )); + + for ( size_t i = 1; i < _params.size(); ++i ) + _params[i] /= _params.back(); + + + _endPoint1._vertex = _maEdges.front()->vertex1(); + _endPoint2._vertex = _maEdges.back ()->vertex0(); + + if ( endType.count( _endPoint1._vertex )) + _endPoint1._type = endType[ _endPoint1._vertex ]; + if ( endType.count( _endPoint2._vertex )) + _endPoint2._type = endType[ _endPoint2._vertex ]; +} + +//================================================================================ +/*! + * \brief fill BranchEnd::_branches of its ends + */ +//================================================================================ + +void SMESH_MAT2d::Branch::setBranchesToEnds( const vector< Branch >& branches ) +{ + for ( size_t i = 0; i < branches.size(); ++i ) + { + if ( this->_endPoint1._vertex == branches[i]._endPoint1._vertex || + this->_endPoint1._vertex == branches[i]._endPoint2._vertex ) + this->_endPoint1._branches.push_back( &branches[i] ); + + if ( this->_endPoint2._vertex == branches[i]._endPoint1._vertex || + this->_endPoint2._vertex == branches[i]._endPoint2._vertex ) + this->_endPoint2._branches.push_back( &branches[i] ); + } +} + +//================================================================================ +/*! + * \brief Returns points on two EDGEs, equidistant from a given point of this Branch + * \param [in] param - [0;1] normalized param on the Branch + * \param [out] bp1 - BoundaryPoint on EDGE with a lower index + * \param [out] bp2 - BoundaryPoint on EDGE with a higher index + * \return bool - true if the BoundaryPoint's found + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::getBoundaryPoints(double param, + BoundaryPoint& bp1, + BoundaryPoint& bp2 ) const +{ + if ( param < _params[0] || param > _params.back() ) + return false; + + // look for an index of a MA edge by param + double ip = param * _params.size(); + size_t i = size_t( Min( int( _maEdges.size()-1), int( ip ))); + + while ( param < _params[i ] ) --i; + while ( param > _params[i+1] ) ++i; + + double r = ( param - _params[i] ) / ( _params[i+1] - _params[i] ); + + return getBoundaryPoints( i, r, bp1, bp2 ); +} + +//================================================================================ +/*! + * \brief Returns points on two EDGEs, equidistant from a given point of this Branch + * \param [in] iMAEdge - index of a MA edge within this Branch + * \param [in] maEdgeParam - [0;1] normalized param on the \a iMAEdge + * \param [out] bp1 - BoundaryPoint on EDGE with a lower index + * \param [out] bp2 - BoundaryPoint on EDGE with a higher index + * \return bool - true if the BoundaryPoint's found + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::getBoundaryPoints(std::size_t iMAEdge, + double maEdgeParam, + BoundaryPoint& bp1, + BoundaryPoint& bp2 ) const +{ + if ( iMAEdge > _maEdges.size() ) + return false; + + size_t iGeom1 = getGeomEdge( _maEdges[ iMAEdge ] ); + size_t iGeom2 = getGeomEdge( _maEdges[ iMAEdge ]->twin() ); + size_t iSeg1 = getBndSegment( _maEdges[ iMAEdge ] ); + size_t iSeg2 = getBndSegment( _maEdges[ iMAEdge ]->twin() ); + + return ( _boundary->getPoint( iGeom1, iSeg1, maEdgeParam, bp1 ) && + _boundary->getPoint( iGeom2, iSeg2, maEdgeParam, bp2 )); +} + +//================================================================================ +/*! + * \brief Returns points on two EDGEs, equidistant from a given point of this Branch + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::getBoundaryPoints(const BranchPoint& p, + BoundaryPoint& bp1, + BoundaryPoint& bp2 ) const +{ + return ( p._branch ? p._branch : this )->getBoundaryPoints( p._iEdge, p._edgeParam, bp1, bp2 ); +} + +//================================================================================ +/*! + * \brief Return a parameter of a BranchPoint normalized within this Branch + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::getParameter(const BranchPoint & p, double & u ) const +{ + if ( p._iEdge > _params.size()-1 ) + return false; + + u = ( _params[ p._iEdge ] * ( 1 - p._edgeParam ) + + _params[ p._iEdge+1 ] * p._edgeParam ); + + return true; +} + +//================================================================================ +/*! + * \brief Check type of both ends + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::hasEndOfType(BranchEndType type) const +{ + return ( _endPoint1._type == type || _endPoint2._type == type ); +} + +//================================================================================ +/*! + * \brief Returns MA points + * \param [out] points - the 2d points + * \param [in] scale - the scale that was used to scale the 2d space of MA + */ +//================================================================================ + +void SMESH_MAT2d::Branch::getPoints( std::vector< gp_XY >& points, + const double scale[2]) const +{ + points.resize( _maEdges.size() + 1 ); + + points[0].SetCoord( _maEdges[0]->vertex1()->x() / scale[0], // CCW order! -> vertex1 not vertex0 + _maEdges[0]->vertex1()->y() / scale[1] ); + + for ( size_t i = 0; i < _maEdges.size(); ++i ) + points[i+1].SetCoord( _maEdges[i]->vertex0()->x() / scale[0], + _maEdges[i]->vertex0()->y() / scale[1] ); +} + +//================================================================================ +/*! + * \brief Return indices of EDGEs equidistant from this branch + */ +//================================================================================ + +void SMESH_MAT2d::Branch::getGeomEdges( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2 ) const +{ + edgeIDs1.push_back( getGeomEdge( _maEdges[0] )); + edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() )); + + for ( size_t i = 1; i < _maEdges.size(); ++i ) + { + size_t ie1 = getGeomEdge( _maEdges[i] ); + size_t ie2 = getGeomEdge( _maEdges[i]->twin() ); + + if ( edgeIDs1.back() != ie1 ) edgeIDs1.push_back( ie1 ); + if ( edgeIDs2.back() != ie2 ) edgeIDs2.push_back( ie2 ); + } +} + +//================================================================================ +/*! + * \brief Looks for a BranchPoint position around a concave VERTEX + */ +//================================================================================ + +bool SMESH_MAT2d::Branch::addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2, + std::vector< BranchPoint >& divPoints, + const vector& maEdges, + const vector& maEdgesTwin, + size_t & i) const +{ + // if there is a concave vertex between EDGEs + // then position of a dividing BranchPoint is undefined, it is somewhere + // on an arc-shaped part of the Branch around the concave vertex. + // Chose this position by a VERTEX of the opposite EDGE, or put it in the middle + // of the arc if there is no opposite VERTEX. + // All null-length segments around a VERTEX belong to one of EDGEs. + + BranchPoint divisionPnt; + divisionPnt._branch = this; + + size_t ie1 = getGeomEdge( maEdges [i] ); + size_t ie2 = getGeomEdge( maEdgesTwin[i] ); + + size_t iSeg1 = getBndSegment( maEdges[ i-1 ] ); + size_t iSeg2 = getBndSegment( maEdges[ i ] ); + bool isConcaPrev = _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 ); + bool isConcaNext = _boundary->IsConcaveSegment( ie1, iSeg2 ); + if ( !isConcaNext && !isConcaPrev ) + return false; + + bool isConcaveV = false; + + int iPrev = i-1, iNext = i; + if ( isConcaNext ) // all null-length segments follow + { + // look for a VERTEX of the opposite EDGE + ++iNext; // end of null-length segments + while ( iNext < maEdges.size() ) + { + iSeg2 = getBndSegment( maEdges[ iNext ] ); + if ( _boundary->IsConcaveSegment( ie1, iSeg2 )) + ++iNext; + else + break; + } + bool vertexFound = false; + for ( size_t iE = i+1; iE < iNext; ++iE ) + { + ie2 = getGeomEdge( maEdgesTwin[iE] ); + if ( ie2 != edgeIDs2.back() ) + { + // opposite VERTEX found + divisionPnt._iEdge = iE; + divisionPnt._edgeParam = 0; + divPoints.push_back( divisionPnt ); + edgeIDs1.push_back( ie1 ); + edgeIDs2.push_back( ie2 ); + vertexFound = true; + } + } + if ( vertexFound ) + { + i = --iNext; + isConcaveV = true; + } + } + else if ( isConcaPrev ) + { + // all null-length segments passed, find their beginning + while ( iPrev-1 >= 0 ) + { + iSeg1 = getBndSegment( maEdges[ iPrev-1 ] ); + if ( _boundary->IsConcaveSegment( edgeIDs1.back(), iSeg1 )) + --iPrev; + else + break; + } + } + + if ( iPrev < i-1 || iNext > i ) + { + // no VERTEX on the opposite EDGE, put the Branch Point in the middle + double par1 = _params[ iPrev ], par2 = _params[ iNext ]; + double midPar = 0.5 * ( par1 + par2 ); + divisionPnt._iEdge = iPrev; + while ( _params[ divisionPnt._iEdge + 1 ] < midPar ) + ++divisionPnt._iEdge; + divisionPnt._edgeParam = + ( _params[ divisionPnt._iEdge + 1 ] - midPar ) / + ( _params[ divisionPnt._iEdge + 1 ] - _params[ divisionPnt._iEdge ] ); + divPoints.push_back( divisionPnt ); + isConcaveV = true; + } + + return isConcaveV; +} + +//================================================================================ +/*! + * \brief Return indices of opposite parts of EDGEs equidistant from this branch + * \param [out] edgeIDs1 - EDGE index opposite to the edgeIDs2[i]-th EDGE + * \param [out] edgeIDs2 - EDGE index opposite to the edgeIDs1[i]-th EDGE + * \param [out] divPoints - BranchPoint's located between two successive unique + * pairs of EDGE indices. A \a divPoints[i] can separate e.g. two following pairs + * of EDGE indices < 0, 2 > and < 0, 1 >. Number of \a divPoints is one less + * than number of \a edgeIDs + */ +//================================================================================ + +void SMESH_MAT2d::Branch::getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2, + std::vector< BranchPoint >& divPoints) const +{ + edgeIDs1.clear(); + edgeIDs2.clear(); + divPoints.clear(); + + edgeIDs1.push_back( getGeomEdge( _maEdges[0] )); + edgeIDs2.push_back( getGeomEdge( _maEdges[0]->twin() )); + + std::vector twins( _maEdges.size() ); + for ( size_t i = 0; i < _maEdges.size(); ++i ) + twins[i] = _maEdges[i]->twin(); + + // size_t lastConcaE1 = _boundary.nbEdges(); + // size_t lastConcaE2 = _boundary.nbEdges(); + + BranchPoint divisionPnt; + divisionPnt._branch = this; + + for ( size_t i = 0; i < _maEdges.size(); ++i ) + { + size_t ie1 = getGeomEdge( _maEdges[i] ); + size_t ie2 = getGeomEdge( _maEdges[i]->twin() ); + + if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ) + { + bool isConcaveV = false; + if ( edgeIDs1.back() != ie1 && edgeIDs2.back() == ie2 ) + { + isConcaveV = addDivPntForConcaVertex( edgeIDs1, edgeIDs2, divPoints, _maEdges, twins, i ); + } + if ( edgeIDs1.back() == ie1 && edgeIDs2.back() != ie2 ) + { + isConcaveV = addDivPntForConcaVertex( edgeIDs2, edgeIDs1, divPoints, twins, _maEdges, i ); + } + + if ( isConcaveV ) + { + ie1 = getGeomEdge( _maEdges[i] ); + ie2 = getGeomEdge( _maEdges[i]->twin() ); + } + if (( !isConcaveV ) || + ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 )) + { + edgeIDs1.push_back( ie1 ); + edgeIDs2.push_back( ie2 ); + } + if ( divPoints.size() < edgeIDs1.size() - 1 ) + { + divisionPnt._iEdge = i; + divisionPnt._edgeParam = 0; + divPoints.push_back( divisionPnt ); + } + + } // if ( edgeIDs1.back() != ie1 || edgeIDs2.back() != ie2 ) + } // loop on _maEdges +} + +//================================================================================ +/*! + * \brief Store data of boundary segments in TVDEdge + */ +//================================================================================ + +void SMESH_MAT2d::Branch::setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge ) +{ + if ( maEdge ) maEdge->cell()->color( geomIndex ); +} +std::size_t SMESH_MAT2d::Branch::getGeomEdge( const TVDEdge* maEdge ) +{ + return maEdge ? maEdge->cell()->color() : std::string::npos; +} +void SMESH_MAT2d::Branch::setBndSegment( std::size_t segIndex, const TVDEdge* maEdge ) +{ + if ( maEdge ) maEdge->color( segIndex ); +} +std::size_t SMESH_MAT2d::Branch::getBndSegment( const TVDEdge* maEdge ) +{ + return maEdge ? maEdge->color() : std::string::npos; +} + +//================================================================================ +/*! + * \brief Returns a boundary point on a given EDGE + * \param [in] iEdge - index of the EDGE within MedialAxis + * \param [in] iSeg - index of a boundary segment within this Branch + * \param [in] u - [0;1] normalized param within \a iSeg-th segment + * \param [out] bp - the found BoundaryPoint + * \return bool - true if the BoundaryPoint is found + */ +//================================================================================ + +bool SMESH_MAT2d::Boundary::getPoint( std::size_t iEdge, + std::size_t iSeg, + double u, + BoundaryPoint& bp ) const +{ + if ( iEdge >= _pointsPerEdge.size() ) + return false; + if ( iSeg+1 >= _pointsPerEdge[ iEdge ]._params.size() ) + return false; + + // This method is called by Branch that can have an opposite orientation, + // hence u is inverted depending on orientation coded as a sign of _maEdge index + bool isReverse = ( _pointsPerEdge[ iEdge ]._maEdges[ iSeg ].second < 0 ); + if ( isReverse ) + u = 1. - u; + + double p0 = _pointsPerEdge[ iEdge ]._params[ iSeg ]; + double p1 = _pointsPerEdge[ iEdge ]._params[ iSeg+1 ]; + + bp._param = p0 * ( 1. - u ) + p1 * u; + bp._edgeIndex = iEdge; + + return true; +} + diff --git a/src/SMESHUtils/SMESH_MAT2d.hxx b/src/SMESHUtils/SMESH_MAT2d.hxx new file mode 100644 index 000000000..7e1061d05 --- /dev/null +++ b/src/SMESHUtils/SMESH_MAT2d.hxx @@ -0,0 +1,224 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : SMESH_MAT2d.hxx +// Created : Thu May 28 17:49:53 2015 +// Author : Edward AGAPOV (eap) + +#ifndef __SMESH_MAT2d_HXX__ +#define __SMESH_MAT2d_HXX__ + +#include "SMESH_Utils.hxx" + +#include +#include + +#include +#include + +#include +#include + +class Adaptor3d_Curve; + +// Medial Axis Transform 2D +namespace SMESH_MAT2d +{ + class MedialAxis; // MedialAxis is the entry point + class Branch; + class BranchEnd; + class Boundary; + struct BoundaryPoint; + + typedef boost::polygon::voronoi_diagram TVD; + typedef TVD::cell_type TVDCell; + typedef TVD::edge_type TVDEdge; + typedef TVD::vertex_type TVDVertex; + + //------------------------------------------------------------------------------------- + // type of Branch end point + enum BranchEndType { BE_UNDEF, + BE_ON_VERTEX, // branch ends at a convex VRTEX + BE_BRANCH_POINT, // branch meats 2 or more other branches + BE_END // branch end equidistant from several adjacent segments + }; + //------------------------------------------------------------------------------------- + /*! + * \brief End point of MA Branch + */ + struct SMESHUtils_EXPORT BranchEnd + { + const TVDVertex* _vertex; + BranchEndType _type; + std::vector< const Branch* > _branches; + + BranchEnd(): _vertex(0), _type( BE_UNDEF ) {} + }; + //------------------------------------------------------------------------------------- + /*! + * \brief Point on MA Branch + */ + struct SMESHUtils_EXPORT BranchPoint + { + const Branch* _branch; + std::size_t _iEdge; // MA edge index within the branch + double _edgeParam; // normalized param within the MA edge + + BranchPoint(): _branch(0), _iEdge(0), _edgeParam(-1) {} + }; + //------------------------------------------------------------------------------------- + /*! + * \brief Branch is a set of MA edges enclosed between branch points and/or MA ends. + * It's main feature is to return two BoundaryPoint's per a point on it. + */ + class SMESHUtils_EXPORT Branch + { + public: + bool getBoundaryPoints(double param, BoundaryPoint& bp1, BoundaryPoint& bp2 ) const; + bool getBoundaryPoints(std::size_t iMAEdge, double maEdgeParam, + BoundaryPoint& bp1, BoundaryPoint& bp2 ) const; + bool getBoundaryPoints(const BranchPoint& p, + BoundaryPoint& bp1, BoundaryPoint& bp2 ) const; + bool getParameter(const BranchPoint& p, double & u ) const; + + std::size_t nbEdges() const { return _maEdges.size(); } + + const BranchEnd* getEnd(bool the2nd) const { return & ( the2nd ? _endPoint2 : _endPoint1 ); } + + bool hasEndOfType(BranchEndType type) const; + + void getPoints( std::vector< gp_XY >& points, const double scale[2]) const; + + void getGeomEdges( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2 ) const; + + void getOppositeGeomEdges( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2, + std::vector< BranchPoint >& divPoints) const; + + // construction + void init( std::vector& maEdges, + const Boundary* boundary, + std::map< const TVDVertex*, BranchEndType > endType); + void setBranchesToEnds( const std::vector< Branch >& branches); + + static void setGeomEdge( std::size_t geomIndex, const TVDEdge* maEdge ); + static std::size_t getGeomEdge( const TVDEdge* maEdge ); + static void setBndSegment( std::size_t segIndex, const TVDEdge* maEdge ); + static std::size_t getBndSegment( const TVDEdge* maEdge ); + + private: + + bool addDivPntForConcaVertex( std::vector< std::size_t >& edgeIDs1, + std::vector< std::size_t >& edgeIDs2, + std::vector< BranchPoint >& divPoints, + const std::vector& maEdges, + const std::vector& maEdgesTwin, + size_t & i) const; + + // association of _maEdges with boundary segments is stored in this way: + // index of an EDGE: TVDEdge->cell()->color() + // index of a segment on EDGE: TVDEdge->color() + std::vector _maEdges; // MA edges ending at points located at _params + std::vector _params; // params of points on MA, normalized [0;1] within this branch + const Boundary* _boundary; // face boundary + BranchEnd _endPoint1; + BranchEnd _endPoint2; + }; + + //------------------------------------------------------------------------------------- + /*! + * \brief Data of a discretized EDGE allowing to get a point on MA by a parameter on EDGE + */ + struct BndPoints + { + std::vector< double > _params; // params of discretization points on an EDGE + std::vector< std::pair< const Branch*, int > > _maEdges; /* index of TVDEdge in branch; + index sign means orientation; + index == Branch->nbEdges() means + end point of a Branch */ + }; + //------------------------------------------------------------------------------------- + /*! + * \brief Face boundary is discretized so that each its segment to correspond to + * an edge of MA + */ + class Boundary + { + public: + + Boundary( std::size_t nbEdges ): _pointsPerEdge( nbEdges ) {} + BndPoints& getPoints( std::size_t iEdge ) { return _pointsPerEdge[ iEdge ]; } + std::size_t nbEdges() const { return _pointsPerEdge.size(); } + + bool getPoint( std::size_t iEdge, std::size_t iSeg, double u, BoundaryPoint& bp ) const; + + bool getBranchPoint( const std::size_t iEdge, double u, BranchPoint& p ) const; + + bool IsConcaveSegment( std::size_t iEdge, std::size_t iSeg ) const; + + private: + std::vector< BndPoints > _pointsPerEdge; + }; + + //------------------------------------------------------------------------------------- + /*! + * \brief Point on FACE boundary + */ + struct SMESHUtils_EXPORT BoundaryPoint + { + std::size_t _edgeIndex; // index of an EDGE in a sequence passed to MedialAxis() + double _param; // parameter of this EDGE + }; + //------------------------------------------------------------------------------------- + /*! + * \brief Medial axis (MA) is defined as the loci of centres of locally + * maximal balls inside 2D representation of a face. This class + * implements a piecewise approximation of MA. + */ + class SMESHUtils_EXPORT MedialAxis + { + public: + MedialAxis(const TopoDS_Face& face, + const std::vector< TopoDS_Edge >& edges, + const double minSegLen, + const bool ignoreCorners = false ); + const Boundary& getBoundary() const { return _boundary; } + const std::vector< Branch >& getBranches() const { return _branch; } + const std::vector< const BranchEnd* >& getBranchPoints() const { return _branchPnt; } + + void getPoints( const Branch& branch, std::vector< gp_XY >& points) const; + Adaptor3d_Curve* make3DCurve(const Branch& branch) const; + + private: + + private: + TopoDS_Face _face; + TVD _vd; + std::vector< Branch > _branch; + std::vector< const BranchEnd* > _branchPnt; + Boundary _boundary; + double _scale[2]; + }; + +} + +#endif diff --git a/src/SMESH_SWIG/StdMeshersBuilder.py b/src/SMESH_SWIG/StdMeshersBuilder.py index 41e0c2269..3a07810ad 100644 --- a/src/SMESH_SWIG/StdMeshersBuilder.py +++ b/src/SMESH_SWIG/StdMeshersBuilder.py @@ -42,6 +42,8 @@ Hexa = "Hexa_3D" QUADRANGLE = "Quadrangle_2D" ## Algorithm type: Radial Quadrangle 1D-2D algorithm, see StdMeshersBuilder_RadialQuadrangle1D2D RADIAL_QUAD = "RadialQuadrangle_1D2D" +## Algorithm type: Quadrangle (Medial Axis Projection) 1D-2D algorithm, see StdMeshersBuilder_QuadMA_1D2D +QUAD_MA_PROJ = "QuadFromMedialAxis_1D2D" # import items of enums for e in StdMeshers.QuadType._items: exec('%s = StdMeshers.%s'%(e,e)) @@ -455,9 +457,6 @@ class StdMeshersBuilder_Segment_Python(Mesh_Algorithm): algoType = PYTHON ## doc string of the method # @internal - docHelper = "Creates tetrahedron 3D algorithm for solids" - ## doc string of the method - # @internal docHelper = "Creates segment 1D algorithm for edges" ## Private constructor. @@ -856,7 +855,7 @@ class StdMeshersBuilder_Projection1D2D(StdMeshersBuilder_Projection2D): algoType = "Projection_1D2D" ## doc string of the method # @internal - docHelper = "Creates projection 1D-2D algorithm for edges and faces" + docHelper = "Creates projection 1D-2D algorithm for faces" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1115,7 +1114,7 @@ class StdMeshersBuilder_RadialPrism3D(StdMeshersBuilder_Prism3D): algoType = "RadialPrism_3D" ## doc string of the method # @internal - docHelper = "Creates prism 3D algorithm for volumes" + docHelper = "Creates Raial Prism 3D algorithm for volumes" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1147,7 +1146,7 @@ class StdMeshersBuilder_RadialQuadrangle1D2D(Mesh_Algorithm): algoType = RADIAL_QUAD ## doc string of the method # @internal - docHelper = "Creates quadrangle 1D-2D algorithm for triangular faces" + docHelper = "Creates quadrangle 1D-2D algorithm for faces having a shape of disk or a disk segment" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1258,6 +1257,35 @@ class StdMeshersBuilder_RadialQuadrangle1D2D(Mesh_Algorithm): pass # end of StdMeshersBuilder_RadialQuadrangle1D2D class +## Defines a Quadrangle (Medial Axis Projection) 1D-2D algorithm +# +# It is created by calling smeshBuilder.Mesh.Quadrangle(smeshBuilder.QUAD_MA_PROJ,geom=0) +# +# @ingroup l2_algos_quad_ma +class StdMeshersBuilder_QuadMA_1D2D(Mesh_Algorithm): + + ## name of the dynamic method in smeshBuilder.Mesh class + # @internal + meshMethod = "Quadrangle" + ## type of algorithm used with helper function in smeshBuilder.Mesh class + # @internal + algoType = QUAD_MA_PROJ + ## doc string of the method + # @internal + docHelper = "Creates quadrangle 1D-2D algorithm for faces" + + ## Private constructor. + # @param mesh parent mesh object algorithm is assigned to + # @param geom geometry (shape/sub-shape) algorithm is assigned to; + # if it is @c 0 (default), the algorithm is assigned to the main shape + def __init__(self, mesh, geom=0): + Mesh_Algorithm.__init__(self) + self.Create(mesh, geom, self.algoType) + pass + + pass + + ## Defines a Use Existing Elements 1D algorithm # # It is created by calling smeshBuilder.Mesh.UseExisting1DElements(geom=0) @@ -1327,7 +1355,7 @@ class StdMeshersBuilder_UseExistingElements_1D2D(Mesh_Algorithm): isDefault = True ## doc string of the method # @internal - docHelper = "Creates 1D-2D algorithm for edges/faces with reusing of existing mesh elements" + docHelper = "Creates 1D-2D algorithm for faces with reusing of existing mesh elements" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1375,7 +1403,7 @@ class StdMeshersBuilder_Cartesian_3D(Mesh_Algorithm): isDefault = True ## doc string of the method # @internal - docHelper = "Creates body fitting 3D algorithm for volumes" + docHelper = "Creates Body Fitting 3D algorithm for volumes" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1505,7 +1533,7 @@ class StdMeshersBuilder_UseExisting_1D(Mesh_Algorithm): algoType = "UseExisting_1D" ## doc string of the method # @internal - docHelper = "Creates 1D algorithm for edges with reusing of existing mesh elements" + docHelper = "Creates 1D algorithm allowing batch meshing of edges" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to @@ -1533,7 +1561,7 @@ class StdMeshersBuilder_UseExisting_2D(Mesh_Algorithm): algoType = "UseExisting_2D" ## doc string of the method # @internal - docHelper = "Creates 2D algorithm for faces with reusing of existing mesh elements" + docHelper = "Creates 2D algorithm allowing batch meshing of faces" ## Private constructor. # @param mesh parent mesh object algorithm is assigned to diff --git a/src/StdMeshers/CMakeLists.txt b/src/StdMeshers/CMakeLists.txt index 46b756b75..8714c1b4f 100644 --- a/src/StdMeshers/CMakeLists.txt +++ b/src/StdMeshers/CMakeLists.txt @@ -130,6 +130,7 @@ SET(StdMeshers_HEADERS StdMeshers_Projection_1D2D.hxx StdMeshers_CartesianParameters3D.hxx StdMeshers_Cartesian_3D.hxx + StdMeshers_QuadFromMedialAxis_1D2D.hxx ) IF(SALOME_SMESH_ENABLE_MEFISTO) @@ -193,6 +194,7 @@ SET(StdMeshers_SOURCES StdMeshers_CartesianParameters3D.cxx StdMeshers_Cartesian_3D.cxx StdMeshers_Adaptive1D.cxx + StdMeshers_QuadFromMedialAxis_1D2D.cxx ) IF(SALOME_SMESH_ENABLE_MEFISTO) diff --git a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx index 0616b630d..264189170 100644 --- a/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx +++ b/src/StdMeshers/StdMeshers_MEFISTO_2D.cxx @@ -696,7 +696,7 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, double xmax = -1.e300; double ymin = 1.e300; double ymax = -1.e300; - int nbp = 23; + const int nbp = 23; scalex = 1; scaley = 1; @@ -720,13 +720,8 @@ void StdMeshers_MEFISTO_2D::ComputeScaleOnFace(SMESH_Mesh & aMesh, ymin = p.Y(); if (p.Y() > ymax) ymax = p.Y(); - // MESSAGE(" "<< f<<" "< maxratio) { - SCRUTE( scaley ); scaley *= xyratio / maxratio; - SCRUTE( scaley ); } else if (xyratio < 1./maxratio) { - SCRUTE( scalex ); scalex *= 1 / xyratio / maxratio; - SCRUTE( scalex ); } - ASSERT(scalex); - ASSERT(scaley); } // namespace diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx new file mode 100644 index 000000000..a6d372d5f --- /dev/null +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.cxx @@ -0,0 +1,1251 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : StdMeshers_QuadFromMedialAxis_1D2D.cxx +// Created : Wed Jun 3 17:33:45 2015 +// Author : Edward AGAPOV (eap) + +#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx" + +#include "SMESH_Block.hxx" +#include "SMESH_Gen.hxx" +#include "SMESH_MAT2d.hxx" +#include "SMESH_Mesh.hxx" +#include "SMESH_MesherHelper.hxx" +#include "SMESH_ProxyMesh.hxx" +#include "SMESH_subMesh.hxx" +#include "StdMeshers_FaceSide.hxx" +#include "StdMeshers_Regular_1D.hxx" +#include "StdMeshers_ViscousLayers2D.hxx" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +//================================================================================ +/*! + * \brief 1D algo + */ +class StdMeshers_QuadFromMedialAxis_1D2D::Algo1D : public StdMeshers_Regular_1D +{ +public: + Algo1D(int studyId, SMESH_Gen* gen): + StdMeshers_Regular_1D( gen->GetANewId(), studyId, gen ) + { + } + void SetSegmentLength( double len ) + { + _value[ BEG_LENGTH_IND ] = len; + _value[ PRECISION_IND ] = 1e-7; + _hypType = LOCAL_LENGTH; + } +}; + +//================================================================================ +/*! + * \brief Constructor sets algo features + */ +//================================================================================ + +StdMeshers_QuadFromMedialAxis_1D2D::StdMeshers_QuadFromMedialAxis_1D2D(int hypId, + int studyId, + SMESH_Gen* gen) + : StdMeshers_Quadrangle_2D(hypId, studyId, gen), + _regular1D( 0 ) +{ + _name = "QuadFromMedialAxis_1D2D"; + _shapeType = (1 << TopAbs_FACE); + _onlyUnaryInput = true; // FACE by FACE so far + _requireDiscreteBoundary = false; // make 1D by myself + _supportSubmeshes = true; // make 1D by myself + _neededLowerHyps[ 1 ] = true; // suppress warning on hiding a global 1D algo + _neededLowerHyps[ 2 ] = true; // suppress warning on hiding a global 2D algo + _compatibleHypothesis.clear(); + //_compatibleHypothesis.push_back("ViscousLayers2D"); +} + +//================================================================================ +/*! + * \brief Destructor + */ +//================================================================================ + +StdMeshers_QuadFromMedialAxis_1D2D::~StdMeshers_QuadFromMedialAxis_1D2D() +{ + delete _regular1D; + _regular1D = 0; +} + +//================================================================================ +/*! + * \brief Check if needed hypotheses are present + */ +//================================================================================ + +bool StdMeshers_QuadFromMedialAxis_1D2D::CheckHypothesis(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus) +{ + return true; // does not require hypothesis +} + +namespace +{ + //================================================================================ + /*! + * \brief Temporary mesh + */ + struct TmpMesh : public SMESH_Mesh + { + TmpMesh() + { + _myMeshDS = new SMESHDS_Mesh(/*id=*/0, /*isEmbeddedMode=*/true); + } + }; + + //================================================================================ + /*! + * \brief Select two EDGEs from a map, either mapped to least values or to max values + */ + //================================================================================ + + // template< class TVal2EdgesMap > + // void getTwo( bool least, + // TVal2EdgesMap& map, + // vector& twoEdges, + // vector& otherEdges) + // { + // twoEdges.clear(); + // otherEdges.clear(); + // if ( least ) + // { + // TVal2EdgesMap::iterator i = map.begin(); + // twoEdges.push_back( i->second ); + // twoEdges.push_back( ++i->second ); + // for ( ; i != map.end(); ++i ) + // otherEdges.push_back( i->second ); + // } + // else + // { + // TVal2EdgesMap::reverse_iterator i = map.rbegin(); + // twoEdges.push_back( i->second ); + // twoEdges.push_back( ++i->second ); + // for ( ; i != map.rend(); ++i ) + // otherEdges.push_back( i->second ); + // } + // TopoDS_Vertex v; + // if ( TopExp::CommonVertex( twoEdges[0], twoEdges[1], v )) + // { + // twoEdges.clear(); // two EDGEs must not be connected + // otherEdges.clear(); + // } + // } + + //================================================================================ + /*! + * \brief Finds out a minimal segment length given EDGEs will be divided into. + * This length is further used to discretize the Medial Axis + */ + //================================================================================ + + double getMinSegLen(SMESH_MesherHelper& theHelper, + const vector& theEdges) + { + TmpMesh tmpMesh; + SMESH_Mesh* mesh = theHelper.GetMesh(); + + vector< SMESH_Algo* > algos( theEdges.size() ); + for ( size_t i = 0; i < theEdges.size(); ++i ) + { + SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] ); + algos[i] = sm->GetAlgo(); + } + + const int nbSegDflt = mesh->GetGen()->GetDefaultNbSegments(); + double minSegLen = Precision::Infinite(); + + for ( size_t i = 0; i < theEdges.size(); ++i ) + { + SMESH_subMesh* sm = mesh->GetSubMesh( theEdges[i] ); + if ( SMESH_Algo::IsStraight( theEdges[i], /*degenResult=*/true )) + continue; + // get algo + size_t iOpp = ( theEdges.size() == 4 ? (i+2)%4 : i ); + SMESH_Algo* algo = sm->GetAlgo(); + if ( !algo ) algo = algos[ iOpp ]; + // get hypo + SMESH_Hypothesis::Hypothesis_Status status = SMESH_Hypothesis::HYP_MISSING; + if ( algo ) + { + if ( !algo->CheckHypothesis( *mesh, theEdges[i], status )) + algo->CheckHypothesis( *mesh, theEdges[iOpp], status ); + } + // compute + if ( status != SMESH_Hypothesis::HYP_OK ) + { + minSegLen = Min( minSegLen, SMESH_Algo::EdgeLength( theEdges[i] ) / nbSegDflt ); + } + else + { + tmpMesh.Clear(); + tmpMesh.ShapeToMesh( TopoDS_Shape()); + tmpMesh.ShapeToMesh( theEdges[i] ); + try { + mesh->GetGen()->Compute( tmpMesh, theEdges[i], true, true ); // make nodes on VERTEXes + if ( !algo->Compute( tmpMesh, theEdges[i] )) + continue; + } + catch (...) { + continue; + } + SMDS_EdgeIteratorPtr segIt = tmpMesh.GetMeshDS()->edgesIterator(); + while ( segIt->more() ) + { + const SMDS_MeshElement* seg = segIt->next(); + double len = SMESH_TNodeXYZ( seg->GetNode(0) ).Distance( seg->GetNode(1) ); + minSegLen = Min( minSegLen, len ); + } + } + } + if ( Precision::IsInfinite( minSegLen )) + minSegLen = mesh->GetShapeDiagonalSize() / nbSegDflt; + + return minSegLen; + } + + //================================================================================ + /*! + * \brief Returns EDGEs located between two VERTEXes at which given MA branches end + * \param [in] br1 - one MA branch + * \param [in] br2 - one more MA branch + * \param [in] allEdges - all EDGEs of a FACE + * \param [out] shortEdges - the found EDGEs + * \return bool - is OK or not + */ + //================================================================================ + + bool getConnectedEdges( const SMESH_MAT2d::Branch* br1, + const SMESH_MAT2d::Branch* br2, + const vector& allEdges, + vector& shortEdges) + { + vector< size_t > edgeIDs[4]; + br1->getGeomEdges( edgeIDs[0], edgeIDs[1] ); + br2->getGeomEdges( edgeIDs[2], edgeIDs[3] ); + + // EDGEs returned by a Branch form a connected chain with a VERTEX where + // the Branch ends at the chain middle. One of end EDGEs of the chain is common + // with either end EDGE of the chain of the other Branch, or the chains are connected + // at a common VERTEX; + + // Get indices of end EDGEs of the branches + bool vAtStart1 = ( br1->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX ); + bool vAtStart2 = ( br2->getEnd(0)->_type == SMESH_MAT2d::BE_ON_VERTEX ); + size_t iEnd[4] = { + vAtStart1 ? edgeIDs[0].back() : edgeIDs[0][0], + vAtStart1 ? edgeIDs[1].back() : edgeIDs[1][0], + vAtStart2 ? edgeIDs[2].back() : edgeIDs[2][0], + vAtStart2 ? edgeIDs[3].back() : edgeIDs[3][0] + }; + + set< size_t > connectedIDs; + TopoDS_Vertex vCommon; + // look for the same EDGEs + for ( int i = 0; i < 2; ++i ) + for ( int j = 2; j < 4; ++j ) + if ( iEnd[i] == iEnd[j] ) + { + connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() ); + connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() ); + i = j = 4; + } + if ( connectedIDs.empty() ) + // look for connected EDGEs + for ( int i = 0; i < 2; ++i ) + for ( int j = 2; j < 4; ++j ) + if ( TopExp::CommonVertex( allEdges[ iEnd[i]], allEdges[ iEnd[j]], vCommon )) + { + connectedIDs.insert( edgeIDs[i].begin(), edgeIDs[i].end() ); + connectedIDs.insert( edgeIDs[j].begin(), edgeIDs[j].end() ); + i = j = 4; + } + if ( connectedIDs.empty() || // nothing + allEdges.size() - connectedIDs.size() < 2 ) // too many + return false; + + // set shortEdges in the order as in allEdges + if ( connectedIDs.count( 0 ) && + connectedIDs.count( allEdges.size()-1 )) + { + size_t iE = allEdges.size()-1; + while ( connectedIDs.count( iE-1 )) + --iE; + for ( size_t i = 0; i < connectedIDs.size(); ++i ) + { + shortEdges.push_back( allEdges[ iE ]); + iE = ( iE + 1 ) % allEdges.size(); + } + } + else + { + set< size_t >::iterator i = connectedIDs.begin(); + for ( ; i != connectedIDs.end(); ++i ) + shortEdges.push_back( allEdges[ *i ]); + } + return true; + } + + //================================================================================ + /*! + * \brief Find EDGEs to discretize using projection from MA + * \param [in] theFace - the FACE to be meshed + * \param [in] theWire - ordered EDGEs of the FACE + * \param [out] theSinuEdges - the EDGEs to discretize using projection from MA + * \param [out] theShortEdges - other EDGEs + * \return bool - OK or not + * + * Is separate all EDGEs into four sides of a quadrangle connected in the order: + * theSinuEdges[0], theShortEdges[0], theSinuEdges[1], theShortEdges[1] + */ + //================================================================================ + + bool getSinuousEdges( SMESH_MesherHelper& theHelper, + const TopoDS_Face& theFace, + list& theWire, + vector theSinuEdges[2], + vector theShortEdges[2]) + { + theSinuEdges[0].clear(); + theSinuEdges[1].clear(); + theShortEdges[0].clear(); + theShortEdges[1].clear(); + + vector allEdges( theWire.begin(), theWire.end() ); + const size_t nbEdges = allEdges.size(); + if ( nbEdges < 4 ) + return false; + + // create MedialAxis to find short edges by analyzing MA branches + double minSegLen = getMinSegLen( theHelper, allEdges ); + SMESH_MAT2d::MedialAxis ma( theFace, allEdges, minSegLen ); + + // in an initial request case, theFace represents a part of a river with almost parallel banks + // so there should be two branch points + using SMESH_MAT2d::BranchEnd; + using SMESH_MAT2d::Branch; + const vector< const BranchEnd* >& braPoints = ma.getBranchPoints(); + if ( braPoints.size() < 2 ) + return false; + TopTools_MapOfShape shortMap; + size_t nbBranchPoints = 0; + for ( size_t i = 0; i < braPoints.size(); ++i ) + { + vector< const Branch* > vertBranches; // branches with an end on VERTEX + for ( size_t ib = 0; ib < braPoints[i]->_branches.size(); ++ib ) + { + const Branch* branch = braPoints[i]->_branches[ ib ]; + if ( branch->hasEndOfType( SMESH_MAT2d::BE_ON_VERTEX )) + vertBranches.push_back( branch ); + } + if ( vertBranches.size() != 2 || braPoints[i]->_branches.size() != 3) + continue; + + // get common EDGEs of two branches + if ( !getConnectedEdges( vertBranches[0], vertBranches[1], + allEdges, theShortEdges[ nbBranchPoints > 0 ] )) + return false; + + for ( size_t iS = 0; iS < theShortEdges[ nbBranchPoints ].size(); ++iS ) + shortMap.Add( theShortEdges[ nbBranchPoints ][ iS ]); + + ++nbBranchPoints; + } + + if ( nbBranchPoints != 2 ) + return false; + + // add to theSinuEdges all edges that are not theShortEdges + vector< vector > sinuEdges(1); + TopoDS_Vertex vCommon; + for ( size_t i = 0; i < allEdges.size(); ++i ) + { + if ( !shortMap.Contains( allEdges[i] )) + { + if ( !sinuEdges.back().empty() ) + if ( !TopExp::CommonVertex( sinuEdges.back().back(), allEdges[ i ], vCommon )) + sinuEdges.resize( sinuEdges.size() + 1 ); + + sinuEdges.back().push_back( allEdges[i] ); + } + } + if ( sinuEdges.size() == 3 ) + { + if ( !TopExp::CommonVertex( sinuEdges.back().back(), sinuEdges[0][0], vCommon )) + return false; + vector& last = sinuEdges.back(); + last.insert( last.end(), sinuEdges[0].begin(), sinuEdges[0].end() ); + sinuEdges[0].swap( last ); + sinuEdges.resize( 2 ); + } + if ( sinuEdges.size() != 2 ) + return false; + + theSinuEdges[0].swap( sinuEdges[0] ); + theSinuEdges[1].swap( sinuEdges[1] ); + + if ( !TopExp::CommonVertex( theSinuEdges[0].back(), theShortEdges[0][0], vCommon ) || + !vCommon.IsSame( theHelper.IthVertex( 1, theSinuEdges[0].back() ))) + theShortEdges[0].swap( theShortEdges[1] ); + + return ( theShortEdges[0].size() > 0 && theShortEdges[1].size() > 0 && + theSinuEdges [0].size() > 0 && theSinuEdges [1].size() > 0 ); + + // the sinuous EDGEs can be composite and C0 continuous, + // therefor we use a complex criterion to find TWO short non-sinuous EDGEs + // and the rest EDGEs will be treated as sinuous. + // A short edge should have the following features: + // a) straight + // b) short + // c) with convex corners at ends + // d) far from the other short EDGE + + // vector< double > isStraightEdge( nbEdges, 0 ); // criterion value + + // // a0) evaluate continuity + // const double contiWgt = 0.5; // weight of continuity in the criterion + // multimap< int, TopoDS_Edge > continuity; + // for ( size_t i = 0; i < nbEdges; ++I ) + // { + // BRepAdaptor_Curve curve( allEdges[i] ); + // GeomAbs_Shape C = GeomAbs_CN; + // try: + // C = curve.Continuity(); // C0, G1, C1, G2, C2, C3, CN + // catch ( Standard_Failure ) {} + // continuity.insert( make_pair( C, allEdges[i] )); + // isStraight[i] += double( C ) / double( CN ) * contiWgt; + // } + + // // try to choose by continuity + // int mostStraight = (int) continuity.rbegin()->first; + // int lessStraight = (int) continuity.begin()->first; + // if ( mostStraight != lessStraight ) + // { + // int nbStraight = continuity.count( mostStraight ); + // if ( nbStraight == 2 ) + // { + // getTwo( /*least=*/false, continuity, theShortEdges, theSinuEdges ); + // } + // else if ( nbStraight == 3 && nbEdges == 4 ) + // { + // theSinuEdges.push_back( continuity.begin()->second ); + // vector::iterator it = + // std::find( allEdges.begin(), allEdges.end(), theSinuEdges[0] ); + // int i = std::distance( allEdges.begin(), it ); + // theSinuEdges .push_back( allEdges[( i+2 )%4 ]); + // theShortEdges.push_back( allEdges[( i+1 )%4 ]); + // theShortEdges.push_back( allEdges[( i+3 )%4 ]); + // } + // if ( theShortEdges.size() == 2 ) + // return true; + // } + + // // a) curvature; evaluate aspect ratio + // { + // const double curvWgt = 0.5; + // for ( size_t i = 0; i < nbEdges; ++I ) + // { + // BRepAdaptor_Curve curve( allEdges[i] ); + // double curvature = 1; + // if ( !curve.IsClosed() ) + // { + // const double f = curve.FirstParameter(), l = curve.LastParameter(); + // gp_Pnt pf = curve.Value( f ), pl = curve.Value( l ); + // gp_Lin line( pf, pl.XYZ() - pf.XYZ() ); + // double distMax = 0; + // for ( double u = f; u < l; u += (l-f)/30. ) + // distMax = Max( distMax, line.SquareDistance( curve.Value( u ))); + // curvature = Sqrt( distMax ) / ( pf.Distance( pl )); + // } + // isStraight[i] += curvWgt / ( curvature + 1e-20 ); + // } + // } + // // b) length + // { + // const double lenWgt = 0.5; + // for ( size_t i = 0; i < nbEdges; ++I ) + // { + // double length = SMESH_Algo::Length( allEdges[i] ); + // if ( length > 0 ) + // isStraight[i] += lenWgt / length; + // } + // } + // // c) with convex corners at ends + // { + // const double cornerWgt = 0.25; + // for ( size_t i = 0; i < nbEdges; ++I ) + // { + // double convex = 0; + // int iPrev = SMESH_MesherHelper::WrapIndex( int(i)-1, nbEdges ); + // int iNext = SMESH_MesherHelper::WrapIndex( int(i)+1, nbEdges ); + // TopoDS_Vertex v = helper.IthVertex( 0, allEdges[i] ); + // double angle = SMESH_MesherHelper::GetAngle( allEdges[iPrev], allEdges[i], theFace, v ); + // if ( angle < M_PI ) // [-PI; PI] + // convex += ( angle + M_PI ) / M_PI / M_PI; + // v = helper.IthVertex( 1, allEdges[i] ); + // angle = SMESH_MesherHelper::GetAngle( allEdges[iNext], allEdges[i], theFace, v ); + // if ( angle < M_PI ) // [-PI; PI] + // convex += ( angle + M_PI ) / M_PI / M_PI; + // isStraight[i] += cornerWgt * convex; + // } + // } + } + + //================================================================================ + /*! + * \brief Creates an EDGE from a sole branch of MA + */ + //================================================================================ + + TopoDS_Edge makeEdgeFromMA( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA ) + { + if ( theMA.getBranches().size() != 1 ) + return TopoDS_Edge(); + + vector< gp_XY > uv; + theMA.getPoints( theMA.getBranches()[0], uv ); + if ( uv.size() < 2 ) + return TopoDS_Edge(); + + TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); + Handle(Geom_Surface) surface = BRep_Tool::Surface( face ); + + // cout << "from salome.geom import geomBuilder" << endl; + // cout << "geompy = geomBuilder.New(salome.myStudy)" << endl; + Handle(TColgp_HArray1OfPnt) points = new TColgp_HArray1OfPnt(1, uv.size()); + for ( size_t i = 0; i < uv.size(); ++i ) + { + gp_Pnt p = surface->Value( uv[i].X(), uv[i].Y() ); + points->SetValue( i+1, p ); + //cout << "geompy.MakeVertex( "<< p.X()<<", " << p.Y()<<", " << p.Z()<<" )" << endl; + } + + GeomAPI_Interpolate interpol( points, /*isClosed=*/false, gp::Resolution()); + interpol.Perform(); + if ( !interpol.IsDone()) + return TopoDS_Edge(); + + TopoDS_Edge branchEdge = BRepBuilderAPI_MakeEdge(interpol.Curve()); + return branchEdge; + } + + //================================================================================ + /*! + * \brief Returns a type of shape, to which a hypothesis used to mesh a given edge is assigned + */ + //================================================================================ + + TopAbs_ShapeEnum getHypShape( SMESH_Mesh* mesh, const TopoDS_Shape& edge ) + { + TopAbs_ShapeEnum shapeType = TopAbs_SHAPE; + + SMESH_subMesh* sm = mesh->GetSubMesh( edge ); + SMESH_Algo* algo = sm->GetAlgo(); + if ( !algo ) return shapeType; + + const list & hyps = + algo->GetUsedHypothesis( *mesh, edge, /*ignoreAuxiliary=*/true ); + if ( hyps.empty() ) return shapeType; + + TopoDS_Shape shapeOfHyp = + SMESH_MesherHelper::GetShapeOfHypothesis( hyps.front(), edge, mesh); + + return SMESH_MesherHelper::GetGroupType( shapeOfHyp, /*woCompound=*/true); + } + + //================================================================================ + /*! + * \brief Discretize a sole branch of MA an returns parameters of divisions on MA + */ + //================================================================================ + + bool divideMA( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA, + const vector& theSinuEdges, + const size_t theSinuSide0Size, + SMESH_Algo* the1dAlgo, + vector& theMAParams ) + { + // check if all EDGEs of one size are meshed, then MA discretization is not needed + SMESH_Mesh* mesh = theHelper.GetMesh(); + size_t nbComputedEdges[2] = { 0, 0 }; + for ( size_t i = 1; i < theSinuEdges.size(); ++i ) + { + bool isComputed = ( ! mesh->GetSubMesh( theSinuEdges[i] )->IsEmpty() ); + nbComputedEdges[ i < theSinuSide0Size ] += isComputed; + } + if ( nbComputedEdges[0] == theSinuSide0Size || + nbComputedEdges[1] == theSinuEdges.size() - theSinuSide0Size ) + return true; // discretization is not needed + + + TopoDS_Edge branchEdge = makeEdgeFromMA( theHelper, theMA ); + if ( branchEdge.IsNull() ) + return false; + + // const char* file = "/misc/dn25/salome/eap/salome/misc/tmp/MAedge.brep"; + // BRepTools::Write( branchEdge, file); + // cout << "Write " << file << endl; + + // look for a most local hyps assigned to theSinuEdges + TopoDS_Edge edge = theSinuEdges[0]; + int mostSimpleShape = (int) getHypShape( mesh, edge ); + for ( size_t i = 1; i < theSinuEdges.size(); ++i ) + { + int shapeType = (int) getHypShape( mesh, theSinuEdges[i] ); + if ( shapeType > mostSimpleShape ) + edge = theSinuEdges[i]; + } + + SMESH_Algo* algo = the1dAlgo; + if ( mostSimpleShape != TopAbs_SHAPE ) + { + algo = mesh->GetSubMesh( edge )->GetAlgo(); + SMESH_Hypothesis::Hypothesis_Status status; + if ( !algo->CheckHypothesis( *mesh, edge, status )) + algo = the1dAlgo; + } + + TmpMesh tmpMesh; + tmpMesh.ShapeToMesh( branchEdge ); + try { + mesh->GetGen()->Compute( tmpMesh, branchEdge, true, true ); // make nodes on VERTEXes + if ( !algo->Compute( tmpMesh, branchEdge )) + return false; + } + catch (...) { + return false; + } + return SMESH_Algo::GetNodeParamOnEdge( tmpMesh.GetMeshDS(), branchEdge, theMAParams ); + } + + //================================================================================ + /*! + * \brief Modifies division parameters on MA to make them coincide with projections + * of VERTEXes to MA for a given pair of opposite EDGEs + * \param [in] theEdgePairInd - index of the EDGE pair + * \param [in] theDivPoints - the BranchPoint's dividing MA into parts each + * corresponding to a unique pair of opposite EDGEs + * \param [in,out] theMAParams - the MA division parameters to modify + * \param [in,out] theParBeg - index of the 1st division point for the given EDGE pair + * \param [in,out] theParEnd - index of the last division point for the given EDGE pair + * \return bool - is OK + */ + //================================================================================ + + bool getParamsForEdgePair( const size_t theEdgePairInd, + const vector< SMESH_MAT2d::BranchPoint >& theDivPoints, + const vector& theMAParams, + vector& theSelectedMAParams) + { + if ( theDivPoints.empty() ) + { + theSelectedMAParams = theMAParams; + return true; + } + if ( theEdgePairInd > theDivPoints.size() ) + return false; + + // TODO + return false; + } + + //-------------------------------------------------------------------------------- + // node or node parameter on EDGE + struct NodePoint + { + const SMDS_MeshNode* _node; + double _u; + int _edgeInd; // index in theSinuEdges vector + + NodePoint(const SMDS_MeshNode* n=0 ): _node(n), _u(0), _edgeInd(-1) {} + NodePoint(double u, size_t iEdge) : _node(0), _u(u), _edgeInd(iEdge) {} + NodePoint(const SMESH_MAT2d::BoundaryPoint& p) : _node(0), _u(p._param), _edgeInd(p._edgeIndex) {} + }; + + //================================================================================ + /*! + * \brief Finds a VERTEX corresponding to a point on EDGE, which is also filled + * with a node on the VERTEX, present or created + * \param [in,out] theNodePnt - the node position on the EDGE + * \param [in] theSinuEdges - the sinuous EDGEs + * \param [in] theMeshDS - the mesh + * \return bool - true if the \a theBndPnt is on VERTEX + */ + //================================================================================ + + bool findVertex( NodePoint& theNodePnt, + const vector& theSinuEdges, + SMESHDS_Mesh* theMeshDS) + { + if ( theNodePnt._edgeInd >= theSinuEdges.size() ) + return false; + + double f,l; + BRep_Tool::Range( theSinuEdges[ theNodePnt._edgeInd ], f,l ); + + TopoDS_Vertex V; + if ( Abs( f - theNodePnt._u )) + V = SMESH_MesherHelper::IthVertex( 0, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); + else if ( Abs( l - theNodePnt._u )) + V = SMESH_MesherHelper::IthVertex( 1, theSinuEdges[ theNodePnt._edgeInd ], /*CumOri=*/false); + + if ( !V.IsNull() ) + { + theNodePnt._node = SMESH_Algo::VertexNode( V, theMeshDS ); + if ( !theNodePnt._node ) + { + gp_Pnt p = BRep_Tool::Pnt( V ); + theNodePnt._node = theMeshDS->AddNode( p.X(), p.Y(), p.Z() ); + theMeshDS->SetNodeOnVertex( theNodePnt._node, V ); + } + return true; + } + return false; + } + + //================================================================================ + /*! + * \brief Add to the map of NodePoint's those on VERTEXes + * \param [in,out] theHelper - the helper + * \param [in] theMA - Medial Axis + * \param [in] theDivPoints - projections of VERTEXes to MA + * \param [in] theSinuEdges - the sinuous EDGEs + * \param [in] theSideEdgeIDs - indices of sinuous EDGEs per side + * \param [in] theIsEdgeComputed - is sinuous EGDE is meshed + * \param [in,out] thePointsOnE - the map to fill + */ + //================================================================================ + + bool projectVertices( SMESH_MesherHelper& theHelper, + const SMESH_MAT2d::MedialAxis& theMA, + const vector< SMESH_MAT2d::BranchPoint >& theDivPoints, + const vector& theSinuEdges, + //const vector< int > theSideEdgeIDs[2], + const vector< bool >& theIsEdgeComputed, + map< double, pair< NodePoint, NodePoint > > & thePointsOnE) + { + if ( theDivPoints.empty() ) + return true; + + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + + double uMA; + SMESH_MAT2d::BoundaryPoint bp[2]; + const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + + for ( size_t i = 0; i < theDivPoints.size(); ++i ) + { + if ( !branch.getParameter( theDivPoints[i], uMA )) + return false; + if ( !branch.getBoundaryPoints( theDivPoints[i], bp[0], bp[1] )) + return false; + + NodePoint np[2] = { NodePoint( bp[0] ), + NodePoint( bp[1] ) }; + bool isVertex[2] = { findVertex( np[0], theSinuEdges, meshDS ), + findVertex( np[1], theSinuEdges, meshDS )}; + + map< double, pair< NodePoint, NodePoint > >::iterator u2NP = + thePointsOnE.insert( make_pair( uMA, make_pair( np[0], np[1]))).first; + + if ( isVertex[0] && isVertex[1] ) + continue; + + bool isOppComputed = theIsEdgeComputed[ np[ isVertex[0] ]._edgeInd ]; + + if ( !isOppComputed ) + continue; + + // a VERTEX is projected on a meshed EDGE; there are two options: + // - a projected point is joined with a closet node if a strip between this and neighbor + // projection is wide enough; joining is done by setting the same node to the BoundaryPoint + // - a neighbor projection is merged this this one if it too close; a node of deleted + // projection is set to the BoundaryPoint of this projection + + + } + return true; + } + + //================================================================================ + /*! + * \brief Divide the sinuous EDGEs by projecting the division point of Medial + * Axis to the EGDEs + * \param [in] theHelper - the helper + * \param [in] theMA - the Medial Axis + * \param [in] theMAParams - parameters of division points of \a theMA + * \param [in] theSinuEdges - the EDGEs to make nodes on + * \param [in] theSinuSide0Size - the number of EDGEs in the 1st sinuous side + * \return bool - is OK or not + */ + //================================================================================ + + bool computeSinuEdges( SMESH_MesherHelper& theHelper, + SMESH_MAT2d::MedialAxis& theMA, + vector& theMAParams, + const vector& theSinuEdges, + const size_t theSinuSide0Size) + { + if ( theMA.getBranches().size() != 1 ) + return false; + + // normalize theMAParams + for ( size_t i = 0; i < theMAParams.size(); ++i ) + theMAParams[i] /= theMAParams.back(); + + + SMESH_Mesh* mesh = theHelper.GetMesh(); + SMESHDS_Mesh* meshDS = theHelper.GetMeshDS(); + double f,l; + + vector< Handle(Geom_Curve) > curves ( theSinuEdges.size() ); + vector< int > edgeIDs( theSinuEdges.size() ); + vector< bool > isComputed( theSinuEdges.size() ); + //bool hasComputed = false; + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) + { + curves[i] = BRep_Tool::Curve( theSinuEdges[i], f,l ); + if ( !curves[i] ) + return false; + SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); + edgeIDs [i] = sm->GetId(); + isComputed[i] = ( !sm->IsEmpty() ); + // if ( isComputed[i] ) + // hasComputed = true; + } + + const SMESH_MAT2d::Branch& branch = theMA.getBranches()[0]; + SMESH_MAT2d::BoundaryPoint bp[2]; + + vector< std::size_t > edgeIDs1, edgeIDs2; + vector< SMESH_MAT2d::BranchPoint > divPoints; + branch.getOppositeGeomEdges( edgeIDs1, edgeIDs2, divPoints ); + for ( size_t i = 0; i < edgeIDs1.size(); ++i ) + if ( isComputed[ edgeIDs1[i]] && + isComputed[ edgeIDs2[i]]) + return false; + + // map param on MA to parameters of nodes on a pair of theSinuEdges + typedef map< double, pair< NodePoint, NodePoint > > TMAPar2NPoints; + TMAPar2NPoints pointsOnE; + vector maParams; + + // compute params of nodes on EDGEs by projecting division points from MA + //const double tol = 1e-5 * theMAParams.back(); + size_t iEdgePair = 0; + while ( iEdgePair < edgeIDs1.size() ) + { + if ( isComputed[ edgeIDs1[ iEdgePair ]] || + isComputed[ edgeIDs2[ iEdgePair ]]) + { + // "projection" from one side to the other + + size_t iEdgeComputed = edgeIDs1[iEdgePair], iSideComputed = 0; + if ( !isComputed[ iEdgeComputed ]) + ++iSideComputed, iEdgeComputed = edgeIDs2[iEdgePair]; + + map< double, const SMDS_MeshNode* > nodeParams; // params of existing nodes + if ( !SMESH_Algo::GetSortedNodesOnEdge( meshDS, theSinuEdges[ iEdgeComputed ], /*skipMedium=*/true, nodeParams )) + return false; + + SMESH_MAT2d::BoundaryPoint& bndPnt = bp[ 1-iSideComputed ]; + SMESH_MAT2d::BranchPoint brp; + NodePoint npN, npB; + NodePoint& np0 = iSideComputed ? npB : npN; + NodePoint& np1 = iSideComputed ? npN : npB; + + double maParam1st, maParamLast, maParam; + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.begin()->first, brp )) + return false; + branch.getParameter( brp, maParam1st ); + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, nodeParams.rbegin()->first, brp )) + return false; + branch.getParameter( brp, maParamLast ); + + map< double, const SMDS_MeshNode* >::iterator u2n = nodeParams.begin(), u2nEnd = --nodeParams.end(); + TMAPar2NPoints::iterator pos, end = pointsOnE.end(); + TMAPar2NPoints::iterator & hint = (maParamLast > maParam1st) ? end : pos; + for ( ++u2n; u2n != u2nEnd; ++u2n ) + { + if ( !theMA.getBoundary().getBranchPoint( iEdgeComputed, u2n->first, brp )) + return false; + if ( !branch.getBoundaryPoints( brp, bp[0], bp[1] )) + return false; + if ( !branch.getParameter( brp, maParam )) + return false; + + npN = NodePoint( u2n->second ); + npB = NodePoint( bndPnt ); + pos = pointsOnE.insert( hint, make_pair( maParam, make_pair( np0, np1 ))); + } + + // move iEdgePair forward + while ( iEdgePair < edgeIDs1.size() ) + if ( edgeIDs1[ iEdgePair ] == bp[0]._edgeIndex && + edgeIDs2[ iEdgePair ] == bp[1]._edgeIndex ) + break; + else + ++iEdgePair; + } + else + { + // projection from MA + maParams.clear(); + if ( !getParamsForEdgePair( iEdgePair, divPoints, theMAParams, maParams )) + return false; + + for ( size_t i = 1; i < maParams.size()-1; ++i ) + { + if ( !branch.getBoundaryPoints( maParams[i], bp[0], bp[1] )) + return false; + + pointsOnE.insert( pointsOnE.end(), make_pair( maParams[i], make_pair( NodePoint(bp[0]), + NodePoint(bp[1])))); + } + } + ++iEdgePair; + } + + if ( !projectVertices( theHelper, theMA, divPoints, theSinuEdges, isComputed, pointsOnE )) + return false; + + // create nodes + TMAPar2NPoints::iterator u2np = pointsOnE.begin(); + for ( ; u2np != pointsOnE.end(); ++u2np ) + { + NodePoint* np[2] = { & u2np->second.first, & u2np->second.second }; + for ( int iSide = 0; iSide < 2; ++iSide ) + { + if ( np[ iSide ]->_node ) continue; + size_t iEdge = np[ iSide ]->_edgeInd; + double u = np[ iSide ]->_u; + gp_Pnt p = curves[ iEdge ]->Value( u ); + np[ iSide ]->_node = meshDS->AddNode( p.X(), p.Y(), p.Z() ); + meshDS->SetNodeOnEdge( np[ iSide ]->_node, edgeIDs[ iEdge ], u ); + } + } + + // create mesh segments on EDGEs + theHelper.SetElementsOnShape( false ); + TopoDS_Face face = TopoDS::Face( theHelper.GetSubShape() ); + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) + { + SMESH_subMesh* sm = mesh->GetSubMesh( theSinuEdges[i] ); + if ( sm->GetSubMeshDS() && sm->GetSubMeshDS()->NbElements() > 0 ) + continue; + + StdMeshers_FaceSide side( face, theSinuEdges[i], mesh, + /*isFwd=*/true, /*skipMediumNodes=*/true ); + vector nodes = side.GetOrderedNodes(); + for ( size_t in = 1; in < nodes.size(); ++in ) + { + const SMDS_MeshElement* seg = theHelper.AddEdge( nodes[in-1], nodes[in], 0, false ); + meshDS->SetMeshElementOnShape( seg, edgeIDs[ i ] ); + } + } + + // update sub-meshes on VERTEXes + for ( size_t i = 0; i < theSinuEdges.size(); ++i ) + { + mesh->GetSubMesh( theHelper.IthVertex( 0, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + mesh->GetSubMesh( theHelper.IthVertex( 1, theSinuEdges[i] )) + ->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + } + + return true; + } + + //================================================================================ + /*! + * \brief Mesh short EDGEs + */ + //================================================================================ + + bool computeShortEdges( SMESH_MesherHelper& theHelper, + const vector& theShortEdges, + SMESH_Algo* the1dAlgo ) + { + for ( size_t i = 0; i < theShortEdges.size(); ++i ) + { + theHelper.GetGen()->Compute( *theHelper.GetMesh(), theShortEdges[i], true, true ); + + SMESH_subMesh* sm = theHelper.GetMesh()->GetSubMesh(theShortEdges[i] ); + if ( sm->IsEmpty() ) + { + try { + if ( !the1dAlgo->Compute( *theHelper.GetMesh(), theShortEdges[i] )) + return false; + } + catch (...) { + return false; + } + sm->ComputeStateEngine( SMESH_subMesh::CHECK_COMPUTE_STATE ); + if ( sm->IsEmpty() ) + return false; + } + } + return true; + } + + inline double area( const UVPtStruct& p1, const UVPtStruct& p2, const UVPtStruct& p3 ) + { + gp_XY v1 = p2.UV() - p1.UV(); + gp_XY v2 = p3.UV() - p1.UV(); + return v2 ^ v1; + } + + bool ellipticSmooth( FaceQuadStruct::Ptr quad, int nbLoops ) + { + //nbLoops = 10; + if ( quad->uv_grid.empty() ) + return true; + + int nbhoriz = quad->iSize; + int nbvertic = quad->jSize; + + const double dksi = 0.5, deta = 0.5; + const double dksi2 = dksi*dksi, deta2 = deta*deta; + double err = 0., g11, g22, g12; + int nbErr = 0; + + FaceQuadStruct& q = *quad; + UVPtStruct pNew; + + double refArea = area( q.UVPt(0,0), q.UVPt(1,0), q.UVPt(1,1) ); + + for ( int iLoop = 0; iLoop < nbLoops; ++iLoop ) + { + err = 0; + for ( int i = 1; i < nbhoriz - 1; i++ ) + for ( int j = 1; j < nbvertic - 1; j++ ) + { + g11 = ( (q.U(i,j+1) - q.U(i,j-1))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 + + (q.V(i,j+1) - q.V(i,j-1))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/4; + + g22 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i+1,j) - q.U(i-1,j))/dksi2 + + (q.V(i+1,j) - q.V(i-1,j))*(q.V(i+1,j) - q.V(i-1,j))/deta2 )/4; + + g12 = ( (q.U(i+1,j) - q.U(i-1,j))*(q.U(i,j+1) - q.U(i,j-1))/dksi2 + + (q.V(i+1,j) - q.V(i-1,j))*(q.V(i,j+1) - q.V(i,j-1))/deta2 )/(4*dksi*deta); + + pNew.u = dksi2/(2*(g11+g22)) * (g11*(q.U(i+1,j) + q.U(i-1,j))/dksi2 + + g22*(q.U(i,j+1) + q.U(i,j-1))/dksi2 + - 0.5*g12*q.U(i+1,j+1) + 0.5*g12*q.U(i-1,j+1) + + - 0.5*g12*q.U(i-1,j-1) + 0.5*g12*q.U(i+1,j-1)); + + pNew.v = deta2/(2*(g11+g22)) * (g11*(q.V(i+1,j) + q.V(i-1,j))/deta2 + + g22*(q.V(i,j+1) + q.V(i,j-1))/deta2 + - 0.5*g12*q.V(i+1,j+1) + 0.5*g12*q.V(i-1,j+1) + + - 0.5*g12*q.V(i-1,j-1) + 0.5*g12*q.V(i+1,j-1)); + + // if (( refArea * area( q.UVPt(i-1,j-1), q.UVPt(i,j-1), pNew ) > 0 ) && + // ( refArea * area( q.UVPt(i+1,j-1), q.UVPt(i+1,j), pNew ) > 0 ) && + // ( refArea * area( q.UVPt(i+1,j+1), q.UVPt(i,j+1), pNew ) > 0 ) && + // ( refArea * area( q.UVPt(i-1,j), q.UVPt(i-1,j-1), pNew ) > 0 )) + { + err += sqrt(( q.U(i,j) - pNew.u ) * ( q.U(i,j) - pNew.u ) + + ( q.V(i,j) - pNew.v ) * ( q.V(i,j) - pNew.v )); + q.U(i,j) = pNew.u; + q.V(i,j) = pNew.v; + } + // else if ( ++nbErr < 10 ) + // { + // cout << i << ", " << j << endl; + // cout << "x = [" + // << "[ " << q.U(i-1,j-1) << ", " < theSinuEdges[2], + const vector theShortEdges[2]) +{ + SMESH_Mesh* mesh = theHelper.GetMesh(); + SMESH_ProxyMesh::Ptr proxyMesh = StdMeshers_ViscousLayers2D::Compute( *mesh, theFace ); + if ( !proxyMesh ) + return false; + + StdMeshers_Quadrangle_2D::myProxyMesh = proxyMesh; + StdMeshers_Quadrangle_2D::myHelper = &theHelper; + StdMeshers_Quadrangle_2D::myNeedSmooth = false; + StdMeshers_Quadrangle_2D::myCheckOri = false; + StdMeshers_Quadrangle_2D::myQuadList.clear(); + + // fill FaceQuadStruct + + list< TopoDS_Edge > side[4]; + side[0].insert( side[0].end(), theShortEdges[0].begin(), theShortEdges[0].end() ); + side[1].insert( side[1].end(), theSinuEdges[1].begin(), theSinuEdges[1].end() ); + side[2].insert( side[2].end(), theShortEdges[1].begin(), theShortEdges[1].end() ); + side[3].insert( side[3].end(), theSinuEdges[0].begin(), theSinuEdges[0].end() ); + + FaceQuadStruct::Ptr quad( new FaceQuadStruct ); + quad->side.resize( 4 ); + quad->face = theFace; + for ( int i = 0; i < 4; ++i ) + { + quad->side[i] = StdMeshers_FaceSide::New( theFace, side[i], mesh, i < QUAD_TOP_SIDE, + /*skipMediumNodes=*/true, proxyMesh ); + } + int nbNodesShort0 = quad->side[0].NbPoints(); + int nbNodesShort1 = quad->side[2].NbPoints(); + + // compute UV of internal points + myQuadList.push_back( quad ); + if ( !StdMeshers_Quadrangle_2D::setNormalizedGrid( quad )) + return false; + + // elliptic smooth of internal points to get boundary cell normal to the boundary + ellipticSmooth( quad, 1 ); + + // create quadrangles + bool ok; + if ( nbNodesShort0 == nbNodesShort1 ) + ok = StdMeshers_Quadrangle_2D::computeQuadDominant( *mesh, theFace, quad ); + else + ok = StdMeshers_Quadrangle_2D::computeTriangles( *mesh, theFace, quad ); + + StdMeshers_Quadrangle_2D::myProxyMesh.reset(); + StdMeshers_Quadrangle_2D::myHelper = 0; + + return ok; +} + +//================================================================================ +/*! + * \brief Generate quadrangle mesh + */ +//================================================================================ + +bool StdMeshers_QuadFromMedialAxis_1D2D::Compute(SMESH_Mesh& theMesh, + const TopoDS_Shape& theShape) +{ + SMESH_MesherHelper helper( theMesh ); + helper.SetSubShape( theShape ); + + TopoDS_Face F = TopoDS::Face( theShape ); + if ( F.Orientation() >= TopAbs_INTERNAL ) F.Orientation( TopAbs_FORWARD ); + + list< TopoDS_Edge > edges; + list< int > nbEdgesInWire; + int nbWire = SMESH_Block::GetOrderedEdges (F, edges, nbEdgesInWire); + + vector< TopoDS_Edge > sinuSide[2], shortSide[2]; + if ( nbWire == 1 && getSinuousEdges( helper, F, edges, sinuSide, shortSide )) + { + vector< TopoDS_Edge > sinuEdges = sinuSide[0]; + sinuEdges.insert( sinuEdges.end(), sinuSide[1].begin(), sinuSide[1].end() ); + if ( sinuEdges.size() > 2 ) + return error(COMPERR_BAD_SHAPE, "Not yet supported case" ); + + double minSegLen = getMinSegLen( helper, sinuEdges ); + SMESH_MAT2d::MedialAxis ma( F, sinuEdges, minSegLen, /*ignoreCorners=*/true ); + + if ( !_regular1D ) + _regular1D = new Algo1D( _studyId, _gen ); + _regular1D->SetSegmentLength( minSegLen ); + + vector maParams; + if ( ! divideMA( helper, ma, sinuEdges, sinuSide[0].size(), _regular1D, maParams )) + return error(COMPERR_BAD_SHAPE); + + if ( !computeShortEdges( helper, shortSide[0], _regular1D ) || + !computeShortEdges( helper, shortSide[1], _regular1D )) + return error("Failed to mesh short edges"); + + if ( !computeSinuEdges( helper, ma, maParams, sinuEdges, sinuSide[0].size() )) + return error("Failed to mesh sinuous edges"); + + return computeQuads( helper, F, sinuSide, shortSide ); + } + + return error(COMPERR_BAD_SHAPE, "Not implemented so far"); +} + +//================================================================================ +/*! + * \brief Predict nb of elements + */ +//================================================================================ + +bool StdMeshers_QuadFromMedialAxis_1D2D::Evaluate(SMESH_Mesh & theMesh, + const TopoDS_Shape & theShape, + MapShapeNbElems& theResMap) +{ + return StdMeshers_Quadrangle_2D::Evaluate(theMesh,theShape,theResMap); +} + diff --git a/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx new file mode 100644 index 000000000..df8a84a4c --- /dev/null +++ b/src/StdMeshers/StdMeshers_QuadFromMedialAxis_1D2D.hxx @@ -0,0 +1,65 @@ +// Copyright (C) 2007-2015 CEA/DEN, EDF R&D, OPEN CASCADE +// +// Copyright (C) 2003-2007 OPEN CASCADE, EADS/CCR, LIP6, CEA/DEN, +// CEDRAT, EDF R&D, LEG, PRINCIPIA R&D, BUREAU VERITAS +// +// This library is free software; you can redistribute it and/or +// modify it under the terms of the GNU Lesser General Public +// License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// +// This library is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU +// Lesser General Public License for more details. +// +// You should have received a copy of the GNU Lesser General Public +// License along with this library; if not, write to the Free Software +// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +// +// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com +// +// File : StdMeshers_QuadFromMedialAxis_1D2D.hxx +// Created : Wed Jun 3 17:22:35 2015 +// Author : Edward AGAPOV (eap) + + +#ifndef __StdMeshers_QuadFromMedialAxis_1D2D_HXX__ +#define __StdMeshers_QuadFromMedialAxis_1D2D_HXX__ + +#include "StdMeshers_Quadrangle_2D.hxx" + +#include + +/*! + * \brief Quadrangle mesher using Medial Axis + */ +class STDMESHERS_EXPORT StdMeshers_QuadFromMedialAxis_1D2D: public StdMeshers_Quadrangle_2D +{ + public: + StdMeshers_QuadFromMedialAxis_1D2D(int hypId, int studyId, SMESH_Gen* gen); + virtual ~StdMeshers_QuadFromMedialAxis_1D2D(); + + virtual bool CheckHypothesis(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape, + Hypothesis_Status& aStatus); + + virtual bool Compute(SMESH_Mesh& aMesh, + const TopoDS_Shape& aShape); + + virtual bool Evaluate(SMESH_Mesh & aMesh, + const TopoDS_Shape & aShape, + MapShapeNbElems& aResMap); + + private: + + bool computeQuads( SMESH_MesherHelper& theHelper, + const TopoDS_Face& theFace, + const std::vector theSinuEdges[2], + const std::vector theShortEdges[2]); + + class Algo1D; + Algo1D* _regular1D; +}; + +#endif diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx index af3d51d76..0e35c560e 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.cxx @@ -1599,7 +1599,7 @@ void StdMeshers_Quadrangle_2D::shiftQuad(FaceQuadStruct::Ptr& quad, const int nu //================================================================================ /*! - * \brief Rotate sides of a quad by given nb of quartes + * \brief Rotate sides of a quad CCW by given nb of quartes * \param nb - number of rotation quartes * \param ori - to keep orientation of sides as in an unit quad or not * \param keepGrid - if \c true Side::grid is not changed, Side::from and Side::to @@ -1611,6 +1611,8 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) { if ( nb == 0 ) return; + nb = nb % NB_QUAD_SIDES; + vector< Side > newSides( side.size() ); vector< Side* > sidePtrs( side.size() ); for (int i = QUAD_BOTTOM_SIDE; i < NB_QUAD_SIDES; ++i) @@ -1640,7 +1642,33 @@ void FaceQuadStruct::shift( size_t nb, bool ori, bool keepGrid ) } newSides.swap( side ); - uv_grid.clear(); + if ( keepGrid && !uv_grid.empty() ) + { + if ( nb == 2 ) // "PI" + { + std::reverse( uv_grid.begin(), uv_grid.end() ); + } + else + { + FaceQuadStruct newQuad; + newQuad.uv_grid.resize( uv_grid.size() ); + newQuad.iSize = jSize; + newQuad.jSize = iSize; + int i, j, iRev, jRev; + int *iNew = ( nb == 1 ) ? &jRev : &j; + int *jNew = ( nb == 1 ) ? &i : &iRev; + for ( i = 0, iRev = iSize-1; i < iSize; ++i, --iRev ) + for ( j = 0, jRev = jSize-1; j < jSize; ++j, --jRev ) + newQuad.UVPt( *iNew, *jNew ) = UVPt( i, j ); + + std::swap( iSize, jSize ); + std::swap( uv_grid, newQuad.uv_grid ); + } + } + else + { + uv_grid.clear(); + } } //======================================================================= diff --git a/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx b/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx index 8687bbfda..6a06c9a46 100644 --- a/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx +++ b/src/StdMeshers/StdMeshers_Quadrangle_2D.hxx @@ -119,6 +119,8 @@ struct FaceQuadStruct FaceQuadStruct ( const TopoDS_Face& F = TopoDS_Face(), const std::string& nm="main" ); UVPtStruct& UVPt( int i, int j ) { return uv_grid[ i + j * iSize ]; } + double& U( int i, int j ) { return UVPt( i, j ).u; } + double& V( int i, int j ) { return UVPt( i, j ).v; } void shift ( size_t nb, bool keepUnitOri, bool keepGrid=false ); int & nbNodeOut( int iSide ) { return side[ iSide ].nbNodeOut; } bool findCell ( const gp_XY& uv, int & i, int & j ); diff --git a/src/StdMeshersGUI/StdMeshers_images.ts b/src/StdMeshersGUI/StdMeshers_images.ts index b94521ab3..e744cf34e 100644 --- a/src/StdMeshersGUI/StdMeshers_images.ts +++ b/src/StdMeshersGUI/StdMeshers_images.ts @@ -143,6 +143,10 @@ ICON_SMESH_TREE_ALGO_RadialQuadrangle_1D2D mesh_tree_algo_radial_quadrangle_1D2D.png + + ICON_SMESH_TREE_ALGO_QuadFromMedialAxis_1D2D + mesh_tree_algo_quad.png + ICON_SMESH_TREE_ALGO_Prism_3D mesh_tree_algo_prism.png diff --git a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx index 2af7cca75..cf497681c 100644 --- a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx +++ b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.cxx @@ -32,6 +32,8 @@ #include "Utils_CorbaException.hxx" #include "utilities.h" +#include "StdMeshers_QuadFromMedialAxis_1D2D.hxx" + using namespace std; //============================================================================= @@ -98,3 +100,39 @@ CORBA::Boolean StdMeshers_Quadrangle_2D_i::IsApplicable( const TopoDS_Shape &S, return ::StdMeshers_Quadrangle_2D::IsApplicable( S, toCheckAll ); } +//============================================================================= +/*! + * StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i + * + * Constructor + */ +//============================================================================= + +StdMeshers_QuadFromMedialAxis_1D2D_i:: +StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA, + int theStudyId, + ::SMESH_Gen* theGenImpl ) + : SALOME::GenericObj_i( thePOA ), + SMESH_Hypothesis_i( thePOA ), + SMESH_Algo_i( thePOA ), + SMESH_2D_Algo_i( thePOA ) +{ + MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::StdMeshers_QuadFromMedialAxis_1D2D_i" ); + myBaseImpl = new ::StdMeshers_QuadFromMedialAxis_1D2D( theGenImpl->GetANewId(), + theStudyId, + theGenImpl ); +} + +//============================================================================= +/*! + * StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i + * + * Destructor + * + */ +//============================================================================= + +StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i() +{ + MESSAGE( "StdMeshers_QuadFromMedialAxis_1D2D_i::~StdMeshers_QuadFromMedialAxis_1D2D_i" ); +} diff --git a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx index 4fa55b135..145d31e3c 100644 --- a/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx +++ b/src/StdMeshers_I/StdMeshers_Quadrangle_2D_i.hxx @@ -62,4 +62,27 @@ class STDMESHERS_I_EXPORT StdMeshers_Quadrangle_2D_i: static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll); }; +// ====================================================== +// Quadrangle (Medial Axis Projection) 2d algorithm +// ====================================================== +class STDMESHERS_I_EXPORT StdMeshers_QuadFromMedialAxis_1D2D_i: + public virtual POA_StdMeshers::StdMeshers_QuadFromMedialAxis_1D2D, + public virtual SMESH_2D_Algo_i +{ + public: + // Constructor + StdMeshers_QuadFromMedialAxis_1D2D_i( PortableServer::POA_ptr thePOA, + int theStudyId, + ::SMESH_Gen* theGenImpl ); + + // Destructor + virtual ~StdMeshers_QuadFromMedialAxis_1D2D_i(); + + // Get implementation + //::StdMeshers_Quadrangle_2D* GetImpl(); + + // Return true if the algorithm is applicable to a shape + //static CORBA::Boolean IsApplicable(const TopoDS_Shape &S, CORBA::Boolean toCheckAll); +}; + #endif diff --git a/src/StdMeshers_I/StdMeshers_i.cxx b/src/StdMeshers_I/StdMeshers_i.cxx index c7ae1a76a..3e9eb70b7 100644 --- a/src/StdMeshers_I/StdMeshers_i.cxx +++ b/src/StdMeshers_I/StdMeshers_i.cxx @@ -216,6 +216,8 @@ STDMESHERS_I_EXPORT #endif else if (strcmp(aHypName, "Quadrangle_2D") == 0) aCreator = new StdHypothesisCreator_i; + else if (strcmp(aHypName, "QuadFromMedialAxis_1D2D") == 0) + aCreator = new StdHypothesisCreator_i; else if (strcmp(aHypName, "Hexa_3D") == 0) aCreator = new StdHypothesisCreator_i; else if (strcmp(aHypName, "Projection_1D") == 0) -- 2.39.2