From f35a157d9ca15dff56cb76024c785bfb879d7f38 Mon Sep 17 00:00:00 2001 From: nivesh Date: Tue, 29 May 2018 11:21:23 +0200 Subject: [PATCH] added Helper function to create FECollection for overlapping sub-domains with individual enrichment functions --- doc/doxygen/images/3source_fe_indices.png | Bin 0 -> 15043 bytes doc/news/changes/major/20180420NiveshD | 9 + include/deal.II/fe/fe_enriched.h | 314 +++ include/deal.II/fe/fe_enriched.templates.h | 212 ++ source/fe/fe_enriched.cc | 370 +++ source/fe/fe_enriched.inst.in | 55 + tests/fe/fe_enriched_color_01.cc | 148 ++ tests/fe/fe_enriched_color_01.output | 53 + tests/fe/fe_enriched_color_02.cc | 163 ++ .../fe_enriched_color_02.with_mpi=true.output | 13 + tests/fe/fe_enriched_color_03.cc | 175 ++ .../fe_enriched_color_03.with_mpi=true.output | 17 + tests/fe/fe_enriched_color_04.cc | 176 ++ .../fe_enriched_color_04.with_mpi=true.output | 18 + tests/fe/fe_enriched_color_05.cc | 264 +++ .../fe_enriched_color_05.with_mpi=true.output | 23 + tests/fe/fe_enriched_color_06.cc | 142 ++ ...olor_06.with_gsl=true.with_mpi=true.output | 18 + tests/fe/fe_enriched_color_07.cc | 1988 +++++++++++++++++ ...ith_gsl=true.with_mpi=true.mpirun=2.output | 14 + 20 files changed, 4172 insertions(+) create mode 100644 doc/doxygen/images/3source_fe_indices.png create mode 100644 doc/news/changes/major/20180420NiveshD create mode 100644 include/deal.II/fe/fe_enriched.templates.h create mode 100644 tests/fe/fe_enriched_color_01.cc create mode 100644 tests/fe/fe_enriched_color_01.output create mode 100644 tests/fe/fe_enriched_color_02.cc create mode 100644 tests/fe/fe_enriched_color_02.with_mpi=true.output create mode 100644 tests/fe/fe_enriched_color_03.cc create mode 100644 tests/fe/fe_enriched_color_03.with_mpi=true.output create mode 100644 tests/fe/fe_enriched_color_04.cc create mode 100644 tests/fe/fe_enriched_color_04.with_mpi=true.output create mode 100644 tests/fe/fe_enriched_color_05.cc create mode 100644 tests/fe/fe_enriched_color_05.with_mpi=true.output create mode 100644 tests/fe/fe_enriched_color_06.cc create mode 100644 tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output create mode 100644 tests/fe/fe_enriched_color_07.cc create mode 100644 tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output diff --git a/doc/doxygen/images/3source_fe_indices.png b/doc/doxygen/images/3source_fe_indices.png new file mode 100644 index 0000000000000000000000000000000000000000..f07c0b9c2aab0445f34e0d1b5bc964ccfb9b5652 GIT binary patch literal 15043 zcmeHuXH-;Kv-WO+m@o_?0%AKuR75~fL_+rrqYNNHQIU*-fFObd5lOAi7zs+2tSE>G zsN~cQI)Y>rNllg{5+&13ci%dI&dmGHyVkuw?poja?(n0}>C-#Zu3h!iQ?(Db^mNp@ zR|~902ytut^7B!IIN+AuXXOgG+zIk1gbRnuPa4No!k^bltD6XIM;br>cf=r{AP)T}!nny3yL9 z)qN>momwxMt!^BcbnLT3O3|mCXQg9}7H~gu_-`i!-TnFM6)QIOpZ1hj4b!=-BFcvk ziv`Cp+9w%gGyCeR3@=JBqlH~buGUkvb~3ZUF3SicBptPv&W9`Nddq`dh@302OFSF4 z7^H_|!!Az{4!ay9VV`Ixaq+^t9lMZLZ@@19m%7Cqc7s}s-nY9;@6P8a zeZY<>;{UpG{>5G_ZeN-x7k4S)qxy*c-H5{JYmw0>1i8xM^5QI{AP(i`FzLNckezlJ7PYK|hgh z7UM6d(BI)XR^688Qs3GdC7E96Jl9cfQZSO;8c8#P-LWmqHcCAqnZI(~c?*l}LPA1k z?d>bs*DK$x++DCMP-Uh$f`0k(<(C&nsL{?9J6>BxDLB@=JX86;rA0B>zBR{@@KXBx zZ9}q`_d8KAa;<#&=$nJ##rrvs{##w2iQt7-(hEP4Cu7EXQ_LB5VbbQR{ju5E*=yIW zODtkcHwaz3c5QHYcq?;Ot*GDk-aE{v^x{x_R;b;o;(m8Z!Kd#^-Dn2nJrC&X@XeJ8yNk&xQ4oV*(r0dF$+7%&R zhYgKJd*9Fy$}4%I_qyQuUknVI9;InV?Mppke)Z_{RofBzV8c_~c;YGQ>UJk*fGV@5 z?9CHJV<9jEs}_;gqBIdBRr-xg`iACbN-ni?<(+KSl8!FD}84#1=Lrgmuyj z$3D@{_+vPd3U(JyKFWy>EZ(X9F#z+*m4{mJkQ8?8DI#)VKYtL%K1r)343NzJldh}a zOI?sAFWSEX!^h9CWfs`8Z|(kyn$U1+{n|~YnD=Bp`ZgumnG;>Rc5qQBspMn0V$wQn zj3+4|9@yieF=3alpriK`6>nR7Zl3&gMRC~oKevd_DEfp2t@yp6lIW8}J~pSV_LAR{ zV8Sr0J<;^>YLxH#?RFQocoJG4;vE)Nnl%(Oic zw%at)_jMM#&nYngEsvde`m8)#pi~U3_y=I&#|Hwm}WF@Al7R#s-ysVi z$aSOtM3H+ouj#8?ADVqI{b-!)dkw5denCO7h;-rjPdWs~RN?rmH8;G*9*3JnV$gO^ zeVc}^u5MbvXugT}bi;juuTNu|X_UT>;+eTxYN3Zqn3QqYNNbMn)6thk{{~W;SzJUz z1zcO^S$*rO2-$_X5-xmIkf@fJ+ZX_-bGc}EFN|D~s2%0pawhz7L4i_hwnODr&NYV^ zhT?Hv&KHN9-*p)P;t)4-Angy0`(CX3jS;oykj^LHx!3r9(1@PuY+h{xqwQKdz3bX~ zYij0MjUXPK{oltir;;-48lH-1c+x-delH{>#pL7!K z4s6XGRQ!z|40U%?7ACv~wahDi<1^b2I}ZTk1@dik8|w)q#QwI1&$^?~qXj_yiTzMx zs+O*9Z8tBt5m5_VwT?}Rfxo4Xm=n^hJqm43i~CnNU+uWfh4C~cEwZkN&+O8v#&q*I zrO{mBEF&qUi^Fd#LM5@qnOu0oYnbPJ=CY?}aN*_XE^USO5WuR#WG9f(=uol_?$Y^# z_Hk@ZA3O21@XG@=G4}y-0b=*qH$FSk)R^w;ra*&W@xd-@ReOQkP@{=U^522myi)_;9MjC4{alyweP#Z4WI=0uk|KV>^?IQuww}+ z)8MY}G8O>pWcnia3}HW$@+){z$WQ~^5?_bdSD2(>mW8L$sf*3tomDF{@iN`kdOamwf*;j-N9dK@*lwU*OgaJ+rjVg zrz-ZI$zH_6&pQo}-k+AUjM4u|n5>6u_TS82RCWAQ4EU${OiSkatCo=X>C$Q^^-lC@QMN_3FCj9f)lZnbqbkCeq{N)(=6 zpE*<3O=~Vw4SREkFDW(Ep^d)nzYm)|W*8TNxwwce;i%1PvN~_RKWe=IYIE=O&STPs zCcgywdVQG7YBG!Zcbc{b=!SL>yBWsR0-Eh1`ot?zAtOos^n*a3?ik+pv`J>DqkUC) zerl>z55uNl>b+LRg`QI7iMPUogM;M&ui+RcoE)gFrDoZ?RP0Dv+gac-*}uy;dsj?M z3|C>cH+|;%2h+`!;px)(uurxc?wv)spT662V;$)`BOp7w^Lm5YV4N1f$K3_PFJ0Z; zZ2$lbpIekX@zhnF?;UP_A*@ka1b$HbK@8bGyp}aO=8~(;zg~6D0zflNM6vrZ*;eKb#pvkXcvo3+*Iy2fH1$eYB5D3>UUF9!7 z2U$5eIeqw@OkqtY@RUNQ;1iyrp0kA!53&>F9H`GSXRfLsiS(nKD76Duw>D6fW~(}# ztlL*<&YT!;{t9pOnUO_!bFyg9cf*uZS67Ei#}z-|e;dok{9f0Ab~<*1VRq4ff!TJ! zIl$OSyW>*n?XB&GPha0e@qMA6l^oxdEeH)-qV>H1dIHSZHPVGLLrv+mLhfz%fRG6( zU9{?}2sO;IJD!&N@%A&5ygfxlMH)c>No}t6+J-MqzsYR^bhq>NV3xC_%qgs%mZ@Q% zV3r8K*8!!R1)ovRao^s)2S4rPAXV%HGK$Yl$%Q5WCQ78OR-2WirR=zbaR&=^6AOoT)|rlG>Uvlx3f>nWfKU zJbcq)O!@@?1fxy-$8nH?nZ5 zt=#C(PMieHZFb-xB8g-4XIjoR^WlXvwz@NgQzPeEFIir8?}&@3C4}7@s0?|sY9mVg z6*eDmp0rseb8itLUqw)XVR|u}#$SZ=z6bDK@WY|s1gO(}$`i%vlT~T=62z@)G$=u_ zJu4;?)LA>F`HpWz$4v4T*`QuDHEiLNCq=<-f?&_Y_a*t-9Sj)#{ENdZ-{j#74LD-o zGM(FGJNCMzQQ@>HtIGc6lr|-{g}b@mw>nZTNQBRPq4!A{eT#zG;|rwnolCTK(^*Q- zB_-=5%ir}J;uXlp*@+jZcB0g)vmaRRFOEm8&OCD*({Z-fJPmlWt+v`wMz`@;1+Dua zv#j@5Ime1Ym+KSnGkk-!jme?y!so}LXX7Rs6?SPf`Mnv8Z9GmaLGsk)H2cWYbqKh3 zJaNSFZL&>~g0AmCO_LZ++KSH8WJSlF$6uPEwXx3ap+EDjo=@7{O!2j7W{$8>aW5$T z(hFACEJEsR}%zVk3gH5aS z?`&9W8v?rkVUqOF&5HWMpm#XO^j^l0iOP(FuJw=gfL=jCoan{x81=@2L7gJTut_Tf z(?>k2#mVizYj(ZOyEHuyME)2F18bv*EEX{#1n4)U0mPPO#Ek?JAl_h+L-Wwe74n)( z8@3dems1lgJ`Sq$nc`5L1vM5IFALCr5`dFfzOsaLXotaKf)NtHU5DajvA9(C>Q{@! zNuGa;Qgw<6OA$a6s6WQGA193PFHNliJtX+Q-%=+5Rt4moIQX|xElYr9Db9Zm_n30t zznfW${EsY$98qBoc7exBi;@r$jK>izVZqp>-vq9c{TFc%as$Dzgd!v$+XmVCEzj{y zutdD6D<>jeYx!G^EoB!NCs|HhzAnYqa|Njs=w3X8EK<CX)NEkJJX{7-im9UTtS@=M#=lL3y**1}K@7eb7p>&EWo;6pd?%e@_*t2-x zmp$DBTqynp2iOKy13b~G==h52GoJ8JYm33iu)f%|hl+}dVZlWQNP(tnJO4;B0_l0a z<6ORX2kVmtEb;V{LA~#_q^nE)1t;GM2R=$l5^pkfH%%kcK3ByNSftu_P(m3TkWT2U zY{D;yApz?`?z%>=kFC|axdLw9F^DsFT~)VXIV7u>A3#Q4JjEK+YzO1byH3j34fh7? z)&vMCl`;y$%`WeVV=W9!)~T)M4oo%9*&S};tTCEDl!hf};ieZ2A%$8`+ZycQgGjGF ztO&A+(_LXkH*LyH2;#BfHLe#?1mxR2G0BANfH;ZtaEV zU|$j@tnA(ZXC0={Z>!`rvEB7tRczIckwh${zdOkH_^%;dq}+x;Z}IAvi{mG=AQ!Ci z7>Eh4*HLo#5C=pk^O<3W7{(BMShK!__V#ux@4`eDQS5R7HsDWRKu#hp$X{bnbq7B+ zY22$gZG)u(#Xn!zVCXL!KwhIbF#j+$+-M57()*Ff!N?S+TOn}UlBy3fhTt2{(S7byyk}5BD(NN5}!4{DT z5Fy&mu>p~;*^cRK4!C3@2!KAEWVHz85`qr3hcz8!z9f8KjG-#V*tdJ^MG_iZq-A`8 z6eZsE9I5mX0J+SeayH~!$&T(GmBMXU!l-x=?19sz>i}rOSSbaA>0g3EY}TVcL4LN4 z*c*sF7`}e;e2gKyE$21Q;0F)y*slF6qWf#K*q@asW@>Imb9c5Ol1Q2PU1hHy-yN6# z6(rgO0z#Fa2!W=3sS2d>7X*d1A>tYyOtiGc`p$pA_3@8sOLd^I4`!~|n8o)Do~|R! z%mhI6kMA2P?`Ut8(~Al_E2r~|y79YdCz-`~2112=?n8Ck)jS&q=7m_1BT7u0i|5Rv z!nuNaQ#~e?$aiDUeh1alRam@h@H2!sVaX6t!h7g$^Gfrba%w#?curwHZuL3R3W!55 zQT#D!F4FZ+zZlHc_I_$9YN1!~wCv4woUfx++?Z5I0%d7CO+x zFlkh+pZ)AqCi_ld1v;jYY+VJqkU1sff(4aO2po&U3Mcn8H?ff(tmQ7etHPe;szcZv z<8buKiKT$LG%io}?1OX65|x%03}X+EF>(q2eZ|?PiR*Gqwar!cd5)mE)r%q8uC%QM zD8qLsV&_F)(Wp3nxd_uDmT=%yn#P!-d6nWM-=m_N)GcpT*rK{Ej~Anx%msCu#&ELSHwr`|u*MR}|$zcVS zh%{y2zUMvoiI2sta!omWD4s)eTG$^PzWO&jLlT6KkNI%*uQb)^w%32rJ2^A9B@ zcrnWb2-KW*M*D-PRB@-OPCr+-{jlj0H;z_EfL9(iiMi#Ss*pa+|0UP*l(jBlZZpY8 z$=%97hn*Iw>s8RwClmWPoX_;f9g7hFI!WP>#mCYWNUeZr6WhM`V}NBxj3PM;A4aXDvG4WN*zH&yK!u*!>>wgn&csA3Jucfc1SKfwtalcUA2z$3iKaxS23# z?_iF^UD!$=(^?L}3|5(L(RU3qI-`4Z=K6N+_weXHJjr+zy2Qnf5xSc>Hjvn(?pxN* z>+Dud`g8!T^T(hhh!Z`&tO)j5C-XQ`7OFV5qVAL-zlT~9zh&=DXJ&LXxlFQ{mD*{E zkvfzwO>l5US|QQk*8(+Ki8!jmK?8*P^p-Xqp=S&(fy5;X4$*DRN@wbeBeXSk=HVLK zM>REr*f9cw!N~hkDwtta9<<9SQ$j0zcePYbTTn?!iCVt<%O)7LCUURM+Bgp6v|n>t za`NRx^vta2fS|(pDj+dYc`gIkfq{XU$C>Q&?UAkL3mjg}irxQ^#fZZgT5{4cx&D{u z+2(S5={Xn`>92{262Yox%wAolNM)FmpqwIptfQu;M#FUR)bgE9LtwUDB z=jmv6uFR$6%IJ4;SR1=To!e|d=J2+&Q7nr2?)+ah_cKb(+$Uf!vb>W_{ZSUNI(sHz9bix31!MBL4dMc2@m6c z!(4tjpe;tt$mGVjF;&Jf)!77x`4gSGFD+#>4%E9$y_nC7)t_u#qg$ct2IYBCBq{>^ zjg0Q=;I_172_Qj0~QpYNYJ_{dbyuQ2{ng3YYlRHxblw|dyqE5yJDNw%9>!ldTidwMlEF~?N=hYOC}&LlMfO&Of##~(8s(xf z{h3jnPQ9RjM7}jd4cJIh*a7w~eOA+t4SE9I23dwB?6s=`Xl+Gt;|e#KR! zll1lS0no?xEL{ghP-Bq8;Gkjjoo${ntU^a#;>csKZu)HBn(6a#7wpdOPgOG2ezt)$ zQGap;?7C3BZiRT5RD~iJebBz~YTIRIOJ00sI?JQs8_q6j@t|>l3iH{*Yw4=W&a}?E zl=j9b^^4qA5*pN`zD&ikK3_$TQ)$X)&R<@*y@Di61ZnDiHmwh=q<6TQbj&6~-VtvH zjz`prqS~`?@os$MH<4@o1?Ecsg7hUIVe|?O=ddtgR{Akhe{#Dr1 zaH{d&QW}-4k7hc#iZ8j8%{uoa&llpv_Tf~TJ=r8E1!HC=uB-giz{~ze@BLu!#-L2{gDbuB z20`s3s zv)QSWccLyCOpixq&gq`_IXBhJ$tNDz8KTU>2NO>2r%}y)rCN5tUWZ|(d#HtZY>zXu zlPohIbthQPMzpL-Ta_kRUbwM|9;fXiUpKrHM^eK8^pFQ_@Lgc1leP25n{ivH2WDD- zww2(wJ9G*apL(%sJCdoWEmUeDFHE~JOzYLhrYPO_tcGruIh2X1_&kK!23q$b8O^4t z-Z}Q@5ruKDiF@h(+vx6&S~{*^)tSneovunA(O^!Ri`vilZBgDmU+buGV0FrB6d?hY z>*VAaT)Qya+hcraeL7iLk%ijX;0zpdfu5%$8miNIy`R+&1yaTof>}dTKhF}-VQz?p zU8L(~I?VBzUNvEHmP7d%nJ{0OBVzj7@fB$6PRvK$MyRnS-kuXpO6qI3g)hw0gJ);w zs-gE@=xyQW5GYESXt5gUXe!$*t;HN|3X03_$D>3M5NSN$Xd=c}JB(%dUdh?ZdZM(b z9PfQs+B$6q-4!#+QSA(AtLG?LNMy+$BK6_zEwqqX)=rk6_!~ zBz)eUz=ccpF?R|kDrzX@3A5h9(sA_13m;ynD-@?-@5PB1T*bBF4nxc-x#)b zI9BoVX!gD}puzaiL8xA*jJyptZ#vJo|C0A6=jzyngCj$0dU_{-|zj7A? zu-+atGo2N8k7vN|Fxpn0C^dAHz_z8Tcyxp-7tY!;*>FhDy!DHsPdDp;oej%0Dy>RQ z5Y7_R=#6Y2_T&z=woJUe-9~SxRE|FQVKpak+E=5d#kkDJ_M>Bc&pga)XX8dZ9fsdb zrLn5Lycxq)5qlLJ*aoEZc&MQ%)ncfHC<|nBm2xZjAsReG1?hvv=}zwks-D^ z98uDksXNWWcXpO^@DS(Qy=Ol@)UwzLlwG@rOKnT#?|5Xp8T@fMH0#A@+n1gHZY31h zmi=xRhvfaUSm3=|g|KFNxC4lkcwsg2>w|C`Z^(g&m4=hI9QjPy4p&7sF)xrnlzUU3CIRR4fU1la`#vK2^zTjB+b7XexL z1RUME^MFkZ8b2@w6>!7l=dh@FB5(i{c@S%ecrDhW__NR4Rv^+x=z|Rg8!GRwk1Tfm z)9)X>{EzJgVC~E?!Czv%Rr#M3(JoPsl;6s@ow^EK?HGHDe3RDPJ_^ z-OAV1yq+yuaC`T!*S`2*X3m9(hMA&@Out@&*M7rE&+@Ye`CwGZcYvWQiFpfk*Nm04pJD9Le6Wmt2XA?<`YR^@i{ECiSmQ2X^$BA*yg9aew z8YnWN)xDdRun0!f!(rfNnJNe673K>U|N#o zf4iTaSXK$&dy1#=%-Uh8{7v!eF!Uw8el98PNIw?a?4%w`FV}jV;2_r#6kD$vEvKtO zsAo_&8eN6pDqtUji?^~I;2_ZJyD9W`XxvB%?z8txNgj5pSm3vkf zK=UGNt!L%%Cqo*Wx2(Xxn{M^30(fu5NOH91NmC=36Q*OVq zjHFBLf|t+eE6>TXjs?NNZ5M_{@r)xrQ!iJd+hUkH*R>R{VW|*0xfq>d4_f(?l2_b&lI3hK$y?3&;W3(u}WD!gU zZ9B&zC^Prz+WRirSp2dEZGr}M^7KZgx(4+7p^M9%av#wJ14TP+uKZpp$;`SOfNS z!}dgeV64K&FOnYfkA)H(%S{$5IhHMsifl~ULL5a5v-y0DNoiJm>CZdBEt8v z%(gWv&om})e$E8=;cEW4Fb@Gh|UZ3n(kBAun0azS{B%OGr@tLA!uF?wB z^~<7h5ILd8q8Fg}{Y!8Bu`NNw;z9u(mN^O=nq##VSoUyS-Hxsv+|8l(mP2i~xK;pp zwQj*r^$w?FCTpI%uX+Axt+rgfGFop%~fPU~l zLYLUlmYy^KnPFKN4I!KuBFSc>;?OZ1ia5R42K#TElU6{}Wfn;a2|2>4RU5`pR2Cg@ zdoc%+Dcy{6CsZKUUX~Vyx^oUjOvKO{ZRp|#=?iKQdR_*j!Mw17(q=YfK<0uKFG&j; zwzww+@-0RsG_zc`q{ z@4_rlyX6rO4=BhObyduQ&i+{MUV^zQiz)tL%txBF4viZxEunZM+9i8EF0S{>B?mmL zg!yP-lI+dV0ZN6xrT +(Nivesh Dommaraju, 2018/04/20) diff --git a/include/deal.II/fe/fe_enriched.h b/include/deal.II/fe/fe_enriched.h index 55e66f6720..67335a223e 100644 --- a/include/deal.II/fe/fe_enriched.h +++ b/include/deal.II/fe/fe_enriched.h @@ -24,9 +24,14 @@ #include #include +#include #include #include +#include +#include + +#include #include #include #include @@ -691,6 +696,315 @@ private: &output_data) const; }; + + +/** + * This namespace consists of a class needed to create a collection + * of FE_Enriched finite elements (hp::FECollection) + * to be used with hp::DoFHandler in a domain with multiple, possibly + * overlapping, sub-domains with individual enrichment functions. + * + * To create hp::FECollection, a graph coloring algorithm is used to assign + * colors to enrichment functions before creating hp::FECollection. Hence the + * name. + */ +namespace ColorEnriched +{ + /** + * An alias template for predicate function which returns a + * boolean for a Triangulation::cell_iterator object. + * + * This is used by helper functions and in the implementation of + * ColorEnriched::Helper class. + */ + template + using predicate_function = std::function::cell_iterator &)>; + + /** + * ColorEnriched::Helper class creates a collection of FE_Enriched finite + * elements (hp::FECollection) to be used with hp::DoFHandler in a domain + * with multiple, possibly overlapping, sub-domains with individual + * enrichment functions. Note that the overlapping regions may have + * multiple enrichment functions associated with them. This is implemented + * using a general constructor of FE_Enriched object which allows different + * enrichment functions. + * + * Consider a domain with multiple enriched sub-domains + * which are disjoint i.e. not connected with each other. + * To ensure $C^0$ continuity at the interface between + * the enriched sub-domain (characterized by a single enrichment + * function) and the non-enriched domain, we can use an FE_Enriched + * object in the enriched sub-domain and in the non-enriched domain + * a standard finite element (eg: FE_Q) wrapped into an FE_Enriched + * object (which internally uses a dominating FE_Nothing object). + * Refer to the documentation on FE_Enriched for more + * information on this. It is to be noted that an FE_Enriched + * object is constructed using a base FE + * (FiniteElement objects) and one or more + * enriched FEs. FE_Nothing is a dummy enriched FE. + * + * The situation becomes more + * complicated when two enriched sub-domains + * share an interface. When the number of enrichment functions are + * same for the sub-domains, FE_Enriched object of one sub-domain + * is constructed such that each enriched FE is paired (figuratively) with a + * FE_Nothing in the FE_Enriched object of the other sub-domain. + * For example, let the FEs fe_enr1 and fe_enr2, which will be + * used with enrichment functions, correspond + * to the two sub-domains. Then the FE_Enriched objects of the two + * sub-domains are built using + * [fe_base, fe_enr1, fe_nothing] and + * [fe_base, fe_nothing, fe_enr2] respectively. + * Note that the size of the vector of enriched FEs + * (used in FE_Enriched constructor) is equal to 2, the + * same as the number of enrichment functions. When the number of enrichment + * functions is not the same, additional enriched FEs are paired + * with FE_Nothing. This ensures that the enriched DOF's at the interface + * are set to zero by the DoFTools::make_hanging_node_constraints() function. + * Using these two strategies, we construct the appropriate FE_Enriched + * using the general constructor. Note that this is + * done on a mesh without hanging nodes. + * + * Now consider a domain with multiple sub-domains which may share + * an interface with each other. As discussed previously, + * the number of enriched FEs in the FE_Enriched object of each + * sub-domain needs to be equal to the number of sub-domains. This is because + * we are not using the information of how the domains are connected + * and any sub-domain may share interface with any other sub-domain (not + * considering overlaps for now!). However, in general, a given sub-domain + * shares an interface only with a few sub-domains. This warrants + * the use of a graph coloring algorithm to reduce + * the size of the vector of enriched FEs + * (used in the FE_Enriched constructor). By giving the sub-domains + * that share no interface the same color, a single 'std::function' + * that returns different enrichment functions for each + * sub-domain can be constructed. Then the size of the vector of enriched + * FEs is equal to the number of different colors + * used for predicates (or sub-domains). + * + * @note The graph coloring function, SparsityTools::color_sparsity_pattern, + * used for assigning colors to the sub-domains + * needs MPI (use Utilities::MPI::MPI_InitFinalize to initialize MPI + * and the necessary Zoltan setup). + * The coloring function, based on Zoltan, is a parallel coloring + * algorithm but is used in serial by SparsityTools::color_sparsity_pattern. + * + * Construction of the Helper class needs a base FiniteElement @p fe_base, + * an enriched FiniteElement @p fe_enriched (used for all the + * enrichment functions), a vector of predicate + * functions (used to define sub-domains) as well as the corresponding + * enrichment functions. The FECollection object, a collection of FE_Enriched + * objects to be used with an hp::DoFHandler object, can be retrieved + * using the member function build_fe_collection which also modifies the + * active FE indices of the hp::DoFHandler object (provided as an argument + * to the build_fe_collection function). + * + *

Simple example

+ * Consider a domain with three sub-domains defined by predicate functions. + * Different cells are associated with FE indices as shown in the following + * image. The three equal-sized square-shaped sub-domains 'a', 'b' + * and 'c' can be seen. The predicates associated with these sub-domains + * are also labeled 'a', 'b' and 'c'. + * The subdomains 'a' and 'b' intersect with cell labeled with FE + * index 3. The cells in 'c' are labeled with FE + * index 1. As can be seen, connections exist between 'a' and 'b', + * 'b' and 'c' but 'a' and 'c' are not connected. + * + * \htmlonly \endhtmlonly + * @image html 3source_fe_indices.png "Active FE indices" + * + * As discussed before, the colors of predicates are alloted using + * the graph coloring algorithm. Each predicate is a node in the graph and if + * two sub-domains share an interface, the corresponding predicates + * should be given different colors. + * (The predicate colors are different from what is shown + * in the image. The colors in the image are as per FE indices). + * Predicates 'a' and 'c' can be given the same color since they + * are not connected but the color given to 'b' has to be different from + * 'a' and 'c'. + * + * The name of finite element at an index (i) of @p fe_collection + * (hp::FECollection) can be obtained by + * fe_collection[index].get_name() and is + * show in the table below. Note that all the FE_Enriched elements + * are of the same size and FE_Nothing<2>(dominating) is used as + * discussed before. + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + * + *
Active FE indexFiniteElement name
0FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)]
1FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)]
2FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)]
3FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)]
+ * + * The internal data members used by this class need to be available when the + * problem is solved. This can be ensured by declaring the object static, + * which is deallocated only when the program terminates. An alternative + * would be to use it as a data member of the containing class. Since vector + * of predicates and enrichment functions may not be available while + * constructing the Helper, a 'std::shared_ptr' to Helper object can be used + * and constructed when the predicates and enrichment functions are + * available. + * + *

Example usage:

+ * @code + * FE_Q fe_base(2); + * FE_Q fe_enriched(1); + * std::vector< predicate_function > predicates; + * std::vector< std::shared_ptr> > enrichments; + * + * Triangulation triangulation; + * hp::DoFHandler dof_handler(triangulation); + * + * static ColorEnriched::Helper FE_helper(fe_base, + * fe_enriched, + * predicates, + * enrichments); + * const hp::FECollection& + * fe_collection(FE_helper.build_fe_collection(dof_handler)); + * @endcode + */ + template + struct Helper + { + /** + * Constructor for Helper class. + * + * @param fe_base A base FiniteElement + * @param fe_enriched An enriched FiniteElement + * @param predicates std::vector of predicates defining the sub-domains. + * @p predicates[i] returns true for a cell if it + * belongs to a sub-domain with index (i). + * @param enrichments std::vector of enrichment functions + */ + Helper(const FE_Q & fe_base, + const FE_Q & fe_enriched, + const std::vector> & predicates, + const std::vector>> &enrichments); + + /** + * Prepares an hp::DoFHandler object. The active FE indices of + * mesh cells are initialized to work with + * ColorEnriched::Helper::fe_collection. + * + * Returns an hp::FECollection object. + * + * @param dof_handler an hp::DoFHandler object + * @return hp::FECollection, a collection of + * finite elements needed by @p dof_handler. + */ + const hp::FECollection & + build_fe_collection(hp::DoFHandler &dof_handler); + + private: + /** + * Contains a collection of FiniteElement objects needed by an + * hp::DoFHandler object. + */ + hp::FECollection fe_collection; + + /** + * A base FiniteElement used for constructing FE_Enriched + * object required by ColorEnriched::Helper::fe_collection. + */ + const FE_Q &fe_base; + + /** + * An enriched FiniteElement used for constructing FE_Enriched + * object required by ColorEnriched::Helper::fe_collection. + */ + const FE_Q &fe_enriched; + + /** + * A finite element with zero degrees of freedom used for + * constructing FE_Enriched object required by + * ColorEnriched::Helper::fe_collection + */ + const FE_Nothing fe_nothing; + + /** + * std::vector of predicates defining the sub-domains. + * @p predicates[i] returns true for a cell if it + * belongs to a sub-domain with index (i). + */ + const std::vector> predicates; + + /** + * std::vector of enrichment functions corresponding + * to the predicates. These are needed while constructing + * ColorEnriched::Helper::fe_collection. + */ + const std::vector>> enrichments; + + /** + * An alias template for any callable target such as functions, lambda + * expressions, function objects that take a + * Triangulation::cell_iterator + * and return a pointer to Function. This is used to define + * Helper::color_enrichments + * which returns an enrichment function + * for a cell in Triangulation. + */ + using cell_iterator_function = std::function *( + const typename Triangulation::cell_iterator &)>; + + /** + * std::vector of functions that take in a cell + * and return a function pointer. These are needed while constructing + * fe_collection. + * + * color_enrichments[i](cell_iterator) calls the correct enrichment function + * (i.e. whose corresponding predicate has the color i) for the cell. + * Note that this call to cell_iterator returns the enrichment function + * which is a pointer to class derived from Function. + */ + std::vector color_enrichments; + + /** + * std::vector of colors (unsigned int) associated + * with each sub-domain. No two connected sub-domains (i.e. sub-domains that + * share a vertex) have the same color. + */ + std::vector predicate_colors; + + /** + * Total number of different colors in predicate_colors + */ + unsigned int num_colors; + + /** + * A map of maps used to associate + * a cell with a map that in turn associates colors of active predicates in + * the cell with corresponding predicate ids. + */ + std::map> + cellwise_color_predicate_map; + + /** + * A vector of different possible color sets for a given set of + * predicates and hp::DoFHandler object + */ + std::vector> fe_sets; + }; +} // namespace ColorEnriched //} DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/fe/fe_enriched.templates.h b/include/deal.II/fe/fe_enriched.templates.h new file mode 100644 index 0000000000..674aebef12 --- /dev/null +++ b/include/deal.II/fe/fe_enriched.templates.h @@ -0,0 +1,212 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +#ifndef dealii_fe_enriched_templates_h +#define dealii_fe_enriched_templates_h + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace ColorEnriched +{ + namespace internal + { + /** + * Returns true if there is a connection between subdomains in the mesh + * associated with @p hp::DoFHandler i.e., if the subdomains share at least + * a vertex. The two subdomains are defined by predicates provided by + * @p predicate_1 and @p predicate_2. A predicate is a function (or + * object of a type with an operator()) which takes in a cell iterator and + * gives a boolean. It is said to be active in a cell if it returns true. + * + * An example of a custom predicate is one that checks the distance from a + * fixed point. Note that the operator() takes in a cell iterator. Using the + * constructor, the fixed point and the distance can be chosen. + * @code + * + * struct predicate + * { + * predicate(const Point p, const int radius) + * :p(p),radius(radius){} + * + * template + * bool operator () (const Iterator &i) + * { + * return ( (i->center() - p).norm() < radius); + * } + * + * private: + * Point p; + * int radius; + * + * }; + * @endcode + * and then the function can be used as follows to find if the subdomains + * are connected. + * @code + * find_connection_between_subdomains + * (dof_handler, + * predicate(Point(0,0), 1) + * predicate(Point(2,2), 1)); + * @endcode + * + * @param[in] hp::DoFHandler object + * @param[in] predicate_1 A function (or object of a type with an + * operator()) defining the subdomain 1. The function takes in a cell and + * returns a boolean. + * @param[in] predicate_2 Same as @p predicate_1 but defines subdomain 2. + * @return A boolean "true" if the subdomains share atleast a vertex. + */ + template + bool + find_connection_between_subdomains( + const hp::DoFHandler & dof_handler, + const predicate_function &predicate_1, + const predicate_function &predicate_2); + + /** + * Assign colors to subdomains using Graph coloring algorithm where each + * subdomain is considered as a graph node. Subdomains which are + * connected i.e share atleast a vertex have different color. Each subdomain + * is defined using a predicate function of @p predicates. + * + * @param[in] dof_handler a hp::DoFHandler object + * @param[in] predicates predicates defining the subdomains + * @param[out] predicate_colors Colors (unsigned int) associated with each + * subdomain. + */ + template + unsigned int + color_predicates( + const hp::DoFHandler & dof_handler, + const std::vector> &predicates, + std::vector & predicate_colors); + + /** + * Used to construct data members @p cellwise_color_predicate_map and + * @p fe_sets of Helper class. Inputs are hp::DoFHandler object, + * vector of predicates and colors associated with them. Before calling + * this function, colors can be assigned to predicates (i.e subdomains) + * using the function color_predicates. + * + * Each active FE index has a set of colors associated with it. + * A cell with an active FE index i has a set of colors given by + * fe_sets[i]. An active FE index with color {a,b} + * means that the cell has two active predicates (i.e they return true + * for the cell) of color a and b. + * + * Eg: fe_sets = { {}, {1}, {2}, {1,2} } means + * Cells with active FE index 0 have no predicates associated. + * Cells with index 1 have a active predicate with color 1. + * Cells with index 2 have a active predicate with color 2. + * Cells with index 3 have active predicates with color 1 and color 2. + * + * A map of maps cellwise_color_predicate_map is used to associate + * predicate colors in cells with predicate ids. For this purpose, each + * cell is given a unique id which is stored in material id for now. + * When the grid is refined, material id is inherited to the children, so + * map which associates material id with color map will still be relevant. + * + * Now the color map can be explained with an example. If the cell with + * material id 100 has active predicates 4 (color = 1) and 5 (color = 2), + * the map will insert pairs (1, 4) and (2, 5) at key 100 (i.e unique id + * of cell is mapped with a map which associates color with predicate id). + * + * @param[in] dof_handler hp::DoFHandler object + * @param[in] predicates vector of predicates defining the subdomains. + * @p predicates[i] returns true for a cell if it + * belongs to subdomain with index i. + * @param[in] predicate_colors vector of colors (unsigned int) associated + * with each subdomain. + * @param[out] cellwise_color_predicate_map A map of maps used to associate + * predicate colors in cells with predicate ids. + * @param[out] fe_sets a vector of color lists + */ + template + void + set_cellwise_color_set_and_fe_index( + hp::DoFHandler & dof_handler, + const std::vector> &predicates, + const std::vector & predicate_colors, + std::map> + & cellwise_color_predicate_map, + std::vector> &fe_sets); + + /** + * A function that returns a vector of enrichment functions corresponding + * to a color. The size of the vector is equal to total number of different + * colors associated with predicates (i.e subdomains). + * + * Assume that a cell has a active predicates with ids 4 (color = 1) and + * 5 (color = 2). cellwise_color_predicate_map has this information + * provided we know the material id. + * + * Now a call to color_enrichment[1](cell) should in turn call + * enrichments[4](cell). + * + * @param[in] num_colors number of colors for predicates + * @param[in] enrichments vector of enrichment functions + * @param[in] cellwise_color_predicate_map A map of maps used to associate + * predicate colors in cells with predicate ids. + * @param[out] color_enrichments A vector of functions that take in cell + * and return a function pointer. + */ + template + void + make_colorwise_enrichment_functions( + const unsigned int & num_colors, + const std::vector>> &enrichments, + const std::map> + &cellwise_color_predicate_map, + std::vector *( + const typename Triangulation::cell_iterator &)>> + &color_enrichments); + + + /** + * Creates a hp::FECollection object constructed using FE_Enriched + * elements which itself is constructed using color enrichment functions + * and is of size equal to number of colors. + * + * @param[in] fe_sets a vector of color lists + * @param[in] color_enrichments A vector of functions that take in cell + * and return a function pointer. + * @param[in] fe_base base FiniteElement + * @param[in] fe_enriched enriched FiniteElements + * @param[in] fe_nothing a finite element with zero degrees of freedom + * @param[out] fe_collection a collection of + * finite elements + */ + template + void + make_fe_collection_from_colored_enrichments( + const unsigned int &num_colors, + const std::vector> + &fe_sets, // total list of color sets possible + const std::vector *( + const typename Triangulation::cell_iterator &)>> + &color_enrichments, // color wise enrichment functions + const FE_Q &fe_base, // basic FE element + const FE_Q + &fe_enriched, // FE element multiplied by enrichment function + const FE_Nothing &fe_nothing, + hp::FECollection &fe_collection); + } // namespace internal +} // namespace ColorEnriched + +DEAL_II_NAMESPACE_CLOSE + +#endif // dealii_fe_enriched_templates_h diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc index c08b19a20d..1798052138 100644 --- a/source/fe/fe_enriched.cc +++ b/source/fe/fe_enriched.cc @@ -17,8 +17,10 @@ #include #include +#include #include +#include DEAL_II_NAMESPACE_OPEN @@ -1048,6 +1050,374 @@ FE_Enriched::InternalData::get_fe_output_object( } +namespace ColorEnriched +{ + namespace internal + { + template + bool + find_connection_between_subdomains( + const hp::DoFHandler & dof_handler, + const predicate_function &predicate_1, + const predicate_function &predicate_2) + { + // Use a vector to mark vertices + std::vector vertices_subdomain_1( + dof_handler.get_triangulation().n_vertices(), false); + + // Mark vertices that belong to cells in subdomain 1 + for (auto cell : dof_handler.active_cell_iterators()) + if (predicate_1(cell)) // True ==> part of subdomain 1 + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + vertices_subdomain_1[cell->vertex_index(v)] = true; + + // Find if cells in subdomain 2 and subdomain 1 share vertices. + for (auto cell : dof_handler.active_cell_iterators()) + if (predicate_2(cell)) // True ==> part of subdomain 2 + for (unsigned int v = 0; v < GeometryInfo::vertices_per_cell; + ++v) + if (vertices_subdomain_1[cell->vertex_index(v)] == true) + { + return true; + } + return false; + } + + + + template + unsigned int + color_predicates( + const hp::DoFHandler & mesh, + const std::vector> &predicates, + std::vector & predicate_colors) + { + const unsigned int num_indices = predicates.size(); + + // Use sparsity pattern to represent connections between subdomains. + // Each predicate (i.e a subdomain) is a node in the graph. + DynamicSparsityPattern dsp; + dsp.reinit(num_indices, num_indices); + + /* + * Find connections between subdomains taken two at a time. + * If the connection exists, add it to a graph object represented + * by dynamic sparsity pattern. + */ + for (unsigned int i = 0; i < num_indices; ++i) + for (unsigned int j = i + 1; j < num_indices; ++j) + if (internal::find_connection_between_subdomains( + mesh, predicates[i], predicates[j])) + dsp.add(i, j); + + dsp.symmetrize(); + + // Copy dynamic sparsity pattern to sparsity pattern needed by + // coloring function + SparsityPattern sp_graph; + sp_graph.copy_from(dsp); + predicate_colors.resize(num_indices); + + // Assign each predicate with a color and return number of colors + return SparsityTools::color_sparsity_pattern(sp_graph, predicate_colors); + } + + + + template + void + set_cellwise_color_set_and_fe_index( + hp::DoFHandler & dof_handler, + const std::vector> &predicates, + const std::vector & predicate_colors, + std::map> + & cellwise_color_predicate_map, + std::vector> &fe_sets) + { + // clear output variables first + fe_sets.clear(); + cellwise_color_predicate_map.clear(); + + /* + * Add first element of fe_sets which is empty by default. This means that + * the default, FE index = 0 is associated with an empty set, since no + * predicate is active in these regions. + */ + fe_sets.resize(1); + + /* + * Loop through cells and find set of predicate colors associated + * with the cell. As an example, a cell with an FE index associated + * with colors {a,b} means that predicates active in the cell have + * colors a or b. + * + * Create new active FE index in case of the color + * set is not already listed in fe_sets. If the set already exists, + * find index of the set in fe_sets. In either case, use the id in + * fe_sets to modify cell->active_fe_index. + * + * Associate each cell_id with a set of pairs. The pair represents + * predicate color and the active predicate with that color. + * Each color can only correspond to a single predicate since + * predicates with the same color correspond to disjoint domains. + * This is what the graph coloring in color_predicates + * function ensures. The number of pairs is equal to the number + * of predicates active in the given cell. + */ + unsigned int map_index(0); + auto cell = dof_handler.begin_active(); + auto endc = dof_handler.end(); + for (; cell != endc; ++cell) + { + // set default FE index ==> no enrichment and no active predicates + cell->set_active_fe_index(0); + + // Give each cell a unique id, which the cellwise_color_predicate_map + // can later use to associate a colors of active predicates with + // the actual predicate id. + // + // When the grid is refined, material id is inherited to + // children cells. So, the cellwise_color_predicate_map stays + // relevant. + cell->set_material_id(map_index); + std::set color_list; + + // loop through active predicates for the cell and insert map. + // Eg: if the cell with material id 100 has active + // predicates 4 (color = 1) and 5 (color = 2), the map will insert + // pairs (1, 4) and (2, 5) at key 100 (i.e unique id of cell is + // mapped with a map which associates color with predicate id) + // Note that color list for the cell would be {1,2}. + for (unsigned int i = 0; i < predicates.size(); ++i) + { + if (predicates[i](cell)) + { + /* + * create a pair predicate color and predicate id and add it + * to a map associated with each enriched cell + */ + auto ret = cellwise_color_predicate_map[map_index].insert( + std::pair(predicate_colors[i], + i)); + + AssertThrow( + ret.second == 1, + ExcMessage("Only one enrichment function per color")); + + color_list.insert(predicate_colors[i]); + } + } + + + /* + * check if color combination is already added. + * If already added, set the active FE index based on + * its index in the fe_sets. If the combination doesn't + * exist, add the set to fe_sets and once again set the + * active FE index as last index in fe_sets. + * + * Eg: if cell has color list {1,2} associated and + * fe_sets = { {}, {1}, {2} } for now, we need to add a new set {1,2} + * to fe_sets and a new active FE index 3 because 0 to 2 FE indices + * are already taken by existing sets in fe_sets. + */ + bool found = false; + if (!color_list.empty()) + { + for (unsigned int j = 0; j < fe_sets.size(); ++j) + { + if (fe_sets[j] == color_list) + { + found = true; + cell->set_active_fe_index(j); + break; + } + } + + if (!found) + { + fe_sets.push_back(color_list); + cell->set_active_fe_index(fe_sets.size() - 1); + } + } + ++map_index; // map_index should be unique to cells + } + } + + + + template + void + make_colorwise_enrichment_functions( + const unsigned int & num_colors, + const std::vector>> &enrichments, + const std::map> + &cellwise_color_predicate_map, + std::vector *( + const typename Triangulation::cell_iterator &)>> + &color_enrichments) + { + color_enrichments.clear(); + + // Each color should be associated with a single enrichment function + // called color enrichment function which calls the correct enrichment + // function for a given cell. + // + // Assume that a cell has a active predicates with ids 4 (color = 1) and + // 5 (color = 2). cellwise_color_predicate_map has this information + // provided we know the material id. + // + // Now a call to color_enrichment[1](cell) should in turn call + // enrichments[4](cell). That is just what this lambda is doing. + color_enrichments.resize(num_colors); + for (unsigned int i = 0; i < num_colors; ++i) + { + color_enrichments[i] = + [&, i](const typename Triangulation::cell_iterator + &cell) { + unsigned int id = cell->material_id(); + + /* + * i'th color_enrichment function corresponds to (i+1)'th color. + * Since FE_Enriched takes function pointers, we return a + * function pointer. + */ + return enrichments[cellwise_color_predicate_map.at(id).at(i + 1)] + .get(); + }; + } + } + + + + template + void + make_fe_collection_from_colored_enrichments( + const unsigned int & num_colors, + const std::vector> &fe_sets, + const std::vector *( + const typename Triangulation::cell_iterator &)>> + & color_enrichments, + const FE_Q & fe_base, + const FE_Q & fe_enriched, + const FE_Nothing &fe_nothing, + hp::FECollection &fe_collection) + { + // define dummy function which is associated with FE_Nothing + using cell_function = std::function *( + const typename Triangulation::cell_iterator &)>; + cell_function dummy_function; + dummy_function = + [=](const typename Triangulation::cell_iterator &) + -> const Function * { + AssertThrow(false, + ExcMessage("Called enrichment function for FE_Nothing")); + return nullptr; + }; + + using EnrichmentFunctions_2DVector = + std::vector *( + const typename Triangulation::cell_iterator &)>>>; + + // loop through color sets and create FE_enriched element for each + // of them provided before calling this function, we have color + // enrichment function associated with each color. + for (unsigned int color_set_id = 0; color_set_id != fe_sets.size(); + ++color_set_id) + { + std::vector *> vec_fe_enriched( + num_colors, &fe_nothing); + EnrichmentFunctions_2DVector functions(num_colors, {dummy_function}); + + for (auto it = fe_sets[color_set_id].begin(); + it != fe_sets[color_set_id].end(); + ++it) + { + // Given a color id ( = *it), corresponding color enrichment + // function is at index id-1 because color_enrichments is + // is indexed from zero. + const unsigned int ind = *it - 1; + + AssertIndexRange(ind, vec_fe_enriched.size()); + AssertIndexRange(ind, functions.size()); + AssertIndexRange(ind, color_enrichments.size()); + + // Assume a active predicate colors {1,2} for a cell. + // We then need to create a vector of FE enriched elements + // with vec_fe_enriched[0] = vec_fe_enriched[1] = &fe_enriched + // which can later be associated with enrichment functions. + vec_fe_enriched[ind] = &fe_enriched; + + // color_set_id'th color function is (color_set_id-1) + // element of color wise enrichments + functions[ind].assign(1, color_enrichments[ind]); + } + + AssertDimension(vec_fe_enriched.size(), functions.size()); + + FE_Enriched fe_component( + &fe_base, vec_fe_enriched, functions); + fe_collection.push_back(fe_component); + } + } + } // namespace internal + + + + template + Helper::Helper( + const FE_Q & fe_base, + const FE_Q & fe_enriched, + const std::vector> & predicates, + const std::vector>> &enrichments) : + fe_base(fe_base), + fe_enriched(fe_enriched), + fe_nothing(1, true), + predicates(predicates), + enrichments(enrichments) + {} + + + + template + const hp::FECollection & + Helper::build_fe_collection( + hp::DoFHandler &dof_handler) + { + // color the predicates based on connections between corresponding + // subdomains + if (predicates.size() != 0) + num_colors = + internal::color_predicates(dof_handler, predicates, predicate_colors); + else + num_colors = 0; + + // create color maps and color list for each cell + internal::set_cellwise_color_set_and_fe_index(dof_handler, + predicates, + predicate_colors, + cellwise_color_predicate_map, + fe_sets); + // setup color wise enrichment functions + // i'th function corresponds to (i+1) color! + internal::make_colorwise_enrichment_functions( + num_colors, enrichments, cellwise_color_predicate_map, color_enrichments); + + // make FE_Collection + internal::make_fe_collection_from_colored_enrichments(num_colors, + fe_sets, + color_enrichments, + fe_base, + fe_enriched, + fe_nothing, + fe_collection); + + return fe_collection; + } +} // namespace ColorEnriched + + // explicit instantiations #include "fe_enriched.inst" diff --git a/source/fe/fe_enriched.inst.in b/source/fe/fe_enriched.inst.in index f80860e9a0..41f5a67429 100644 --- a/source/fe/fe_enriched.inst.in +++ b/source/fe/fe_enriched.inst.in @@ -18,5 +18,60 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) { #if deal_II_dimension <= deal_II_space_dimension template class FE_Enriched; + + namespace ColorEnriched \{ + namespace internal \{ + template + bool find_connection_between_subdomains + (const hp::DoFHandler &dof_handler, + const predicate_function &predicate_1, + const predicate_function &predicate_2); + + + template + unsigned int color_predicates + (const hp::DoFHandler &dof_handler, + const std::vector> &, + std::vector &); + + + template + void set_cellwise_color_set_and_fe_index + (hp::DoFHandler &dof_handler, + const std::vector> + &predicates, + const std::vector &predicate_colors, + std::map > + &cellwise_color_predicate_map, + std::vector > &fe_sets); + + + template + void make_colorwise_enrichment_functions + (const unsigned int &num_colors, + const std::vector< std::shared_ptr > > + &enrichments, + const std::map > + &cellwise_color_predicate_map, + std::vector* + (const typename Triangulation::cell_iterator &)>> + &color_enrichments); + + + template + void make_fe_collection_from_colored_enrichments + (const unsigned int &num_colors, + const std::vector > &fe_sets, + const std::vector* + (const typename Triangulation::cell_iterator &)>> + &color_enrichments, + const FE_Q &fe_base, + const FE_Q &fe_enriched, + const FE_Nothing &fe_nothing, + hp::FECollection &fe_collection); + \} + template struct Helper; + \} + #endif } diff --git a/tests/fe/fe_enriched_color_01.cc b/tests/fe/fe_enriched_color_01.cc new file mode 100644 index 0000000000..ddee3960d3 --- /dev/null +++ b/tests/fe/fe_enriched_color_01.cc @@ -0,0 +1,148 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2001 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test ColorEnriched::internal::find_connection_between_subdomains(...) +// function. Check if the function correctly finds if two subdomains +// share an edge/node. + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +using namespace dealii; + +/* + * Construct a class template which finds if a cell + * is within region or not, based on distance of cell + * center from region center. + */ +template +struct predicate_template +{ + predicate_template(const Point p, const int radius) : + p(p), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) + { + return ((i->center() - p).norm() < radius); + } + +private: + Point p; + int radius; +}; + + + +template +void +test() +{ + deallog << "dim = " << dim << std::endl; + + // Construct grid + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + GridGenerator::hyper_cube(triangulation, -20, 20); + triangulation.refine_global(4); + + // Construct vector of predicates for 2 and 3 dimensions + Assert(dim == 2 || dim == 3, ExcDimensionMismatch2(dim, 2, 3)); + typedef std::function::active_cell_iterator &)> + predicate_function; + std::vector predicates; + predicates.resize(5); + if (dim == 2) + { + // Radius set such that every region has 2^2 = 4 cells + predicates[1] = predicate_template(Point(7.5, 7.5), 2); + predicates[2] = predicate_template(Point(5, 5), 2); + predicates[0] = predicate_template(Point(), 2); + predicates[3] = predicate_template(Point(-5, -5), 2); + predicates[4] = predicate_template(Point(-10, -10), 2); + } + else if (dim == 3) + { + // Radius set such that every region has 2^3 = 8 cells + predicates[1] = predicate_template(Point(7.5, 7.5, 7.5), 3); + predicates[2] = predicate_template(Point(5, 5, 5), 3); + predicates[0] = predicate_template(Point(), 3); + predicates[3] = predicate_template(Point(-5, -5, -5), 3); + predicates[4] = predicate_template(Point(-10, -10, -10), 3); + } + + // Check pair-wise connections between predicate regions + for (int i = 0; i < 5; ++i) + for (int j = 0; j < 5; ++j) + { + deallog << i << ":" << j << "=" + << ColorEnriched::internal::find_connection_between_subdomains( + dof_handler, predicates[i], predicates[j]) + << std::endl; + } +} + + +int +main() +{ + initlog(); + + try + { + test<2>(); + test<3>(); + } + catch (std::exception &exc) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl + << exc.what() << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl + << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl + << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + }; + return 0; +} diff --git a/tests/fe/fe_enriched_color_01.output b/tests/fe/fe_enriched_color_01.output new file mode 100644 index 0000000000..919b93e47b --- /dev/null +++ b/tests/fe/fe_enriched_color_01.output @@ -0,0 +1,53 @@ + +DEAL::dim = 2 +DEAL::0:0=1 +DEAL::0:1=0 +DEAL::0:2=1 +DEAL::0:3=1 +DEAL::0:4=0 +DEAL::1:0=0 +DEAL::1:1=1 +DEAL::1:2=1 +DEAL::1:3=0 +DEAL::1:4=0 +DEAL::2:0=1 +DEAL::2:1=1 +DEAL::2:2=1 +DEAL::2:3=0 +DEAL::2:4=0 +DEAL::3:0=1 +DEAL::3:1=0 +DEAL::3:2=0 +DEAL::3:3=1 +DEAL::3:4=1 +DEAL::4:0=0 +DEAL::4:1=0 +DEAL::4:2=0 +DEAL::4:3=1 +DEAL::4:4=1 +DEAL::dim = 3 +DEAL::0:0=1 +DEAL::0:1=0 +DEAL::0:2=1 +DEAL::0:3=1 +DEAL::0:4=0 +DEAL::1:0=0 +DEAL::1:1=1 +DEAL::1:2=1 +DEAL::1:3=0 +DEAL::1:4=0 +DEAL::2:0=1 +DEAL::2:1=1 +DEAL::2:2=1 +DEAL::2:3=0 +DEAL::2:4=0 +DEAL::3:0=1 +DEAL::3:1=0 +DEAL::3:2=0 +DEAL::3:3=1 +DEAL::3:4=1 +DEAL::4:0=0 +DEAL::4:1=0 +DEAL::4:2=0 +DEAL::4:3=1 +DEAL::4:4=1 diff --git a/tests/fe/fe_enriched_color_02.cc b/tests/fe/fe_enriched_color_02.cc new file mode 100644 index 0000000000..f0657fd095 --- /dev/null +++ b/tests/fe/fe_enriched_color_02.cc @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Test the function ColorEnriched::internal::color_predicates. + * Check if the function correctly colors vector of predicates + * using graph coloring. + * + * Two predicates are said to be connected if cells belonging to + * different predicates touch each other. + */ + +#include + +#include +#include +#include + +#include + +#include "../tests.h" + +/* + * Predicate function needed by ColorEnriched::internal::color_predicates + * implemented using a struct. + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + GridGenerator::hyper_cube(triangulation, -20, 20); + triangulation.refine_global(4); + hp::DoFHandler dof_handler(triangulation); + + // check the coloring function on different set of predicates. + std::vector> vec_predicates; + std::vector predicate_colors; + { + // case 1: predicates are not connected + vec_predicates.clear(); + vec_predicates.push_back(EnrichmentPredicate(Point(-10, 10), 2)); + vec_predicates.push_back(EnrichmentPredicate(Point(0, 0), 2)); + + predicate_colors.resize(vec_predicates.size()); + + ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + deallog << "Case 1" << std::endl; + for (auto i : predicate_colors) + { + deallog << i << std::endl; + } + } + + { + // case 2: Two predicates that are connected. + vec_predicates.clear(); + vec_predicates.push_back(EnrichmentPredicate(Point(-10, 10), 2)); + vec_predicates.push_back( + EnrichmentPredicate(Point(-7.5, 7.5), 2)); + + predicate_colors.resize(vec_predicates.size()); + + ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + deallog << "Case 2" << std::endl; + for (auto i : predicate_colors) + { + deallog << i << std::endl; + } + } + + { + // case 3: connections between predicates is as follows - + // 0-1 (overlap connection), + // 3-4 (edge connection) + vec_predicates.clear(); + vec_predicates.push_back(EnrichmentPredicate(Point(-10, 10), 2)); + vec_predicates.push_back( + EnrichmentPredicate(Point(-7.5, 7.5), 2)); + vec_predicates.push_back(EnrichmentPredicate(Point(0, 0), 2)); + vec_predicates.push_back( + EnrichmentPredicate(Point(7.5, -7.5), 2)); + vec_predicates.push_back( + EnrichmentPredicate(Point(12.5, -12.5), 2)); + + predicate_colors.resize(vec_predicates.size()); + + ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + deallog << "Case 3" << std::endl; + for (auto i : predicate_colors) + { + deallog << i << std::endl; + } + } + return 0; +} diff --git a/tests/fe/fe_enriched_color_02.with_mpi=true.output b/tests/fe/fe_enriched_color_02.with_mpi=true.output new file mode 100644 index 0000000000..64cb19dc5a --- /dev/null +++ b/tests/fe/fe_enriched_color_02.with_mpi=true.output @@ -0,0 +1,13 @@ + +DEAL:0::Case 1 +DEAL:0::1 +DEAL:0::1 +DEAL:0::Case 2 +DEAL:0::1 +DEAL:0::2 +DEAL:0::Case 3 +DEAL:0::1 +DEAL:0::2 +DEAL:0::1 +DEAL:0::1 +DEAL:0::2 diff --git a/tests/fe/fe_enriched_color_03.cc b/tests/fe/fe_enriched_color_03.cc new file mode 100644 index 0000000000..24b69ad8b1 --- /dev/null +++ b/tests/fe/fe_enriched_color_03.cc @@ -0,0 +1,175 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Test function - ColorEnriched::internal::set_cellwise_color_set_and_fe_index + * for a set of predicates. + * Check for each cell, if appropriate fe index and color-index map is set. + * Color-index map associates different colors of different enrichment + * functions with corresponding enrichment function index. + */ + +#include +#include + +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +/* + * Predicate function needed by ColorEnriched::internal::color_predicates + * implemented using a struct. + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + GridGenerator::hyper_cube(triangulation, -2, 2); + triangulation.refine_global(2); + + // Make predicates. Predicate 0 and 1 overlap. + // Predicate 2 is connected to 0. + std::vector> vec_predicates; + vec_predicates.push_back(EnrichmentPredicate(Point(0, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(-1, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(1.5, -1.5), 1)); + + // Do manual coloring since we are not testing coloring function here! + std::vector predicate_colors; + predicate_colors.resize(vec_predicates.size()); + ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + // Make required objects to call function set_cellwise_color_set_and_fe_index + std::map> + cellwise_color_predicate_map; + std::vector> fe_sets; + ColorEnriched::internal::set_cellwise_color_set_and_fe_index( + dof_handler, + vec_predicates, + predicate_colors, + cellwise_color_predicate_map, + fe_sets); + + /* + * Run through active cells to check fe index, colors of + * enrichment functions associated with it. + * + * A unique color set corresponds to an fe index. + * + * Eg: If an fe index 1 corresponds to color set {2,3}, + * means that a cell with fe index 1 has enrichment functions + * which are colored 2 and 3. Here different enrichment function + * have same color 2 but for a given cell only one of them would + * be relevant. So all additional information we need is which + * enrichment function is relevant that has color 2. We need to + * do the same thing with color 2 as well. + * + * Each cell is assigned unique material id by the function + * set_cellwise_color_set_and_fe_index. Now using material id, + * each cell is associated with a map which assigns a color to a + * particular enrichment function id. + * + */ + auto cell = dof_handler.begin_active(); + auto endc = dof_handler.end(); + for (unsigned int cell_index = 0; cell != endc; ++cell, ++cell_index) + { + // print true predicates for a cell as a binary code + // 0 1 0 indicates predicate with index 2 is true in the cell. + deallog << cell->index() << ":predicates="; + for (auto predicate : vec_predicates) + deallog << predicate(cell) << ":"; + + // print color predicate pairs for a cell + // Note that here material id is used to identify cells + // Here (1,2) indicates predicate 2 of color 1 is relavant for cell. + deallog << "(color, enrichment_func_id):"; + for (auto color_predicate_pair : + cellwise_color_predicate_map[cell->material_id()]) + { + deallog << "(" << color_predicate_pair.first << "," + << color_predicate_pair.second << "):"; + } + + // For a cell, print fe active index and corresponding fe set. + //{1,2} indicates 2 enrichment functions of color 1 and 2 are relevant. + deallog << ":fe_active_index:" << cell->active_fe_index() << ":fe_set:"; + for (auto fe_set_element : fe_sets[cell->active_fe_index()]) + deallog << fe_set_element << ":"; + deallog << std::endl; + } + + return 0; +} diff --git a/tests/fe/fe_enriched_color_03.with_mpi=true.output b/tests/fe/fe_enriched_color_03.with_mpi=true.output new file mode 100644 index 0000000000..38128acdbd --- /dev/null +++ b/tests/fe/fe_enriched_color_03.with_mpi=true.output @@ -0,0 +1,17 @@ + +DEAL:0::0:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::1:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::2:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::3:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::4:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::5:predicates=0:0:1:(color, enrichment_func_id):(1,2)::fe_active_index:1:fe_set:1: +DEAL:0::6:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::7:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::8:predicates=0:1:0:(color, enrichment_func_id):(2,1)::fe_active_index:2:fe_set:2: +DEAL:0::9:predicates=1:1:0:(color, enrichment_func_id):(1,0):(2,1)::fe_active_index:3:fe_set:1:2: +DEAL:0::10:predicates=0:1:0:(color, enrichment_func_id):(2,1)::fe_active_index:2:fe_set:2: +DEAL:0::11:predicates=1:1:0:(color, enrichment_func_id):(1,0):(2,1)::fe_active_index:3:fe_set:1:2: +DEAL:0::12:predicates=1:0:0:(color, enrichment_func_id):(1,0)::fe_active_index:1:fe_set:1: +DEAL:0::13:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: +DEAL:0::14:predicates=1:0:0:(color, enrichment_func_id):(1,0)::fe_active_index:1:fe_set:1: +DEAL:0::15:predicates=0:0:0:(color, enrichment_func_id)::fe_active_index:0:fe_set: diff --git a/tests/fe/fe_enriched_color_04.cc b/tests/fe/fe_enriched_color_04.cc new file mode 100644 index 0000000000..ffd8c5a4db --- /dev/null +++ b/tests/fe/fe_enriched_color_04.cc @@ -0,0 +1,176 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Test function: ColorEnriched::internal::make_colorwise_enrichment_functions + * for a set of predicates. + * + * The container of enrichment functions created by the function is + * called color_enrichments because they return the correct enrichment + * function if cell and color of the function are provided. + */ + +#include + +#include +#include +#include + +#include + +#include + +#include "../tests.h" + +/* + * Predicate function needed by ColorEnriched::internal::color_predicates + * implemented using a struct. + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + GridGenerator::hyper_cube(triangulation, -2, 2); + triangulation.refine_global(2); + + + // Make predicates. Predicate 0 and 1 overlap. + // Predicate 2 is connected to 0. + std::vector> vec_predicates; + vec_predicates.push_back(EnrichmentPredicate(Point(0, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(-1, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(1.5, -1.5), 1)); + + // Do manual coloring since we are not testing coloring function here! + std::vector predicate_colors; + predicate_colors.resize(vec_predicates.size()); + unsigned int num_colors = ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + + // Make required objects to call function set_cellwise_color_set_and_fe_index + std::map> + cellwise_color_predicate_map; + std::vector> fe_sets; + ColorEnriched::internal::set_cellwise_color_set_and_fe_index( + dof_handler, + vec_predicates, + predicate_colors, + cellwise_color_predicate_map, + fe_sets); + + // Construct vector of enrichment functions + std::vector>> vec_enrichments; + vec_enrichments.reserve(vec_predicates.size()); + for (unsigned int i = 0; i < vec_predicates.size(); ++i) + { + // constant function is chosen as enrichment function + ConstantFunction func(i); + vec_enrichments.push_back(std::make_shared>(func)); + } + + // Construct container for color enrichment functions needed + // by function make_colorwise_enrichment_functions + std::vector *( + const typename Triangulation::cell_iterator &)>> + color_enrichments; + + ColorEnriched::internal::make_colorwise_enrichment_functions( + num_colors, + vec_enrichments, + cellwise_color_predicate_map, + color_enrichments); + + + deallog << "color wise enrichment functions:" << std::endl; + auto cell = dof_handler.begin_active(); + auto endc = dof_handler.end(); + for (unsigned int cell_index = 0; cell != endc; ++cell, ++cell_index) + { + // Print ids of predicates active in the cell + unsigned int cell_id = cell->index(); + deallog << cell_id << ":predicates="; + for (auto predicate : vec_predicates) + deallog << predicate(cell) << ":"; + + /* + * Check if a color and enrichment index map exists for the cell. + * If so print the color and corresponding enrichment functions + * value at the cell center. Note that indices of color enrichment + * functions starts with zero not one unlike color indices. + */ + if (cellwise_color_predicate_map.count(cell->material_id()) == 1) + for (unsigned int color = 1; color <= num_colors; ++color) + if (cellwise_color_predicate_map.at(cell->material_id()) + .count(color) == 1) + deallog << ":color:" << color << ":func_value:" + << color_enrichments[color - 1](cell)->value( + cell->center()); + + deallog << std::endl; + } + return 0; +} diff --git a/tests/fe/fe_enriched_color_04.with_mpi=true.output b/tests/fe/fe_enriched_color_04.with_mpi=true.output new file mode 100644 index 0000000000..96671dd233 --- /dev/null +++ b/tests/fe/fe_enriched_color_04.with_mpi=true.output @@ -0,0 +1,18 @@ + +DEAL:0::color wise enrichment functions: +DEAL:0::0:predicates=0:0:0: +DEAL:0::1:predicates=0:0:0: +DEAL:0::2:predicates=0:0:0: +DEAL:0::3:predicates=0:0:0: +DEAL:0::4:predicates=0:0:0: +DEAL:0::5:predicates=0:0:1::color:1:func_value:2.00000 +DEAL:0::6:predicates=0:0:0: +DEAL:0::7:predicates=0:0:0: +DEAL:0::8:predicates=0:1:0::color:2:func_value:1.00000 +DEAL:0::9:predicates=1:1:0::color:1:func_value:0.00000:color:2:func_value:1.00000 +DEAL:0::10:predicates=0:1:0::color:2:func_value:1.00000 +DEAL:0::11:predicates=1:1:0::color:1:func_value:0.00000:color:2:func_value:1.00000 +DEAL:0::12:predicates=1:0:0::color:1:func_value:0.00000 +DEAL:0::13:predicates=0:0:0: +DEAL:0::14:predicates=1:0:0::color:1:func_value:0.00000 +DEAL:0::15:predicates=0:0:0: diff --git a/tests/fe/fe_enriched_color_05.cc b/tests/fe/fe_enriched_color_05.cc new file mode 100644 index 0000000000..339e4a35fe --- /dev/null +++ b/tests/fe/fe_enriched_color_05.cc @@ -0,0 +1,264 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Test function: ColorEnriched::internal + * ::make_fe_collection_from_colored_enrichments for a set of predicates. + * + * The function return FE_Collection which is then printed to test. + */ + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +// uncomment when debugging +//#define DATA_OUT_FE_ENRICHED + +/* + * Predicate function needed by ColorEnriched::internal::color_predicates + * implemented using a struct. + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +template +void +plot_shape_function(hp::DoFHandler &dof_handler, unsigned int patches = 5) +{ + std::cout << "n_cells: " << dof_handler.get_triangulation().n_active_cells() + << std::endl; + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + // output to check if all is good: + std::vector> shape_functions; + std::vector names; + for (unsigned int s = 0; s < dof_handler.n_dofs(); s++) + { + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[s] = 1.0; + + // if the dof is constrained, first output unconstrained vector + if (constraints.is_constrained(s)) + { + names.push_back(std::string("UN_") + + dealii::Utilities::int_to_string(s, 2)); + shape_functions.push_back(shape_function); + } + + names.push_back(std::string("N_") + + dealii::Utilities::int_to_string(s, 2)); + + // make continuous/constrain: + constraints.distribute(shape_function); + shape_functions.push_back(shape_function); + } + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + + // get material ids: + Vector fe_index(dof_handler.get_triangulation().n_active_cells()); + typename hp::DoFHandler::active_cell_iterator cell = dof_handler + .begin_active(), + endc = dof_handler.end(); + for (unsigned int index = 0; cell != endc; ++cell, ++index) + { + fe_index[index] = cell->active_fe_index(); + } + data_out.add_data_vector(fe_index, "fe_index"); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector(shape_functions[i], names[i]); + + data_out.build_patches(patches); + + std::string filename = + "hp-shape_functions_" + dealii::Utilities::int_to_string(dim) + "D.vtu"; + std::ofstream output(filename.c_str()); + data_out.write_vtu(output); +} + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + GridGenerator::hyper_cube(triangulation, -2, 2); + triangulation.refine_global(2); + + // Make predicates. Predicate 0 and 1 overlap. + std::vector> vec_predicates; + vec_predicates.push_back(EnrichmentPredicate(Point(-1, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(0, 1), 1)); + + // find colors for predicates + std::vector predicate_colors; + predicate_colors.resize(vec_predicates.size()); + unsigned int num_colors = ColorEnriched::internal::color_predicates( + dof_handler, vec_predicates, predicate_colors); + + // Make required objects to call function set_cellwise_color_set_and_fe_index + std::map> + cellwise_color_predicate_map; + std::vector> fe_sets; + + ColorEnriched::internal::set_cellwise_color_set_and_fe_index( + dof_handler, + vec_predicates, + predicate_colors, + cellwise_color_predicate_map, + fe_sets); + + // Construct vector of enrichment functions + std::vector>> vec_enrichments; + vec_enrichments.reserve(vec_predicates.size()); + for (unsigned int i = 0; i < vec_predicates.size(); ++i) + { + // constant function. + ConstantFunction func(10 + i); // constant function + vec_enrichments.push_back(std::make_shared>(func)); + } + + // Construct container for color enrichment functions needed + // by function make_colorwise_enrichment_functions + std::vector *( + const typename Triangulation::cell_iterator &)>> + color_enrichments; + + ColorEnriched::internal::make_colorwise_enrichment_functions( + num_colors, // needs number of colors + vec_enrichments, // enrichment functions based on predicate id + cellwise_color_predicate_map, + color_enrichments); + + // Construct object needed to call make_fe_collection_from_colored_enrichments + FE_Q fe_base(2); + FE_Q fe_enriched(1); + FE_Nothing fe_nothing(1, true); + hp::FECollection fe_collection; + ColorEnriched::internal::make_fe_collection_from_colored_enrichments( + num_colors, + fe_sets, // total list of color sets possible + color_enrichments, // color wise enrichment functions + fe_base, // basic fe element + fe_enriched, // fe element multiplied by enrichment function + fe_nothing, + fe_collection); + + // print all the different fe sets needed by different cells + deallog << "fe sets:" << std::endl; + for (auto fe_set : fe_sets) + { + deallog << "color:"; + for (auto color : fe_set) + deallog << ":" << color; + deallog << std::endl; + } + + // check if fe_collection is correctly constructed by function + deallog << "fe_collection[index] mapping:" << std::endl; + for (unsigned int index = 0; index != fe_collection.size(); ++index) + { + deallog << "name:" << fe_collection[index].get_name() << std::endl; + deallog << "n_blocks:" << fe_collection[index].n_blocks() << std::endl; + deallog << "n_comp:" << fe_collection[index].n_components() << std::endl; + deallog << "n_dofs:" << fe_collection[index].n_dofs_per_cell() + << std::endl; + } + +#ifdef DATA_OUT_FE_ENRICHED + GridTools::partition_triangulation( + Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD), triangulation); + dof_handler.distribute_dofs(*fe_collection); + + plot_shape_function(dof_handler, 5); +#endif + + dof_handler.clear(); + return 0; +} diff --git a/tests/fe/fe_enriched_color_05.with_mpi=true.output b/tests/fe/fe_enriched_color_05.with_mpi=true.output new file mode 100644 index 0000000000..62582688e0 --- /dev/null +++ b/tests/fe/fe_enriched_color_05.with_mpi=true.output @@ -0,0 +1,23 @@ + +DEAL:0::fe sets: +DEAL:0::color: +DEAL:0::color::1 +DEAL:0::color::1:2 +DEAL:0::color::2 +DEAL:0::fe_collection[index] mapping: +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:9 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:13 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:17 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:13 diff --git a/tests/fe/fe_enriched_color_06.cc b/tests/fe/fe_enriched_color_06.cc new file mode 100644 index 0000000000..3b48a1b67b --- /dev/null +++ b/tests/fe/fe_enriched_color_06.cc @@ -0,0 +1,142 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Test function: ColorEnriched::internal + * ::make_fe_collection_from_colored_enrichments for a set of predicates. + * + * The function return FE_Collection which is then printed to test. + */ + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include +#include + +#include + +#include "../tests.h" + +/* + * Testing helper class + */ +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + + +/* + * Type used to defined vector of predicates needed by the function + * ColorEnriched::internal::color_predicates. + */ +template +using predicate_function = + std::function::cell_iterator &)>; + + + +int +main(int argc, char **argv) +{ + // Intialize MPI as required by Zoltan library used for graph coloring by this + // test. + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv); + MPILogInitAll all; + + // Make basic grid + const unsigned int dim = 2; + Triangulation triangulation; + hp::DoFHandler dof_handler(triangulation); + GridGenerator::hyper_cube(triangulation, -2, 2); + triangulation.refine_global(2); + + // Make predicates. Predicate 0 and 1 overlap. + std::vector> vec_predicates; + vec_predicates.push_back(EnrichmentPredicate(Point(-1, 1), 1)); + vec_predicates.push_back(EnrichmentPredicate(Point(0, 1), 1)); + + // Construct vector of enrichment functions + std::vector>> vec_enrichments; + vec_enrichments.reserve(vec_predicates.size()); + for (unsigned int i = 0; i < vec_predicates.size(); ++i) + { + // constant function. + ConstantFunction func(10 + i); // constant function + vec_enrichments.push_back(std::make_shared>(func)); + } + + // Construct helper class to construct fe collection + FE_Q fe_base(2); + FE_Q fe_enriched(1); + static ColorEnriched::Helper fe_space( + fe_base, fe_enriched, vec_predicates, vec_enrichments); + const hp::FECollection &fe_collection( + fe_space.build_fe_collection(dof_handler)); + + // check if fe_collection is correctly constructed by function + deallog << "fe_collection[index] mapping:" << std::endl; + for (unsigned int index = 0; index != fe_collection.size(); ++index) + { + deallog << "name:" << fe_collection[index].get_name() << std::endl; + deallog << "n_blocks:" << fe_collection[index].n_blocks() << std::endl; + deallog << "n_comp:" << fe_collection[index].n_components() << std::endl; + deallog << "n_dofs:" << fe_collection[index].n_dofs_per_cell() + << std::endl; + } + + dof_handler.clear(); + return 0; +} diff --git a/tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output b/tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output new file mode 100644 index 0000000000..bcfe1f9539 --- /dev/null +++ b/tests/fe/fe_enriched_color_06.with_gsl=true.with_mpi=true.output @@ -0,0 +1,18 @@ + +DEAL:0::fe_collection[index] mapping: +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Nothing<2>(dominating)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:9 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Nothing<2>(dominating)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:13 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Q<2>(1)-FE_Q<2>(1)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:17 +DEAL:0::name:FE_Enriched<2>[FE_Q<2>(2)-FE_Nothing<2>(dominating)-FE_Q<2>(1)] +DEAL:0::n_blocks:3 +DEAL:0::n_comp:1 +DEAL:0::n_dofs:13 diff --git a/tests/fe/fe_enriched_color_07.cc b/tests/fe/fe_enriched_color_07.cc new file mode 100644 index 0000000000..083c61ce30 --- /dev/null +++ b/tests/fe/fe_enriched_color_07.cc @@ -0,0 +1,1988 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 - 2017 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, 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. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +/* + * Testing the ColorEnriched::Helper class by solving problems + */ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include +#include + +#include "../tests.h" + +template +class SigmaFunction : public Function +{ + Point center; + FunctionParser func; + +public: + SigmaFunction() : Function(), func(1) + {} + + // to help with resize function. doesn't copy function parser(func)! + SigmaFunction(SigmaFunction &&other) : center(other.center), func(1) + {} + + void + initialize(const Point & center, + const double & sigma, + const std::string &func_expr); + double + value(const Point &p, const unsigned int component = 0) const; + Tensor<1, dim> + gradient(const Point &p, const unsigned int component = 0) const; + virtual void + value_list(const std::vector> &points, + std::vector & value_list) const; +}; + +template +void +SigmaFunction::initialize(const Point & _center, + const double & sigma, + const std::string &func_expr) +{ + center = _center; + std::string variables; + std::map constants = {{"sigma", sigma}, + {"pi", numbers::PI}}; + + AssertThrow(dim == 1 || dim == 2 || dim == 3, + ExcMessage("Dimension not implemented")); + switch (dim) + { + case 1: + variables = "x"; + break; + case 2: + variables = "x, y"; + break; + case 3: + variables = "x, y, z"; + break; + } + + func.initialize(variables, func_expr, constants); +} + +template +inline double +SigmaFunction::value(const Point & p, + const unsigned int component) const +{ + const Point d(p - center); + return func.value(d, component); +} + + +template +inline Tensor<1, dim> +SigmaFunction::gradient(const Point & p, + const unsigned int component) const +{ + const Point d(p - center); + return func.gradient(d, component); +} + +template +void +SigmaFunction::value_list(const std::vector> &points, + std::vector & value_list) const +{ + const unsigned int n_points = points.size(); + + AssertDimension(points.size(), value_list.size()); + + for (unsigned int p = 0; p < n_points; ++p) + value_list[p] = value(points[p]); +} + + + +template +struct EnrichmentPredicate +{ + EnrichmentPredicate(const Point origin, const double radius) : + origin(origin), + radius(radius) + {} + + template + bool + operator()(const Iterator &i) const + { + return ((i->center() - origin).norm_square() < radius * radius); + } + + const Point & + get_origin() + { + return origin; + } + const double & + get_radius() + { + return radius; + } + +private: + const Point origin; + const double radius; +}; + + +template +class SplineEnrichmentFunction : public Function +{ +public: + SplineEnrichmentFunction(const Point & origin, + const std::vector &interpolation_points_1d, + const std::vector &interpolation_values_1d) : + Function(1), + origin(origin), + interpolation_points(interpolation_points_1d), + interpolation_values(interpolation_values_1d), + cspline(interpolation_points, interpolation_values) + {} + + SplineEnrichmentFunction(SplineEnrichmentFunction &&other) : + Function(1), + origin(other.origin), + interpolation_points(other.interpolation_points), + interpolation_values(other.interpolation_values), + cspline(interpolation_points, interpolation_values) + {} + + SplineEnrichmentFunction(const SplineEnrichmentFunction &other) : + Function(1), + origin(other.origin), + interpolation_points(other.interpolation_points), + interpolation_values(other.interpolation_values), + cspline(interpolation_points, interpolation_values) + {} + + + + virtual double + value(const Point &point, const unsigned int component = 0) const + { + Tensor<1, dim> dist = point - origin; + const double r = dist.norm(); + return cspline.value(Point<1>(r), component); + } + + virtual Tensor<1, dim> + gradient(const Point &p, const unsigned int component = 0) const + { + Tensor<1, dim> dist = p - origin; + const double r = dist.norm(); + Assert(r > 0., ExcDivideByZero()); + dist /= r; + Assert(component == 0, ExcMessage("Not implemented")); + return cspline.gradient(Point<1>(r))[0] * dist; + } + +private: + /** + * origin + */ + const Point origin; + std::vector interpolation_points; + std::vector interpolation_values; + // enrichment function as CSpline based on radius + Functions::CSpline<1> cspline; +}; + + + +struct ParameterCollection +{ + ParameterCollection(const std::string &file_name); + + ParameterCollection(const int & dim, + const double & size, + const unsigned int & shape, + const unsigned int & global_refinement, + const unsigned int & cycles, + const unsigned int & fe_base_degree, + const unsigned int & fe_enriched_degree, + const unsigned int & max_iterations, + const double & tolerance, + const std::string & rhs_value_expr, + const std::string & boundary_value_expr, + const std::string & rhs_radial_problem, + const std::string & boundary_radial_problem, + const std::string & exact_soln_expr, + const unsigned int & patches, + const unsigned int & debug_level, + const unsigned int & n_enrichments, + const std::vector &points_enrichments, + const std::vector &radii_predicates, + const std::vector &sigmas); + + void + print(); + + void set_enrichment_point(Point<2> &p, const unsigned int i) + { + AssertDimension(dim, 2); + p(0) = points_enrichments[2 * i]; + p(1) = points_enrichments[2 * i + 1]; + } + void set_enrichment_point(Point<3> &p, const unsigned int i) + { + AssertDimension(dim, 3); + p(0) = points_enrichments[3 * i]; + p(1) = points_enrichments[3 * i + 1]; + p(2) = points_enrichments[3 * i + 2]; + } + + int dim; + double size; + unsigned int shape; // 0 = ball, 1 = cube + unsigned int global_refinement; + unsigned int cycles; + unsigned int fe_base_degree; + unsigned int fe_enriched_degree; + unsigned int max_iterations; + double tolerance; + + // parameters related to exact solution + std::string rhs_value_expr; + std::string boundary_value_expr; + + // value = true ==> estimate exact solution from radial problem + std::string rhs_radial_problem; + std::string boundary_radial_problem; + + std::string exact_soln_expr; + + unsigned int patches; + // debug level = 0(output nothing), + // 1 (print statements) + // 2 (output solution) + // 3 (+ output grid data as well) + // 9 (+ shape functions as well) + unsigned int debug_level; + unsigned int n_enrichments; + std::vector points_enrichments; + std::vector radii_predicates; + std::vector sigmas; +}; + + + +ParameterCollection::ParameterCollection(const std::string &file_name) +{ + // std::cout << "...reading parameters" << std::endl; + + ParameterHandler prm; + + // declare parameters + prm.enter_subsection("geometry"); + prm.declare_entry("dim", "2", Patterns::Integer()); + prm.declare_entry("size", "1", Patterns::Double(0)); + prm.declare_entry("shape", "1", Patterns::Integer(0)); + prm.declare_entry("Global refinement", "1", Patterns::Integer(1)); + prm.declare_entry("cycles", "0", Patterns::Integer(0)); + prm.leave_subsection(); + + + prm.enter_subsection("solver"); + prm.declare_entry("fe base degree", "1", Patterns::Integer(1)); + prm.declare_entry("fe enriched degree", "1", Patterns::Integer(1)); + prm.declare_entry("max iterations", "1000", Patterns::Integer(1)); + prm.declare_entry("tolerance", "1e-8", Patterns::Double(0)); + prm.leave_subsection(); + + + prm.enter_subsection("expressions"); + prm.declare_entry("rhs value", "0", Patterns::Anything()); + prm.declare_entry("boundary value", "0", Patterns::Anything()); + prm.declare_entry("rhs value radial problem", "0", Patterns::Anything()); + prm.declare_entry("boundary value radial problem", "0", Patterns::Anything()); + prm.declare_entry("exact solution expression", "", Patterns::Anything()); + prm.declare_entry("estimate exact solution", "false", Patterns::Bool()); + prm.leave_subsection(); + + + prm.enter_subsection("output"); + prm.declare_entry("patches", "1", Patterns::Integer(1)); + prm.declare_entry("debug level", "0", Patterns::Integer(0, 9)); + prm.leave_subsection(); + + // parse parameter file + prm.parse_input(file_name, "#end-of-dealii parser"); + + + // get parameters + prm.enter_subsection("geometry"); + dim = prm.get_integer("dim"); + size = prm.get_double("size"); + shape = prm.get_integer("shape"); + global_refinement = prm.get_integer("Global refinement"); + cycles = prm.get_integer("cycles"); + prm.leave_subsection(); + + prm.enter_subsection("solver"); + fe_base_degree = prm.get_integer("fe base degree"); + fe_enriched_degree = prm.get_integer("fe enriched degree"); + max_iterations = prm.get_integer("max iterations"); + tolerance = prm.get_double("tolerance"); + prm.leave_subsection(); + + prm.enter_subsection("expressions"); + rhs_value_expr = prm.get("rhs value"); + boundary_value_expr = prm.get("boundary value"); + rhs_radial_problem = prm.get("rhs value radial problem"); + boundary_radial_problem = prm.get("boundary value radial problem"); + exact_soln_expr = prm.get("exact solution expression"); + prm.leave_subsection(); + + prm.enter_subsection("output"); + patches = prm.get_integer("patches"); + debug_level = prm.get_integer("debug level"); + prm.leave_subsection(); + + + // manual parsing + // open parameter file + std::ifstream prm_file(file_name); + + // read lines until "#end-of-dealii parser" is reached + std::string line; + while (getline(prm_file, line)) + if (line == "#end-of-dealii parser") + break; + + AssertThrow( + line == "#end-of-dealii parser", + ExcMessage("line missing in parameter file = \'#end-of-dealii parser\' ")); + + // function to read next line not starting with # or empty + auto read_next_proper_line = [&](std::string &line) { + while (getline(prm_file, line)) + { + if (line.size() == 0 || line[0] == '#' || line[0] == ' ') + continue; + else + break; + } + }; + + std::stringstream s_stream; + + // read num of enrichement points + read_next_proper_line(line); + s_stream.str(line); + s_stream >> n_enrichments; + + // note vector of points + for (unsigned int i = 0; i != n_enrichments; ++i) + { + read_next_proper_line(line); + s_stream.clear(); + s_stream.str(line); + + points_enrichments.resize(dim * n_enrichments); + + if (dim == 2) + { + double x, y; + s_stream >> x >> y; + points_enrichments[2 * i] = x; + points_enrichments[2 * i + 1] = y; + } + else if (dim == 3) + { + double x, y, z; + s_stream >> x >> y >> z; + points_enrichments[3 * i] = x; + points_enrichments[3 * i + 1] = y; + points_enrichments[3 * i + 2] = z; + } + else + AssertThrow(false, ExcMessage("Dimension not implemented")); + } + + // note vector of radii for predicates + for (unsigned int i = 0; i != n_enrichments; ++i) + { + read_next_proper_line(line); + s_stream.clear(); + s_stream.str(line); + + double r; + s_stream >> r; + radii_predicates.push_back(r); + } + + // note vector of sigmas for rhs + for (unsigned int i = 0; i != n_enrichments; ++i) + { + read_next_proper_line(line); + s_stream.clear(); + s_stream.str(line); + + double r; + s_stream >> r; + sigmas.push_back(r); + } +} + + + +ParameterCollection::ParameterCollection( + const int & dim, + const double & size, + const unsigned int & shape, + const unsigned int & global_refinement, + const unsigned int & cycles, + const unsigned int & fe_base_degree, + const unsigned int & fe_enriched_degree, + const unsigned int & max_iterations, + const double & tolerance, + const std::string & rhs_value_expr, + const std::string & boundary_value_expr, + const std::string & rhs_radial_problem, + const std::string & boundary_radial_problem, + const std::string & exact_soln_expr, + const unsigned int & patches, + const unsigned int & debug_level, + const unsigned int & n_enrichments, + const std::vector &points_enrichments, + const std::vector &radii_predicates, + const std::vector &sigmas) : + dim(dim), + size(size), + shape(shape), + global_refinement(global_refinement), + cycles(cycles), + fe_base_degree(fe_base_degree), + fe_enriched_degree(fe_enriched_degree), + max_iterations(max_iterations), + tolerance(tolerance), + rhs_value_expr(rhs_value_expr), + boundary_value_expr(boundary_value_expr), + rhs_radial_problem(rhs_radial_problem), + boundary_radial_problem(boundary_radial_problem), + exact_soln_expr(exact_soln_expr), + patches(patches), + debug_level(debug_level), + n_enrichments(n_enrichments), + points_enrichments(points_enrichments), + radii_predicates(radii_predicates), + sigmas(sigmas) +{} + + + +void +ParameterCollection::print() +{ + std::cout << "Dim : " << dim << std::endl + << "Size : " << size << std::endl + << "Shape : " << shape << std::endl + << "Global refinement : " << global_refinement << std::endl + << "Cycles : " << cycles << std::endl + << "FE base degree : " << fe_base_degree << std::endl + << "FE enriched degree : " << fe_enriched_degree << std::endl + << "Max Iterations : " << max_iterations << std::endl + << "Tolerance : " << tolerance << std::endl + << "rhs - main problem : " << rhs_value_expr << std::endl + << "boundary value - main problem : " << boundary_value_expr + << std::endl + << "rhs of radial problem : " << rhs_radial_problem << std::endl + << "boundary value of radial problem : " << boundary_radial_problem + << std::endl + << "exact solution expr : " << exact_soln_expr << std::endl + << "estimate exact solution using radial problem : " + << "Patches used for output: " << patches << std::endl + << "Debug level: " << debug_level << std::endl + << "Number of enrichments: " << n_enrichments << std::endl; + + std::cout << "Enrichment points : " << std::endl; + for (unsigned int i = 0; i < points_enrichments.size(); i = i + dim) + { + for (int d = 0; d < dim; ++d) + std::cout << points_enrichments[i + d] << " "; + + std::cout << std::endl; + } + + std::cout << "Enrichment radii : " << std::endl; + for (auto r : radii_predicates) + std::cout << r << std::endl; + + std::cout << "Sigma values of different sources : " << std::endl; + for (auto r : sigmas) + std::cout << r << std::endl; +} + + + +/* + * EstimateEnrichmentFunction is used to estimate enrichment function by + * solveing a 1D poisson problem with right hand side and boundary + * expression provided as a string of single variable 'x', to be + * interpreted as distance from @par center. + * + * Eg: For a 2D poisson problem with right hand side expression R and boundary + * expression B given by functions R( x^2 + y^2) and B( x^2 + y^2 ), an + * equivalent radial problem can be solved using this class by setting + * @par rhs_expr = R( x^2) + * @par boundary_expr = B (x^2) + * Note that the original poisson problem is defined by right hand side + * and boundary expression dependent on square of distance (x^2 + y^2) + * from center. + */ +class EstimateEnrichmentFunction +{ +public: + EstimateEnrichmentFunction(const Point<1> & center, + const double & domain_size, + const double & sigma, + const std::string &rhs_expr, + const std::string &boundary_expr, + const double & refinement = 11); + EstimateEnrichmentFunction(const Point<1> & center, + const double & left_bound, + const double & right_bound, + const double & sigma, + const std::string &rhs_expr, + const std::string &boundary_expr, + const double & refinement = 11); + ~EstimateEnrichmentFunction(); + void + run(); + void + evaluate_at_x_values(std::vector &interpolation_points, + std::vector &interpolation_values); + double + value(const Point<1> &p, const unsigned int &component = 0); + +private: + void + make_grid(); + void + setup_system(); + void + assemble_system(); + void + solve(); + void + refine_grid(); + void + output_results() const; + Point<1> center; + double domain_size; + double left_bound, right_bound; + double sigma; + std::string rhs_expr; + std::string boundary_expr; + std::vector rhs_values; + +public: + unsigned int debug_level; + +private: + Triangulation<1> triangulation; + unsigned int refinement; + FE_Q<1> fe; + DoFHandler<1> dof_handler; + SparsityPattern sparsity_pattern; + SparseMatrix system_matrix; + Vector solution; + Vector system_rhs; +}; + +EstimateEnrichmentFunction::EstimateEnrichmentFunction( + const Point<1> & center, + const double & domain_size, + const double & sigma, + const std::string &rhs_expr, + const std::string &boundary_expr, + const double & refinement) : + center(center), + domain_size(domain_size), + sigma(sigma), + rhs_expr(rhs_expr), + boundary_expr(boundary_expr), + debug_level(0), + refinement(refinement), + fe(1), + dof_handler(triangulation) +{ + left_bound = center[0] - domain_size / 2; + right_bound = center[0] + domain_size / 2; +} + + +EstimateEnrichmentFunction::EstimateEnrichmentFunction( + const Point<1> & center, + const double & left_bound, + const double & right_bound, + const double & sigma, + const std::string &rhs_expr, + const std::string &boundary_expr, + const double & refinement) : + center(center), + left_bound(left_bound), + right_bound(right_bound), + sigma(sigma), + rhs_expr(rhs_expr), + boundary_expr(boundary_expr), + debug_level(0), + refinement(refinement), + fe(1), + dof_handler(triangulation) +{ + domain_size = right_bound - left_bound; +} + + +void +EstimateEnrichmentFunction::make_grid() +{ + GridGenerator::hyper_cube(triangulation, left_bound, right_bound); + triangulation.refine_global(refinement); +} + + +void +EstimateEnrichmentFunction::setup_system() +{ + dof_handler.distribute_dofs(fe); + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); + solution.reinit(dof_handler.n_dofs()); + system_rhs.reinit(dof_handler.n_dofs()); +} + + +void +EstimateEnrichmentFunction::assemble_system() +{ + QGauss<1> quadrature_formula(2); + SigmaFunction<1> rhs; + rhs.initialize(center, sigma, rhs_expr); + FEValues<1> fe_values(fe, + quadrature_formula, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const unsigned int n_q_points = quadrature_formula.size(); + FullMatrix cell_matrix(dofs_per_cell, dofs_per_cell); + Vector cell_rhs(dofs_per_cell); + std::vector local_dof_indices(dofs_per_cell); + for (const auto &cell : dof_handler.active_cell_iterators()) + { + fe_values.reinit(cell); + cell_matrix = 0; + cell_rhs = 0; + rhs_values.resize(n_q_points); + rhs.value_list(fe_values.get_quadrature_points(), rhs_values); + + for (unsigned int q_index = 0; q_index < n_q_points; ++q_index) + { + double radius = center.distance(fe_values.quadrature_point(q_index)); + + //-1/r (r*u_r) = f form converts to + // r(u_r, v_r) = (r*f,v) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + cell_matrix(i, j) += + (radius * fe_values.shape_grad(i, q_index) * + fe_values.shape_grad(j, q_index)) * + fe_values.JxW(q_index); + cell_rhs(i) += + radius * (fe_values.shape_value(i, q_index) * + rhs_values[q_index] * fe_values.JxW(q_index)); + } + } + cell->get_dof_indices(local_dof_indices); + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = 0; j < dofs_per_cell; ++j) + system_matrix.add( + local_dof_indices[i], local_dof_indices[j], cell_matrix(i, j)); + system_rhs(local_dof_indices[i]) += cell_rhs(i); + } + } + std::map boundary_values; + SigmaFunction<1> boundary_func; + boundary_func.initialize(center, sigma, boundary_expr); + + VectorTools::interpolate_boundary_values( + dof_handler, 0, boundary_func, boundary_values); + VectorTools::interpolate_boundary_values( + dof_handler, 1, boundary_func, boundary_values); + + MatrixTools::apply_boundary_values( + boundary_values, system_matrix, solution, system_rhs); +} + + +void +EstimateEnrichmentFunction::solve() +{ + SolverControl solver_control(50000, 1e-12, false, false); + SolverCG<> solver(solver_control); + solver.solve(system_matrix, solution, system_rhs, PreconditionIdentity()); +} + + +void +EstimateEnrichmentFunction::refine_grid() +{ + Vector estimated_error_per_cell(triangulation.n_active_cells()); + KellyErrorEstimator<1>::estimate(dof_handler, + QGauss<1 - 1>(3), + typename FunctionMap<1>::type(), + solution, + estimated_error_per_cell); + GridRefinement::refine_and_coarsen_fixed_number( + triangulation, estimated_error_per_cell, 0.2, 0.01); + triangulation.execute_coarsening_and_refinement(); +} + + +void +EstimateEnrichmentFunction::output_results() const +{ + DataOut<1> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(solution, "solution"); + data_out.build_patches(); + std::ofstream output("radial_solution.vtk"); + data_out.write_vtk(output); +} + + +void +EstimateEnrichmentFunction::run() +{ + if (debug_level >= 1) + std::cout << "Solving problem in 1.: " << 1 << " with center: " << center + << ", size: " << domain_size << ", sigma: " << sigma << std::endl; + + make_grid(); + + double old_value = 0, value = 1, relative_change = 1; + bool start = true; + do + { + if (!start) + { + refine_grid(); + ++refinement; + } + + if (debug_level >= 1) + std::cout << "Refinement level: " << refinement << std::endl; + + setup_system(); + assemble_system(); + solve(); + + value = VectorTools::point_value(dof_handler, solution, center); + if (!start) + { + relative_change = fabs((old_value - value) / old_value); + } + start = false; + old_value = value; + } + while (relative_change > 0.005); + + if (debug_level >= 1) + std::cout << "Radial solution at origin = " << value + << " after global refinement " << refinement << std::endl; + + if (debug_level >= 1) + output_results(); +} + + +void +EstimateEnrichmentFunction::evaluate_at_x_values( + std::vector &interpolation_points, + std::vector &interpolation_values) +{ + if (interpolation_values.size() != interpolation_points.size()) + interpolation_values.resize(interpolation_points.size()); + + // x varies from 0 to 2*sigma. + // factor 2 because once a cell is decided to be enriched based on its center, + // its quadrature points can cause x to be twice! + for (unsigned int i = 0; i != interpolation_values.size(); ++i) + { + double value = VectorTools::point_value( + dof_handler, solution, Point<1>(interpolation_points[i])); + interpolation_values[i] = value; + } +} + + +double +EstimateEnrichmentFunction::value(const Point<1> & p, + const unsigned int &component) +{ + return VectorTools::point_value(dof_handler, solution, p); +} + + +EstimateEnrichmentFunction::~EstimateEnrichmentFunction() +{ + triangulation.clear(); +} + + + +template +void +plot_shape_function(hp::DoFHandler &dof_handler, unsigned int patches = 5) +{ + std::cout << "...start plotting shape function" << std::endl; + std::cout << "Patches for output: " << patches << std::endl; + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints(dof_handler, constraints); + constraints.close(); + + // find set of dofs which belong to enriched cells + std::set enriched_cell_dofs; + for (auto cell : dof_handler.active_cell_iterators()) + if (cell->active_fe_index() != 0) + { + unsigned int dofs_per_cell = cell->get_fe().dofs_per_cell; + std::vector local_dof_indices(dofs_per_cell); + cell->get_dof_indices(local_dof_indices); + enriched_cell_dofs.insert(local_dof_indices.begin(), + local_dof_indices.end()); + } + + // output to check if all is good: + std::vector> shape_functions; + std::vector names; + for (auto dof : enriched_cell_dofs) + { + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[dof] = 1.0; + + // if the dof is constrained, first output unconstrained vector + names.push_back(std::string("C_") + + dealii::Utilities::int_to_string(dof, 2)); + shape_functions.push_back(shape_function); + + // names.push_back(std::string("UC_") + + // dealii::Utilities::int_to_string(s,2)); + + // // make continuous/constraint: + // constraints.distribute(shape_function); + // shape_functions.push_back(shape_function); + } + + if (dof_handler.n_dofs() < 100) + { + std::cout << "...start printing support points" << std::endl; + + std::map> support_points; + MappingQ1 mapping; + hp::MappingCollection hp_mapping; + for (unsigned int i = 0; i < dof_handler.get_fe_collection().size(); ++i) + hp_mapping.push_back(mapping); + DoFTools::map_dofs_to_support_points( + hp_mapping, dof_handler, support_points); + + const std::string base_filename = + "DOFs" + dealii::Utilities::int_to_string(dim) + "_p" + + dealii::Utilities::int_to_string(0); + + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics " + << std::endl + << "unset ytics" << std::endl + << "unset grid" << std::endl + << "unset border" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" + << std::endl; + GridOut grid_out; + grid_out.write_gnuplot(dof_handler.get_triangulation(), f); + f << "e" << std::endl; + + DoFTools::write_gnuplot_dof_support_point_info(f, support_points); + + f << "e" << std::endl; + + std::cout << "...finished printing support points" << std::endl; + } + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + + // get material ids: + Vector fe_index(dof_handler.get_triangulation().n_active_cells()); + for (auto cell : dof_handler.active_cell_iterators()) + { + fe_index[cell->active_cell_index()] = cell->active_fe_index(); + } + data_out.add_data_vector(fe_index, "fe_index"); + + for (unsigned int i = 0; i < shape_functions.size(); i++) + data_out.add_data_vector(shape_functions[i], names[i]); + + data_out.build_patches(patches); + + std::string filename = "shape_functions.vtu"; + std::ofstream output(filename.c_str()); + data_out.write_vtu(output); + + std::cout << "...finished plotting shape functions" << std::endl; +} + + + +template +using predicate_function = + std::function::cell_iterator &)>; + + + +/** + * Main class + */ +template +class LaplaceProblem +{ +public: + LaplaceProblem(); + LaplaceProblem(const ParameterCollection &prm); + virtual ~LaplaceProblem(); + void + run(); + +protected: + void + initialize(); + void + build_fe_space(); + virtual void + make_enrichment_functions(); + void + setup_system(); + +private: + void + build_tables(); + void + assemble_system(); + unsigned int + solve(); + void + refine_grid(); + void + output_results(const unsigned int cycle); + void + process_solution(); + +protected: + ParameterCollection prm; + unsigned int n_enriched_cells; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + + std::shared_ptr> fe_collection; + hp::QCollection q_collection; + + FE_Q fe_base; + FE_Q fe_enriched; + FE_Nothing fe_nothing; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + PETScWrappers::MPI::SparseMatrix system_matrix; + PETScWrappers::MPI::Vector solution; + Vector localized_solution; + PETScWrappers::MPI::Vector system_rhs; + + ConstraintMatrix constraints; + + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + + ConditionalOStream pcout; + + std::vector> vec_rhs; + + using cell_iterator_function = std::function *( + const typename hp::DoFHandler::active_cell_iterator &)>; + + std::vector>> vec_enrichments; + std::vector> vec_predicates; + std::vector color_enrichments; + + // output vectors. size triangulation.n_active_cells() + // change to Vector + std::vector> predicate_output; + Vector color_output; + Vector vec_fe_index; + Vector mat_id; +}; + + + +template +LaplaceProblem::LaplaceProblem() : + prm(), + n_enriched_cells(0), + dof_handler(triangulation), + fe_base(prm.fe_base_degree), + fe_enriched(prm.fe_enriched_degree), + fe_nothing(1, true), + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout(std::cout, (this_mpi_process == 0)) +{ + prm.print(); + + pcout << "...default parameters set" << std::endl; +} + + + +template +LaplaceProblem::LaplaceProblem(const ParameterCollection &_par) : + prm(_par), + n_enriched_cells(0), + dof_handler(triangulation), + fe_base(prm.fe_base_degree), + fe_enriched(prm.fe_enriched_degree), + fe_nothing(1, true), + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes(Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process(Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout(std::cout, (this_mpi_process == 0) && (prm.debug_level >= 1)) +{ + AssertThrow(prm.dim == dim, ExcMessage("parameter file dim != problem dim")); + prm.print(); + pcout << "...parameters set" << std::endl; +} + + +/* + * Construct basic grid, vector of predicate functions and + * right hand side function of the problem. + */ +template +void +LaplaceProblem::initialize() +{ + pcout << "...Start initializing" << std::endl; + + /* + * set up basic grid which is a hyper cube or hyper ball based on + * parameter file. Refine as per the global refinement value in the + * parameter file. + * + */ + if (prm.shape == 1) + GridGenerator::hyper_cube(triangulation, -prm.size / 2.0, prm.size / 2.0); + else if (prm.shape == 0) + { + Point center = Point(); + GridGenerator::hyper_ball(triangulation, center, prm.size / 2.0); + triangulation.set_all_manifold_ids_on_boundary(0); + static SphericalManifold spherical_manifold(center); + triangulation.set_manifold(0, spherical_manifold); + } + else + AssertThrow(false, ExcMessage("Shape not implemented.")); + triangulation.refine_global(prm.global_refinement); + + /* + * Ensure that num of radii, sigma and focal points of enrichment domains + * provided matches the number of predicates. Since focal points of + * enrichment domains are stored as vector of numbers, the dimension of + * points_enrichments is dim times the n_enrichments in prm file. + */ + Assert( + prm.points_enrichments.size() / dim == prm.n_enrichments && + prm.radii_predicates.size() == prm.n_enrichments && + prm.sigmas.size() == prm.n_enrichments, + ExcMessage( + "Number of enrichment points, predicate radii and sigmas should be equal")); + + /* + * Construct vector of predicate functions, where a function at index i + * returns true for a given cell if it belongs to enrichment domain i. + * The decision is based on cell center being at a distance within radius + * of the domain from focal point of the domain. + */ + for (unsigned int i = 0; i != prm.n_enrichments; ++i) + { + Point p; + prm.set_enrichment_point(p, i); + vec_predicates.push_back( + EnrichmentPredicate(p, prm.radii_predicates[i])); + } + + /* + * Construct a vector of right hand side functions where the actual right hand + * side of problem is evaluated through summation of contributions from each + * of the functions in the vector. Each function i at a given point is + * evaluated with respect to point's position {x_i, y_i, z_i} relative to + * focal point of enrichment domain i. The value is then a function of x_i, + * y_i, z_i and sigmas[i] given by the rhs_value_expr[i] in parameter file. + */ + vec_rhs.resize(prm.n_enrichments); + for (unsigned int i = 0; i != prm.n_enrichments; ++i) + { + Point p; + prm.set_enrichment_point(p, i); + vec_rhs[i].initialize(p, prm.sigmas[i], prm.rhs_value_expr); + } + + pcout << "...finish initializing" << std::endl; +} + + + +/* + * Create a enrichment function associated with each enrichment domain. + * Here the approximate solution is assumed to be only dependent on + * distance (r) from the focal point of enrichment domain. This is + * true for poisson problem (a linear PDE) provided the right hand side is a + * radial function and drops exponentially for points farther from focal + * point of enrichment domain. + */ +template +void +LaplaceProblem::make_enrichment_functions() +{ + pcout << "!!! Make enrichment function called" << std::endl; + + for (unsigned int i = 0; i < vec_predicates.size(); ++i) + { + /* + * Formulate a 1d/radial problem with center and size appropriate + * to cover the enrichment domain. The function determining + * right hand side and boundary value of the problem is provided in + * parameter file. The sigma governing this function is the same as + * sigma provided for the corresponding enrichment domain and predicate + * function for which the enrichment function is to be estimated. + * + * The center is 0 since the function is evaluated with respect to + * relative position from focal point anyway. + * + * For hexahedral cells, dimension can extend upto sqrt(3) < 2 times! + * So take a factor of 4 as size of the problem. This ensures that + * enrichment function can be evaluated at all points in the enrichment + * domain + */ + double center = 0; + double sigma = prm.sigmas[i]; + double size = prm.radii_predicates[i] * 4; + + /* + * Radial problem is solved only when enrichment domain has a positive + * radius and hence non-empty domain. + */ + if (prm.radii_predicates[i] != 0) + { + EstimateEnrichmentFunction radial_problem( + Point<1>(center), + size, + sigma, + prm.rhs_radial_problem, + prm.boundary_radial_problem); + radial_problem.debug_level = prm.debug_level; // print output + radial_problem.run(); + pcout << "solved problem with " + << "x and sigma : " << center << ", " << sigma << std::endl; + + // make points at which solution needs to interpolated + std::vector interpolation_points, interpolation_values; + double cut_point = 3 * sigma; + unsigned int n1 = 15, n2 = 15; + double radius = size / 2; + double right_bound = center + radius; + double h1 = cut_point / n1, h2 = (radius - cut_point) / n2; + for (double p = center; p < center + cut_point; p += h1) + interpolation_points.push_back(p); + for (double p = center + cut_point; p < right_bound; p += h2) + interpolation_points.push_back(p); + interpolation_points.push_back(right_bound); + + // add enrichment function only when predicate radius is non-zero + radial_problem.evaluate_at_x_values(interpolation_points, + interpolation_values); + + + // construct enrichment function and push + Point p; + prm.set_enrichment_point(p, i); + SplineEnrichmentFunction func( + p, interpolation_points, interpolation_values); + vec_enrichments.push_back( + std::make_shared>(func)); + } + else + { + pcout << "Dummy function added at " << i << std::endl; + ConstantFunction func(0); + vec_enrichments.push_back( + std::make_shared>(func)); + } + } +} + + + +/* + * Since each enrichment domain has different enrichment function + * associated with it and the cells common to different enrichment + * domains need to treated differently, we use helper function in + * ColorEnriched namespace to construct finite element space and necessary + * data structures required by it. The helper function is also used to + * set FE indices of dof handler cells. We also set the cell with a unique + * material id for now, which is used to map cells with a pairs of + * color and corresponding enrichment function. All this is internally + * used to figure out the correct set of enrichment functions to be used + * for a cell. + * + * The quadrature point collection, with size equal to FE collection + * is also constructed here. + */ +template +void +LaplaceProblem::build_fe_space() +{ + pcout << "...building fe space" << std::endl; + + make_enrichment_functions(); + static ColorEnriched::Helper fe_space( + fe_base, fe_enriched, vec_predicates, vec_enrichments); + fe_collection = std::make_shared>( + fe_space.build_fe_collection(dof_handler)); + pcout << "size of fe collection: " << fe_collection->size() << std::endl; + + if (prm.debug_level == 9) + { + if (triangulation.n_active_cells() < 100) + { + pcout << "...start print fe indices" << std::endl; + + // print fe index + const std::string base_filename = + "fe_indices" + dealii::Utilities::int_to_string(dim) + "_p" + + dealii::Utilities::int_to_string(0); + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(triangulation, f); + f << "e" << std::endl; + + for (auto it : dof_handler.active_cell_iterators()) + f << it->center() << " \"" << it->active_fe_index() << "\"\n"; + + f << std::flush << "e" << std::endl; + pcout << "...finished print fe indices" << std::endl; + } + + if (triangulation.n_active_cells() < 100) + { + pcout << "...start print cell indices" << std::endl; + + // print cell ids + const std::string base_filename = + "cell_id" + dealii::Utilities::int_to_string(dim) + "_p" + + dealii::Utilities::int_to_string(0); + const std::string filename = base_filename + ".gp"; + std::ofstream f(filename.c_str()); + + f << "set terminal png size 400,410 enhanced font \"Helvetica,8\"" + << std::endl + << "set output \"" << base_filename << ".png\"" << std::endl + << "set size square" << std::endl + << "set view equal xy" << std::endl + << "unset xtics" << std::endl + << "unset ytics" << std::endl + << "plot '-' using 1:2 with lines notitle, '-' with labels point pt 2 offset 1,1 notitle" + << std::endl; + GridOut().write_gnuplot(triangulation, f); + f << "e" << std::endl; + + for (auto it : dof_handler.active_cell_iterators()) + f << it->center() << " \"" << it->index() << "\"\n"; + + f << std::flush << "e" << std::endl; + + pcout << "...end print cell indices" << std::endl; + } + } + + // q collections the same size as different material identities + q_collection.push_back(QGauss(4)); + for (unsigned int i = 1; i < fe_collection->size(); ++i) + q_collection.push_back(QGauss(10)); + + pcout << "...building fe space" << std::endl; +} + + + +template +void +LaplaceProblem::setup_system() +{ + pcout << "...start setup system" << std::endl; + + GridTools::partition_triangulation(n_mpi_processes, triangulation); + + dof_handler.distribute_dofs(*fe_collection); + + DoFRenumbering::subdomain_wise(dof_handler); + std::vector locally_owned_dofs_per_proc = + DoFTools::locally_owned_dofs_per_subdomain(dof_handler); + locally_owned_dofs = locally_owned_dofs_per_proc[this_mpi_process]; + locally_relevant_dofs.clear(); + DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs); + + constraints.clear(); + constraints.reinit(locally_relevant_dofs); + DoFTools::make_hanging_node_constraints(dof_handler, constraints); + + SigmaFunction boundary_value_func; + Point p; + prm.set_enrichment_point(p, 0); + boundary_value_func.initialize(p, prm.sigmas[0], prm.boundary_value_expr); + + + VectorTools::interpolate_boundary_values( + dof_handler, 0, boundary_value_func, constraints); + constraints.close(); + + // Initialise the stiffness and mass matrices + DynamicSparsityPattern dsp(locally_relevant_dofs); + DoFTools::make_sparsity_pattern(dof_handler, dsp, constraints, false); + std::vector n_locally_owned_dofs(n_mpi_processes); + for (unsigned int i = 0; i < n_mpi_processes; ++i) + n_locally_owned_dofs[i] = locally_owned_dofs_per_proc[i].n_elements(); + + SparsityTools::distribute_sparsity_pattern( + dsp, n_locally_owned_dofs, mpi_communicator, locally_relevant_dofs); + + system_matrix.reinit( + locally_owned_dofs, locally_owned_dofs, dsp, mpi_communicator); + + solution.reinit(locally_owned_dofs, mpi_communicator); + system_rhs.reinit(locally_owned_dofs, mpi_communicator); + + pcout << "...finished setup system" << std::endl; +} + + +template +void +LaplaceProblem::assemble_system() +{ + pcout << "...assemble system" << std::endl; + + system_matrix = 0; + system_rhs = 0; + + FullMatrix cell_system_matrix; + Vector cell_rhs; + + std::vector local_dof_indices; + + std::vector rhs_value, tmp_rhs_value; + + hp::FEValues fe_values_hp(*fe_collection, + q_collection, + update_values | update_gradients | + update_quadrature_points | + update_JxW_values); + + + for (auto cell : dof_handler.active_cell_iterators()) + if (cell->subdomain_id() == this_mpi_process) + { + fe_values_hp.reinit(cell); + const FEValues &fe_values = fe_values_hp.get_present_fe_values(); + + const unsigned int &dofs_per_cell = cell->get_fe().dofs_per_cell; + const unsigned int &n_q_points = fe_values.n_quadrature_points; + + /* + * Initialize rhs values vector to zero. Add values calculated + * from each of different rhs functions (vec_rhs). + */ + rhs_value.assign(n_q_points, 0); + tmp_rhs_value.assign(n_q_points, 0); + for (unsigned int i = 0; i < vec_rhs.size(); ++i) + { + vec_rhs[i].value_list(fe_values.get_quadrature_points(), + tmp_rhs_value); + + // add tmp to the total one at quadrature points + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + rhs_value[q_point] += tmp_rhs_value[q_point]; + } + } + + local_dof_indices.resize(dofs_per_cell); + cell_system_matrix.reinit(dofs_per_cell, dofs_per_cell); + cell_rhs.reinit(dofs_per_cell); + + cell_system_matrix = 0; + cell_rhs = 0; + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + for (unsigned int j = i; j < dofs_per_cell; ++j) + cell_system_matrix(i, j) += + (fe_values.shape_grad(i, q_point) * + fe_values.shape_grad(j, q_point) * fe_values.JxW(q_point)); + + cell_rhs(i) += + (rhs_value[q_point] * fe_values.shape_value(i, q_point) * + fe_values.JxW(q_point)); + } + + // exploit symmetry + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i; j < dofs_per_cell; ++j) + cell_system_matrix(j, i) = cell_system_matrix(i, j); + + cell->get_dof_indices(local_dof_indices); + + constraints.distribute_local_to_global(cell_system_matrix, + cell_rhs, + local_dof_indices, + system_matrix, + system_rhs); + } + + system_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + + pcout << "...finished assemble system" << std::endl; +} + +template +unsigned int +LaplaceProblem::solve() +{ + pcout << "...solving" << std::endl; + SolverControl solver_control(prm.max_iterations, prm.tolerance, false, false); + PETScWrappers::SolverCG cg(solver_control, mpi_communicator); + + // choose preconditioner +#define amg +#ifdef amg + PETScWrappers::PreconditionBoomerAMG::AdditionalData additional_data; + additional_data.symmetric_operator = true; + PETScWrappers::PreconditionBoomerAMG preconditioner(system_matrix, + additional_data); +#else + PETScWrappers::PreconditionJacobi preconditioner(system_matrix); +#endif + + cg.solve(system_matrix, solution, system_rhs, preconditioner); + + Vector local_soln(solution); + + constraints.distribute(local_soln); + solution = local_soln; + pcout << "...finished solving" << std::endl; + return solver_control.last_step(); +} + + + +template +void +LaplaceProblem::refine_grid() +{ + const Vector localized_solution(solution); + Vector local_error_per_cell(triangulation.n_active_cells()); + + hp::QCollection q_collection_face; + for (unsigned int i = 0; i < q_collection.size(); ++i) + q_collection_face.push_back(QGauss(1)); + + KellyErrorEstimator::estimate(dof_handler, + q_collection_face, + typename FunctionMap::type(), + localized_solution, + local_error_per_cell, + ComponentMask(), + nullptr, + n_mpi_processes, + this_mpi_process); + const unsigned int n_local_cells = + GridTools::count_cells_with_subdomain_association(triangulation, + this_mpi_process); + PETScWrappers::MPI::Vector distributed_all_errors( + mpi_communicator, triangulation.n_active_cells(), n_local_cells); + for (unsigned int i = 0; i < local_error_per_cell.size(); ++i) + if (local_error_per_cell(i) != 0) + distributed_all_errors(i) = local_error_per_cell(i); + distributed_all_errors.compress(VectorOperation::insert); + const Vector localized_all_errors(distributed_all_errors); + GridRefinement::refine_and_coarsen_fixed_fraction( + triangulation, localized_all_errors, 0.85, 0); + triangulation.execute_coarsening_and_refinement(); + ++prm.global_refinement; +} + + + +template +void +LaplaceProblem::output_results(const unsigned int cycle) +{ + pcout << "...output results" << std::endl; + pcout << "Patches used: " << prm.patches << std::endl; + + Vector exact_soln_vector, error_vector; + SigmaFunction exact_solution; + + if (prm.exact_soln_expr != "") + { + // create exact solution vector + exact_soln_vector.reinit(dof_handler.n_dofs()); + exact_solution.initialize( + Point(), prm.sigmas[0], prm.exact_soln_expr); + VectorTools::project(dof_handler, + constraints, + q_collection, + exact_solution, + exact_soln_vector); + + // create error vector + error_vector.reinit(dof_handler.n_dofs()); + Vector full_solution(localized_solution); + error_vector += full_solution; + error_vector -= exact_soln_vector; + } + + Assert(cycle < 10, ExcNotImplemented()); + if (this_mpi_process == 0) + { + std::string filename = "solution-"; + filename += Utilities::to_string(cycle); + filename += ".vtk"; + std::ofstream output(filename.c_str()); + + DataOut> data_out; + data_out.attach_dof_handler(dof_handler); + data_out.add_data_vector(localized_solution, "solution"); + if (prm.exact_soln_expr != "") + { + data_out.add_data_vector(exact_soln_vector, "exact_solution"); + data_out.add_data_vector(error_vector, "error_vector"); + } + data_out.build_patches(prm.patches); + data_out.write_vtk(output); + output.close(); + } + pcout << "...finished output results" << std::endl; +} + + + +// use this only when exact solution is known +template +void +LaplaceProblem::process_solution() +{ + Vector difference_per_cell(triangulation.n_active_cells()); + double L2_error, H1_error; + + if (!prm.exact_soln_expr.empty()) + { + pcout << "...using exact solution for error calculation" << std::endl; + + SigmaFunction exact_solution; + exact_solution.initialize( + Point(), prm.sigmas[0], prm.exact_soln_expr); + + VectorTools::integrate_difference(dof_handler, + localized_solution, + exact_solution, + difference_per_cell, + q_collection, + VectorTools::L2_norm); + L2_error = VectorTools::compute_global_error( + triangulation, difference_per_cell, VectorTools::L2_norm); + + VectorTools::integrate_difference(dof_handler, + localized_solution, + exact_solution, + difference_per_cell, + q_collection, + VectorTools::H1_norm); + H1_error = VectorTools::compute_global_error( + triangulation, difference_per_cell, VectorTools::H1_norm); + } + + pcout << "refinement h_smallest Dofs L2_norm H1_norm" << std::endl; + pcout << prm.global_refinement << " " + << prm.size / std::pow(2.0, prm.global_refinement) << " " + << dof_handler.n_dofs() << " " << L2_error << " " << H1_error + << std::endl; +} + + + +template +void +LaplaceProblem::run() +{ + pcout << "...run problem" << std::endl; + double norm_soln_old(0), norm_rel_change_old(1); + + // Run making grids and building fe space only once. + initialize(); + build_fe_space(); + + + if (this_mpi_process == 0) + deallog << "Solving problem with number of sources: " << prm.n_enrichments + << std::endl; + + for (unsigned int cycle = 0; cycle <= prm.cycles; ++cycle) + { + pcout << "Cycle " << cycle << std::endl; + + setup_system(); + + pcout << "Number of active cells: " + << triangulation.n_active_cells() << std::endl + << "Number of degrees of freedom: " << dof_handler.n_dofs() + << std::endl; + + if (prm.debug_level == 9 && this_mpi_process == 0) + plot_shape_function(dof_handler); + + assemble_system(); + auto n_iterations = solve(); + pcout << "Number of iterations: " << n_iterations << std::endl; + localized_solution.reinit(dof_handler.n_dofs()); + localized_solution = solution; + double value = + VectorTools::point_value(dof_handler, localized_solution, Point()); + pcout << "Solution at origin: " << value << std::endl; + + + // calculate L2 norm of solution + if (this_mpi_process == 0) + { + pcout << "calculating L2 norm of soln" << std::endl; + double norm_soln_new(0), norm_rel_change_new(0); + Vector difference_per_cell(triangulation.n_active_cells()); + VectorTools::integrate_difference(dof_handler, + localized_solution, + ZeroFunction(), + difference_per_cell, + q_collection, + VectorTools::H1_norm); + norm_soln_new = VectorTools::compute_global_error( + triangulation, difference_per_cell, VectorTools::H1_norm); + // relative change can only be calculated for cycle > 0 + if (cycle > 0) + { + norm_rel_change_new = + std::abs((norm_soln_new - norm_soln_old) / norm_soln_old); + pcout << "relative change of solution norm " + << norm_rel_change_new << std::endl; + } + + // moniter relative change of norm in later stages + if (cycle > 1) + { + deallog << (norm_rel_change_new < norm_rel_change_old) + << std::endl; + } + + norm_soln_old = norm_soln_new; + + // first sample of relative change of norm comes only cycle = 1 + if (cycle > 0) + norm_rel_change_old = norm_rel_change_new; + + pcout << "End of L2 calculation" << std::endl; + } + + if (prm.debug_level >= 2 && this_mpi_process == 0) + output_results(cycle); + + // Donot refine if loop is at the end + if (cycle != prm.cycles) + refine_grid(); + + pcout << "...step run complete" << std::endl; + } + pcout << "...finished run problem" << std::endl; +} + + + +template +LaplaceProblem::~LaplaceProblem() +{ + dof_handler.clear(); +} + + + +int +main(int argc, char **argv) +{ + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1); + MPILogInitAll all; + { + /* + * Case 1: single source with known solution + */ + { + ParameterCollection prm( + 2, // dimension + 2, // domain size + 1, // cube shape + 3, // global refinement + 5, // num of cycles grid is refined and solved again + 1, // fe base degree + 1, // fe enriched degree + 50000, // max iterations + 1e-9, // tolerance + // rhs value + "-(exp(-(x*x + y*y)/(2*sigma*sigma))*(- 2*sigma*sigma + x*x + y*y))/(2*sigma*sigma*sigma*sigma*sigma*sigma*pi)", + // boundary value + "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))", + // rhs value for radial problem solved to find enrichment function + "-(exp(-(x*x)/(2*sigma*sigma))*(- 2*sigma*sigma + x*x))/(2*sigma*sigma*sigma*sigma*sigma*sigma*pi)", + // boundary value for radial problem + "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))", + // exact solution expression. If null nothing is done + "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))", + 1, // patches + 1, // debug level + 1, // num enrichments + // enrichment points interpreted 2 at a time if dimension is 2 + {0, 0}, + // radii defining different predicates + {0.4}, + // sigmas defining different predicates + {0.1}); + + + if (prm.dim == 2) + { + LaplaceProblem<2> problem(prm); + problem.run(); + } + else if (prm.dim == 3) + { + LaplaceProblem<3> problem(prm); + problem.run(); + } + else + AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3")); + } + + + + /* + * Case 2: 3 sources + */ + { + ParameterCollection prm( + 2, // dimension + 4, // domain size + 1, // cube shape + 3, // global refinement + 4, // num of cycles grid is refined and solved again + 1, // fe base degree + 1, // fe enriched degree + 50000, // max iterations + 1e-9, // tolerance + // rhs value + "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y)/(2*sigma*sigma))", + // boundary value + "0", + // rhs value for radial problem solved to find enrichment function + "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))", + // boundary value for radial problem + "0", + // exact solution expression. If null nothing is done + "", + 1, // patches + 1, // debug level + 3, // num enrichments + // enrichment points interpreted 2 at a time if dimension is 2 + {0.5, 0.5, 0, 0, -1, -1}, + // radii defining different predicates + {0.4, 0.4, 0.4}, + // sigmas defining different predicates + {0.1, 0.1, 0.1}); + + + if (prm.dim == 2) + { + LaplaceProblem<2> problem(prm); + problem.run(); + } + else if (prm.dim == 3) + { + LaplaceProblem<3> problem(prm); + problem.run(); + } + else + AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3")); + } + + + /* + * Case 3: five sources 3d + */ + { + ParameterCollection prm( + 3, // dimension + 8, // domain size + 1, // cube shape + 4, // global refinement + 2, // num of cycles grid is refined and solved again + 1, // fe base degree + 1, // fe enriched degree + 50000, // max iterations + 1e-9, // tolerance + // rhs value + "1.0/(2*pi*sigma*sigma)*exp(-(x*x + y*y + z*z)/(2*sigma*sigma))", + // boundary value + "0", + // rhs value for radial problem solved to find enrichment function + "1.0/(2*pi*sigma*sigma)*exp(-(x*x)/(2*sigma*sigma))", + // boundary value for radial problem + "0", + // exact solution expression. If null nothing is done + "", + 1, // patches + 1, // debug level + 5, // num enrichments + // enrichment points interpreted 3 at a time if dimension is 3 + {1.5, 1.5, 1.5, 1, 1, 1, 0, 0, 0, -1, -1, -1, 1, -1, -1}, + // radii defining different predicates + {0.45, 0.45, 0.45, 0.45, 0.45}, + // sigmas defining different predicates + {0.1, 0.1, 0.1, 0.1, 0.1}); + + + if (prm.dim == 2) + { + LaplaceProblem<2> problem(prm); + problem.run(); + } + else if (prm.dim == 3) + { + LaplaceProblem<3> problem(prm); + problem.run(); + } + else + AssertThrow(false, ExcMessage("Dimension incorect. dim can be 2 or 3")); + } + } +} diff --git a/tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output b/tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output new file mode 100644 index 0000000000..c5bd2f26c5 --- /dev/null +++ b/tests/fe/fe_enriched_color_07.with_gsl=true.with_mpi=true.mpirun=2.output @@ -0,0 +1,14 @@ + +DEAL:0::Solving problem with number of sources: 1 +DEAL:0::1 +DEAL:0::1 +DEAL:0::1 +DEAL:0::1 +DEAL:0::Solving problem with number of sources: 3 +DEAL:0::1 +DEAL:0::1 +DEAL:0::1 +DEAL:0::Solving problem with number of sources: 5 +DEAL:0::1 + + -- 2.39.5