From 3fc80271ca27f808368c276bff5e892526c469be Mon Sep 17 00:00:00 2001 From: Denis Davydov Date: Wed, 28 Sep 2016 11:52:07 +0200 Subject: [PATCH] add FE_Enriched finite element --- doc/doxygen/images/fe_enriched_1d.png | Bin 0 -> 6521 bytes .../images/fe_enriched_h-refinement.png | Bin 0 -> 29424 bytes doc/doxygen/stylesheet.css | 8 + doc/news/changes.h | 6 + include/deal.II/fe/fe_enriched.h | 613 ++++++++++++ include/deal.II/fe/fe_system.h | 5 +- source/fe/CMakeLists.txt | 2 + source/fe/fe_enriched.cc | 926 ++++++++++++++++++ source/fe/fe_enriched.inst.in | 23 + tests/fe/fe_enriched_01.cc | 127 +++ tests/fe/fe_enriched_01.with_cxx11=on.output | 6 + tests/fe/fe_enriched_02.cc | 163 +++ tests/fe/fe_enriched_02.with_cxx11=on.output | 20 + tests/fe/fe_enriched_03.cc | 157 +++ tests/fe/fe_enriched_03.with_cxx11=on.output | 99 ++ tests/fe/fe_enriched_04.cc | 244 +++++ tests/fe/fe_enriched_04.with_cxx11=on.output | 11 + tests/fe/fe_enriched_05.cc | 186 ++++ tests/fe/fe_enriched_05.with_cxx11=on.output | 20 + tests/fe/fe_enriched_06.cc | 223 +++++ tests/fe/fe_enriched_06.with_cxx11=on.output | 22 + tests/fe/fe_enriched_07.cc | 237 +++++ tests/fe/fe_enriched_07.with_cxx11=on.output | 10 + tests/fe/fe_enriched_08.cc | 219 +++++ tests/fe/fe_enriched_08.with_cxx11=on.output | 6 + tests/fe/fe_enriched_08a.cc | 133 +++ tests/fe/fe_enriched_08a.with_cxx11=on.output | 6 + tests/fe/fe_enriched_compare_to_fe_system.cc | 320 ++++++ ..._compare_to_fe_system.with_cxx11=on.output | 5 + .../fe/fe_enriched_compare_to_fe_system_2.cc | 403 ++++++++ ...ompare_to_fe_system_2.with_cxx11=on.output | 5 + tests/fe/fe_enriched_step-36.cc | 779 +++++++++++++++ ...etsc=on.with_slepc=on.with_cxx11=on.output | 30 + tests/fe/fe_enriched_step-36b.cc | 903 +++++++++++++++++ ...etsc=on.with_slepc=on.with_cxx11=on.output | 12 + 35 files changed, 5928 insertions(+), 1 deletion(-) create mode 100644 doc/doxygen/images/fe_enriched_1d.png create mode 100644 doc/doxygen/images/fe_enriched_h-refinement.png create mode 100644 include/deal.II/fe/fe_enriched.h create mode 100644 source/fe/fe_enriched.cc create mode 100644 source/fe/fe_enriched.inst.in create mode 100644 tests/fe/fe_enriched_01.cc create mode 100644 tests/fe/fe_enriched_01.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_02.cc create mode 100644 tests/fe/fe_enriched_02.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_03.cc create mode 100644 tests/fe/fe_enriched_03.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_04.cc create mode 100644 tests/fe/fe_enriched_04.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_05.cc create mode 100644 tests/fe/fe_enriched_05.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_06.cc create mode 100644 tests/fe/fe_enriched_06.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_07.cc create mode 100644 tests/fe/fe_enriched_07.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_08.cc create mode 100644 tests/fe/fe_enriched_08.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_08a.cc create mode 100644 tests/fe/fe_enriched_08a.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_compare_to_fe_system.cc create mode 100644 tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_compare_to_fe_system_2.cc create mode 100644 tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_step-36.cc create mode 100644 tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output create mode 100644 tests/fe/fe_enriched_step-36b.cc create mode 100644 tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output diff --git a/doc/doxygen/images/fe_enriched_1d.png b/doc/doxygen/images/fe_enriched_1d.png new file mode 100644 index 0000000000000000000000000000000000000000..6d1af3172d1077d02fc688fecb9cc9c3ad7de296 GIT binary patch literal 6521 zcmaKQ2{@GB+y9xFXUt5tku~d(5M@c(Yep)R&>||K6v>*1hL6q3k8| zNTjblSyD0*CB@jsnB_nH|G(?{zt?-c*Lyw7xwmtl=f2Om&vl;96<23_$yLg$005GX z4mc3I5yC#m zUy>pR%ZYB_7HE|2N*Q7dWrCh5`jqe?V8En>g)cRU;;*X%+2bh_w_><(xNs+6!?2L* zDwN9r9o<48x`h99Oajmj0K8^R7y!ut2o00nwUV`fkP1$UI3Vl?rBhSFs)nf{I!h>I zyRn}QjHgefvo~!Lg*G-e*4EaBAV_(6c}`9al}h#V^V_?3uc@i2va+%$6%r;gvMbS{;zhfBtzz8e^Fhz)vYy(DGS^z=;U}_2xjVF@-Z2*J85Oo&~v0#y(C!$ah z`<)K}%J_hvh#`bI96kqt)Z4!g>ZXB|3IK}Tjyr4*M&JE2M7VX|YW3xjj;R-tK^>DY z--v;7X>-uqZt_#`<6!;h(~E}=aScq~f23L^Kr3nU=J3B#;QDEhv;#xYRt8$#fOtg= zcGwJd6rf0}CI7z)g%mqer4F8%08ftYVaGAH?-An^X@DfdEQzwXi2v9A0>oo!-}IC~ z_m~9`@cA+}2I1XIs!)oWkM&9JGAF=Ps&GZ4LDoJkg%% z{Kx_@ta2Q`;o(MV0njWfGULC;T-@l4wVe*0lNICSqBXl>V+9GQ5{S-d75w4$0QQ`e z41Rua8~H4pv7U&2eg38mdH|Pfo4~vzaSIbtv{fKR*P(Q!5~d? zNse7Eg+1PdXcA<#@om2HUt1DYACzEE*muR-@-l&EaDoZph*1!CvH~Pf{IM%NqR@uq zrd-64jWp8f3VeJF2nkBi?DE?um!!ILED;?X)s(*tCZ4X`BZxqKxr06%JV43f1Yn}Ih*0nOVD98WbOI&z*kYk1 z3A+9JB{>o-nZF?wRYLKrP}>=mH+bu@(Hn@6ynO7RsGM!E{Q;7WZScC;0|X-aPv&oz zjmRbHK{co{v%^If85QNEa>s~^E(Kpfyh6SM&IWtfYV#`45^{L@iEbxq?irX^WWI@d z1eNDPBkC#Tk}Cm5lwH>ZUM41il6sxnBW}E%g!kL9FgU&GHLx>#lMS7-Z*$o~!GBb& zB`<(ZEZwK-=kkOYp7!z^yHJG$ow?eS9}XXXMN=!0^iR5+A0rxKoeot8Y-e2b6nqr9 zm4Ck6BTYohzx~k(5Y`X{|NUxrJ{G=?woSLl>oLlIDoXg{XnqdB)95BEzF^>E@X$qd z8!j%b#l(D-7@l5ncFq@gEU?ryB6?r{qHe@Ci@fG=_}!R`X8eD~1u&+({=-33`zvG;#6cH_(vE;AlQ%oUd z)-JwkGelZdjyusyH?vJ8iRg?EBb_!#GCEneC?P)>CPrIs5d>ipkmQLOe^s9wNO-B_ z)#ewtpH79vtxqaB9ofT<#WJzuTOwA{a;QkaL_Cx(gzcCHl_(09>$ z6Kxma$hRM|KA+IR%X5UXleA`@CDuW=Yq0XTFYgQbfx_toM&(55oOqc7O^m;cb%@v( z(sGrOEKtR0%blh(9OI`@x7KzQkvqW}sCeD!AEN2(1diC^Jg9y4oK%?{Es<}DbGShI z5U}oaPJS$m#aH|(#lt07<)t6$q)JL~e`4yK0U)JO)K4XC#?Fh~sWHhZP}JMeB`d0r&PT$oFW%C(M{M2RrM|x@ORMK~W9r<7%EoeS zlGHLbDx!1$_IBBJzA`Tk(G;WAmP%TWz5`{s+k&m9Fx$+SJG9JbIrx}60;}`zUg5~u zWR(n4F?bEuYW`gvAuoWJYAMHmyZNxf&kk?!w}@hF@P|<;@%vv|n7zq4551m=ZHiDH zx4ggzkZsdXzsYztF|F2f<>#*@N{?+z`d*qB|E$+ll|1gfUz^*WkQL~4=pn9wpUIoe z>a-MyKWiY%*h0(W-z20**{5eS@*h&~$gfIwBsFRmB@Mrxh@W&S(wAq*)2{L&-HxR) z9Sjsmty`{Z?xx@7M2`M8@1SSi6z_s5rv%(oc3dq=pxDn{FODlZ!my$q>oqh-KsNCX z<#u4>FV1nX2iFe!thSBo(O5)H^cc9_C#zz63f&mHs20*KGWiHZ*Izbqv!bM;_srs8 z9(pZdO5KPn0W7gz*-+~@9Bg+Rx-~8h{;|OBlG$+0TLSszRXy};mz{*~ss}rKXo5GW zZNDC!`6unS@9qcN8HcD(UhUbZAQ4tGdc@k|q{ad5L$u?f#XhStHPgcR-2gtnH~MOF z>R%E&ozXWDTJ!1tDmcs6|e zHMX+63RNt%9=m9AxmB+#iLCYmg%RTDCJ$S6x~`J@PN2*Ls9EugG7$6L&r9ivaZby# zgp+D-cUVZ4yBX4++6~}=xvwH8^hbrqSrqHMQ`8LqV`!Xa*ait7>nm|iYR^}xrjbcu zcvF_wMKrn4d$Gj%alnF$M5J85<~WKK?}Vln3M<4t&mGNh+Nk6cX%HxD^L7)#*@Zl* zP!rjk%pLa_0n?rFVf_XL<^_s1bzDumuz1Hm`>GpTxyvU%fBa*y#4oAk*LzI~;NUHS zo=cJqszEqXEwa#SE<-EdiFF%F5xyM$kxT9SRF2HSdgTR+{k1df1SX ztwYDwzxxz%xHl>$mii*DUc<$#dbm`Vnij@+^{sejYlEBP3hWz!UQ*C-z~9}lK`l46 z%qu^Loc8%ysdrtpFg1{$3SlmOEgpH6)A_=z!R{b>Ns);cnLZP$dDd@Rj=3>i3fcS;0vzNd2e27#p3k<1FY-ZTPj=GRt)>ngC!%Y{&KXA>*N)jR(FCbj!vnGRD z2t*qR<#%bamu8DzD%czv`aJ90Soccbc>I~YJ|yr`4}Nm|Y98uUpf&yZ^u+_V#{MHu zj^}N8FImv}?&DFOwvw=;^#k-x&*G5-Eo;hJJ3QFK@GdSze>bbwJWp{Q5v&zgjOI7v zGSS0(lKKk|()$!jrR-YLwgclT3N2u1h_RN#Js$dXVjC;c!(Td(yk&W)yWrxlz_TFw-T|8Sl3IZb zXAtwH?w$7vzwP2)!AQYHUh;&pm4ZX$=PfrJ%Z~A)UR3;=dpvfV*SK@}5q#|K zkj~_m6uyTFfnw{VzbR`VTV#ib(a{$6HJnoS9abHb6S|u@ryZiM0^xz}YI_ZJB zzJ7svX(AgH=4;RXP&r(_x>T|hglJ%0`oC)!{_9Sa^IRC_`G#ZCQqQryfT6CfRJ@%d z+8so7V9a{vU|jNSZ$3_SYS@Gy7=Dy_b4U`Yr5I8J(9-b(ub7b)dfxn-=`*3u3$ITM z2rr+N2-8YW8}5pvJA#c~;w8Eo!?H2ov#Vmnzb|Bw!lJrc8_r~3+PTbiZ>-2jE6wRL zN0zi(es@l5*?Pg%6Ag!cX7_RN-__q0Qu_9?c@aIKPPbK_6$c25@xfC%4FDg2)3i44 zTS_QZi4v07;ZCB1m%SYaZ862SW%Xt|?{INZSQg`KFw@+)d_ z-4lvE^jt~v;lw90zD^G4Y$mUyTpBSzvBnN5lEVoRJBgsQ>g-H#RPVaTYY>`MdH$AQ zO#Ua01ia2f{Xv;EUPrNRUhAv(F;$fZ%~3DrS7orOVXEbWAx^cO4vO`%=OJ|rpThG$ z1JVo99~MesS9@S6SKl1U8_$I-^Rip|CoBd;3pIM|b#1$`Vp%ybev}jypH0l}AM^uC z`@UGGeauUf!B$IvOUd!0cZwZd#cNIhrNPVQUjKTl94Zc34H6?~s*J+>6C=x4u(0j1 z_~8rV?PW5N8Nm2K=8F1Bt<%XvM!?2srb^{EI`u;=7A{;jym)yh>uXA3Y~&j&V3Yi< z%B}*v@swu>tSpOX)P~tn*-ww?V#QM`k%(un<*M#uF z*i*J7VA#rgvc~=s-aIUGx@Fj3y!R zIByS*HXgQnE6_3umD@-JSC)^}RvnrZ$ub3o-h=ISuJ zDR+)f*`Cl>CB2`wXM{d&aV!ghAwd#5^`>6gX?Q~LlbtHBS9TH(;1!wG&CSZ3R9{;W zj(X%omBXo7pR!avw3(nFZrq@;mTAZtnhcnb@z6(F1r~o^8|?PU0X@stse8&-X?bWP z$pX%gdhaIu1|qBSA=Pq0@|3KeQ<)9y&N($ObJV?A+d~x*7iiB*o*IzNRP<0nqT}~- z6L(uEtQYsGAlt_I3UnI^z%FjaU0(f+Y^J=&YD6zyUu0ghNo#$XDs0k|-?DEyga=`D z8%}^-oC6O-Q*s%1k^MkyrKd2o?HRHOG|%}8Zxx)|U51Bgf=uW5;cUXc_y@TE{I7*% z;Ap=470>^1e^q~P1KLzp`#;A_<=8L90$W*gowJ&}2(<b?|Tx8?&!>{*YwXISt2pBGj8k5F?rnbN1<)K+741> z$F}~d*>_DTqr{Y|vn=_1b(=xDhEHUgYQ|Gad``W}hiIeJ?{Y0nLcRMFlYa{GEvhmm zIyW&j57Be}ZXXjk%|Pk}Kg5^Xw=8*fbzAisKP$<(i{#>|w4Wx9Rg~u&%*Jn8HIE5v^bX)9M&oA*5ty()r90y7Hb9@C z7*x%yyFOD>itdTd1)f}koM9ryR)=N?3|*L+9id5iq|ON?X`m=TrsXYd6=+hRgHIba zIzv}Qw__mY{dw66sjzWrUZ0Gtw&H0OR?wgIoM8Js=`lBp17Mb-0?zjKj~9|*jT+si zuO`zzkRDTBfd(V}XHYx(E(gP8=n*~;Dps+7NLy1!x@ELwPs@=w$5OI=xINhrCs%+N zZQDfOoWa-$@1e>KA8_0!qb+lq$a<&2Id(0HbeVDvG^A;7fu5lf$i=7eI^q>r_PIZb zJUuN-lkAw_0XYhL?X`4PukiC~B#V4E=HBl##ODUZTGVU(VY(lnYj{sUrQ#LnuwQ9B zmF@4}z0JZZqbz3Zg?xtQYWj@u4gqFmJ$$=yrqZzsP2se3M&m8wMgb2uhSwWQ*32IE1L{KADu);Q3tP zy+Ew?%G$^liIaGr)<-DY z$P^JsX3XYQo1-0EWhm@&_kn#g*C|_o34hJT1XhEV?M85RwkZrX$10($`zxKEmKe#S z$csqJpVMtDGh0c}(7ETGm^oG-;Hf26q?&cE>LgBW;4te5V>5+4_es5>S6felDELxU zZe1g%Cr%W2R$p*>LDVAwRtx0rBbS|rNy zF7o=uZwpeGp4)Iolun8bIByed%L$yDWR1;QbVc}XK(U#@B|~IIqpVh1Jo@wK7*F*M zOIU@E8^zW{66yE7=R^F5`m}HN7A>8NTvcsK;^TcVE(?s9iUxM|oMUB_V^$S0pY2G**#0cSf z_t9Q>dIr5xE&b^KfM2GB`kxW*{{{U1Zz{V{l{k;-&hNV5WKr>x8AU5;VEN0REU32T z;=$pU!_QGt;^1)cks!DcO}S-^Nfd)sk82oPfRlLph7lCW&O%yA{lEGAuUNZCC6Fe) Xv##sDBH=7T{DI?6=N%>6edzxO+@tI1 literal 0 HcmV?d00001 diff --git a/doc/doxygen/images/fe_enriched_h-refinement.png b/doc/doxygen/images/fe_enriched_h-refinement.png new file mode 100644 index 0000000000000000000000000000000000000000..909cee39fe65cd0f657825c6c52fc5c7a07783b1 GIT binary patch literal 29424 zcmc$Fdpy(szj&q64M|10btQ>eXw6)bB*`sBxl9PTZCJ+KRwF%*$Ug-#atR=zxjO6dHl}r{LcCP{ycx|@!sS0e!ZUe+w1vyy`#)b41^AxJg{TO z4k5!^H}3A(A;9GSKJOFYk4P!geEE+PPxT)>wSc=k^>^@f-l6*#?&y5V(96Nq`L46W zW0X&;vnGF%gxQ^Y`mK%C7Zea9{q4W0os;7O-IV&8s$w3GC-D9%amNl{nBfiGdjZ3< zsc&DLsVsVFH#^AQf~Oqn^fF8sy{y)Gt&dC|nz&=~ew1^;u-|k&`@(en$G;ZW%o|cN zFF2Z*9|h)Kc#W|*s!vjQeTI2dUlj50kIYE0{(m*=|D61vXZ_b8|9#ef1M>fm`2VKz z-$8CO@!#{ZJ&VuMf1LcEIsV^_|5vmAd&H4o^nw2@y8j5_e}t(2w>17u&i|3h|0s?B z8O^dxcmBUr^glw>wm1Cui2pCm=GW}rI%?Lytv_ub$2tn+J#ABM0qc`=l9sr{xPSg( zeSJ+t*&VrDH0dBXHZ!F0Iy$NUEG+sMJ)5$E9X|W)7K4 z!LW0G30;Uq2y=PmZ<9VC98ILtm!9;-w2R9NKJ_&r-hT zQty15Ni``>oVMUXqTNMgW0FDkHi!rDiY*)EtLjSj*H_shN7m!nc5g8j#W}Vo{~~*1 zUd?=5^@zB0ny*ToU?{B$*9`uGwr__|o&0h3?X_)M14S8L+O-a>s%uaSw!VT$7IJu- zhPXt7_w3|Ez^m2ufwz2%;A3)zYCE&yl+>zsCvIymHa!t8&uS0ADg0C#jNR5TDL>In zo5gGo;9M^9Z48+TG#NY-7`>RlH(&mn0|1d~u&*?r;|kCo^0?xHE0qtN6^Qylwt^VE?JUM|+5)S7f83^O z3MvC+e@#7G9A~dMXVu$Ntt}@}zO8z#2;Z`6{-ORQh?YcE^KhebQx~_FXSW3ro_MzD zZ})zAt6P2a?DDbNE5z_jWdAlryW4t>vf{9g_mV>mL_-mDCy^pwW>4k{ zddzXEE`_hCQ*A*=+Wy5A73(iHas`Q25RrAT#I{^?E z)L6hIe;?*Ja=s-XnttgyI7pv{{KLx*h8&pA$K#LA|7rZfBv%JJv26UVZX<0r z3oU0?tqpuEx^DR-b-J4?=l;#Tt^(_+-Owpm|(5$-i=m| z0{5vR?w8RTkVn`lE*Fu7i6U%7jr(h2N<+%@%ULm4OuU>b$GZWRwQ3f+P5ZIW_#0h~ z4Ffbx8z|ELebCFQvFAIUOWBn%b9A_#<%|b zquY+9gh@QGo|D?nLuM-5yUa;xs8$xm#PRd8qe*^T$R24d6zXyi=}aTPN|o>E`!@o` zO5r_`oaR644?wY~yB)0IzUPsth`3i6i|C#&-=s%SsE3-;<^X^jl2qPGw%{Xsct%f8 zo*5rCXbjV?idUhJ$ZZK*VDxC$kY9aFUTh1lFLuy5o3Owfzz1{cjy`$$0UdF9z&I;5 z{{+Ak_JX!4$O5PXdF}vyl5IH|v>*I*?qhl>BJZT_mx{C|5s&S#7Lyo+Lq#0VBizZz zlWsStvlSK2%0L)vUPtl~a4PTC9?kR0Cg^2_s)(r9>8f>^^#1xKZDBOGBhY<(J8vjB zF`etk#*R=+)8@TJ`b~!rVL-|`_YYs>|E0z7*^40mX2MgbjN8kwg4)J8sZ=dWs}mmF zdY#XJGag^{5eW|%nP6mhNBcBDyHl`P;IYnymxwFIZ-g(htx_>)NrueC=KV`bYBqb{ z#>0FNZ>;psZ`M>C23=BX=S37sE8 zqqagNRr39gU(fSH9`Ajm)+F#^-j!2O@|w|4R;|MD0NA+VoP(qF(Aff9 zPdIGAq=71hHAr~xM8vWxf#8Cm93Jz!S_pN34(G#P7Qmhq;a+ zGRB!uG4AS&mAovx*?dJnV_Vv{QMzg%@mB?71#cDRIZN;N!sY>CLPr}TfI+U~yZJNd zMkqUK*-ni%^}0%6;5+@b;X{HekD#8aH3~yc(^oBOt~EWJ^h|^Zg5@V0xF0SM9ifN4 z+3}|o#~L4?7T3q(xC-|as&z zl>0#vI{vf*w4q*AUbDHh9U?eu=1^1dRr%fn(+1bwj_awv^+}ZAoh%1Aua5b~<{KIP zk?G?0&1LT6*M&gd#kTY4PPiQn7`}2gBLF)@4wzcfiHgJ8qo3uo^KKFTXbLT{Ragb} z{0+|1l3(R*0Lb*k(wY7rdF0tuW1*?NpmgvH%|Nw_yD3_&fAIOl%cfoVshh+XRHIuEuW z4K~g#eAMRv;->C0lioTvq(eWCf%UowPU_l{19&E1^gYifKcTwuY=yJ z?R{#mwTti3SqVDNzhRyDb{e+YTkx_+R-nZqrKdBBlN-P(w`A6!*bMB`TYp9U5TD<} zqJ}LHkw0s~Y?BP|&CAYAP)5%{+p;PZ2mOZ3m(2+WU z_e0#6k2A$KG!^Bx1-CwB{=uU>@WPKkc8j zdYbhp7Q$DID!}ZOS|y;VHU+EYtsjUpfB}QDZ!0bJRd11Mz$Sk@rsaG=2_kDLW6YA9 zd@<0e?K7PyKZIxSIn2G9-g(gdS7O3&Bly6+Pg?`3tsHv@dFrLZuX(OjD?^i?U$=zf z#WBbI-e@109X6=6D~hToU(PQLKcs#IBXJerN-_Nd2bLW?#W1P5%KBBb-*F2ZFON~Z zJzlJyU0Hq!9;#SsTMQkE)(VDp_qKC57E>Z@F=}Da=zTG1Qq6P7)t7uRhA-_3IlVul zSjnAZwDEDBfyY=TuWl;Uo||6U=t)8Zd_C9wOS!Zi6(jYKU%_$h+TmXILBnOC?aHrD zxcV5)1#aFb9BQ>aeT&1e`2EychyYIO|GS~$PJe|+X@p;K3^h{~$7WWC04oaZdb`br zE7fK?v)OHFIUaK!xIeEst(t}WF#g$)B5Jc2wKpZfloV@>x-OJcvRPi8ym#jJm`=m> z)=WoadGeP#Rj80o8x@_T-11(0X6|AS=Yno?Q^M7Kb;UfF-KiMa55B*dx-5X9$-Cia zb}_}1S#Eyv@lxK=@+?Qv-WLegulK5-lDQ^T?=rG#Mj!5iAs;UV#b9LLEtTU6Pi>WF zzVEk2*(qfl@x=w+-r8dhiQw9Q!9*V-=mgigvCJvQel_B5@{4*-*BAJ3=6F9`mX62R zJ*7>4@MLPZ?>EWU9eo>z;hhMazCUDaWh;x2_cjpWc0;r*`0lBblLTl5JKEbggozgz zO?IU`efMz(0xT*Du=lx(@z68fdfV$GN&!8MK(M>+-~B$Bz|Ro^$(+7S*8<`Pj6A!l zf4QqI&Q7yp&((+_IqW6YBsn+xNw3T}Rcq~{Ht7J?Et0d7z|~L;k1C#v_#5u`0Az3* zniPVNIA&bnslroTW^zbD_hlVv-Z2gYI;`12RS2k4keJ?;iU|G9fbJ=}Qa6+@F?#sx zhLdPZdo=sA%k6nLvu8#lDQD|Yi}1<;v^3v_)AdkP67MD?H-+c8RT^s9?IH(1$X=~{ zms*d)xq}nolgyo3d$LxWIyak6ltC6rvhNX&p=zXq%mO-y&E()uJ~eN5Tq&^Lj1Je! zSp5<0#JcG15BvT(GBx6AG1^ji#9e1Qo18ekakN^s-yZ9T^*Oeoca-pMG+0>Ta>a<_ z#47GXGnndpdGU;`Sd}1ifC{>DQUi#-$?=jqobDUyr*Fs_9=!uKZVAxm=;ll3G}ZRb znil<1Z^+-%f#r)@mTS?P#yw}Vz8!6H15nA|o{*}+>qAwBgZe4L1WuPxPs!#Zu6#AT zUnqdpYVvz3Nqs*4G&8my|2Zn^J8{opF>}c584$#Bc)462ci1xYzWYgNm(!Dea-W*f zW2X%jU4us%#t%!5S99Mi+Puf`qel6_(~$fHYaqACmxhtgv244&T86cKh zxdCB_GAl}M=BoRo1pQgpgLIv%eH|x=3Q?S`D3q4xk&C$>-`Bym3{XKK4)>zNWVDnD zvlA(Q4-99`Gz~RlAm`WA#LJ6QVk20OzQNE5{5*)h<0$jx=Wj>NMRE0`er>@I+=*MN z*ZA2o8p_xmJsRu`;*tNd*09b-&6+g|=R)sVENTx#$EaVgHgV{YJ%jLFM+ zfRE)t83xx4+ci>wJX#;%_KDtLTO1%9n_nJYTh4yL_d8ZHF0fqq9CLvFVx{_SwY(&L ztIs!;>AQsF4Vsy3dAiu`;=hOBcp^y4}(O%JbQoZFte~EIK z*9bku09W2G9O1@_J&VDTi&fuxcXw8bcU<R6;CR??uk=j7vW*Nb3@vD5Ip4VU3n%!Fp1 z5~3)^rv?Uf#9(Lg+2-3FFu_AZM#F*N#Ln>d>^^P>fCVk1g)b(S?U-IJ(e{{^V}&+^ zMqt|yg}6xhB_qBaHXx-NaR$lehZrOGIrqYlYpNNX)A#IaPIko~jE#2l@b{R7WrXZP zj3|D{0mgc#J%~KZ?6^fL#EQuNc$*4`zpKtH8W|4Kg-o=gR9{*rpFX~#@y(0Lcb4o7 z{yr3fp#0zm_KQIKKR>^Hv>QP^YT_={WLDhpyXd_FcrN*A=A*Y?uX7V=ycMUKrU7S( z8i?Qv7O@{O3nxm~GH+cfyZM3b9OHFfI=k}xm;Aw~gTi0Lb1*AK?Gg1iT=*@(RiVwj zkcHIiu6%9SU+3Iiy?zNOYcQ9mYQ>ySX0Ol9VxIG2nCe7FGZ%l}tQ1RrpANB+A`l*9 ztjAA>!Nm>u^pv0|PToGj{T%f}64nW`5rtlJn@kCxw{BW}V}WaL-+I#m+^hR2AR+1J z4!ihh;H~Yc7(!x7O%#fp5oGGt)^Sh+hpPbySM)M*>+2aW4XSaTuWOF6mGmM;0YCYP zO-JtW0x<;;g}VAAV%4`xl@;G}<1|e;e@d|YF@3Nf@fHCKm1Iq6Q`mB>71ue<+AMv3 zmeS@s+I4MrLR$7oEBQBhxi}~KNWaP)&OHjhL^EYwHPU94ZnQfa4R#25CS_S^@Z78tVhrO7Q3PTDRkzfXg{lf9o&kB!Yo!(cbHg!5a*z^%_;b@ZP5lw3)1#;^$rm_ep1N z-VUk`Lr~pVVJ*^`8+~#N+xa`BwS7a*ekZpuLr7oYVl#X&w%<>b*>Zcq*0(F8dVOcg zYrDurBS72K(})au2e=7ie|w(Xarz)%3)*cf4=z03t$kp~xHpCp@$>Gz7vFVq1{e>2 z2l3Gs{P@#1MY%{79oI=0uEC*Qi`XR08Ua7^i%H4T z<0N3S#7A=b;vSji#Qy{Wuey>qW?V&m5!5x`TE{J&$u&fb5M`wIMzCBdKX^2|(0mY! zdI9mdizJlz!zta?QI(!i;M8W%irko<=tgov$~;TT-*Nt=#%x{Gx)^3TLQW_Z&GpfwdqT?I7}3-R0{QGd?5AkCbf=sW&K5I0d3Um7C0RY$Fp2&E}58tAkBTzPCo}f<)HmDYxDr5P7+48z+CO zF?BmdwL5m96E{`}cQ>}X0Y}*HdoSeaxIr5nOHWOpTXCR-ev`Xx^5v09#a?soy1ACV z!ZzW9V(bU6V-`~Tv#ZBn{N7N#4$r&@ga6Ji{At(q z7KOU&Tz*X>)B4-em^XA^xzPAZlY$C%Hq;fQ>;3>bSf=%6D{f{xHG_r)|#GAz?<-k2h@Ka`>s;- z-7-@KqJO~=l-qko@P&c?stbfFFdN+7De$VAPp8oz=gzcA{!S?wBX03){%O;-Il%r` zvo0!x6R}s=_ah@_r~lsmlR7%{7%<0NU3xa5xg+o;!qX>a?DhV2D<@%a$u1+@mS}<3 zOKZD>w|kcQQy(da#Sf@cO0LQhn3ZRoFm_aY2WQ{1*;(>&mR~dC9mPbC2P2USOXI6+ z&Od6?8gX82AG0d!yCVO@ept2sbj>%!OX4t~G zBmS-#<3C9kxq6;wc2fx7Z|@ByRdF=@kB7L*G8yFu{Ozuw0wq{gg6yvuB5)w*9R$p~ zOX`|^1yp<(7;iw8iq8Ht^Sin?G}Sxx&2T-y+gqL#y%&FSV_an+j9-gjQ=)MKT_+4k zJ(bBQqZcXrrp#m%>I$*;T@}^IyF-jgMoH!DB*ClbI_?yVCa9%drP0J)a0)12UI(^a zr&>AER_jk4tK2C(db{JMV_MNf#Y5D{{OikX_u6wOFt-_VCkNh7*tvep(Me#V{e5GU zD7h)2(T(+Bnw^~ajTY@G=Kc>%I-t}B1Q2Pa?g??X zczoL?yf&_uA{+u?O<{NF@Ve_q4t_|`75fd7sO?z&I;-yPA86ATSm*3vt6;Gj?BPsP zCyAw+^`^afjL+(AQeVc~A=DDeBz^>BT%dPiJv0!2E%#>A&sI z)@^&W$C=sok3leKnj_2PwRS7`On7L;a>MS;ap1r+`vk3{3%gDX_dLGU0~1!&eQwgJ z99s6&NnZ1ki}**s&rkP(;6fq7Y-0m);`9S~9udP{7J6ncK@9hyT{d7b@aRbYS%&YO z1GH0LXxFX3U`?wKNo5E^-23`t-MbvCwJ^H0QxiR|T!XZGlKC0yh~ z1>f*|t@G6|`wL?vpJd_T41?49`C5Jy6nh{XZ3XGo&%r8ZG^N! z&i&f8%pf){y-y|q7t7*@+w+uMzV4MbP45>)97U(UENquw+5JtQ(_j@XUw*PA2&Nn~ z_KfAAe2qAla`B3b#oX7EwePp((|#=7sI*p&N!gEh8qxhi7t*jEInAVHG<-h7?o}gu zX=jiA`CNglNOVd1Ix4pE7*hO+9;l8JS?7EpKZDl#vOLz!et7xDShrts-P?C~&9Z^^ zc6cOCOI?*!jkBB8{nl(%K@J`!Hxlm|N`vfqyjBcAVVykC1q?XubB%Xh$25Ikl>>Jd z`1NU9UbtE)Kz0OW@z~xvf3&9gGSGD8KInmSIuwKe_uY;EC_?4Z6v`!4h!X2k%<1-^ zPxd^5je$IK!R%)GBYeS1B@%T5SZUN?2pU zeuRiI!>uOvrCIAOEE5tEA8iS)OCxClgm7SmgaL_u^iRU#N=fe(oJQ4I{WHi!1n9A; ziEPS=2Zk%**UlY|KxL;Zlwn5>5~z6m=e?jH0B&tB>p65_QPMyVh}Qcd-8A4$xC=W7 zUqAlv06&(Lu2882A7w>k9ivhtZ_;&LGOdSnC~bkcrjtMV6FmQux<@(_@?*LEFxyHX z^r73;S5*Vj-_`J$YHvm&)IXcG%;Qte;^O3Sg`+9gl={g6*~pggp;kqu{X%V_FzD)_ zl8-8rr^1?yypw-DI1~LfBlqks=)j55C*yss0~c%1g60W~!=I^p=Pqc0Oer~eD;F|- z*kz2)xYi$UQi~G$b56HSh}DU`F*lDOVBRZ&nugEfe*f(ZGxr%PAd+}Fe< z|7J7U{zQ4lvGK#hfU|7&XQFN-N6hKzo&f1%k_cLUmy+odk00~I?4eK6Vi+6Mg=*ZD z(Q0ntI1_8dKi+uur}OpD!mww;Ew39qYe@PI=?wk5z8*32&o*PO+Xkxtaw&PzMM%Ur z2a{W!?eDm;?kCS9cel^VJz3ju#J}>cOR$g!EeMWj0Anl5vtArhjkB!Y$V7UxH-NDT zqoN}c14V*g__bZ)4(N=Y7Q}$J%VpE(lRatkWC}hWsn;8$Yps0U@Ny+N(1KC^ux|;HMDj zmb7w2z?K(+?Xq4r&disM$RZ3LxmcHPF*!_r6%_=73O{1;ORABb(Js1c>p}vKTbRq$ zoUQ39oQ0?CGO0I(E0qs$5vI2q1?jo=;8iTDhxAe1z-e6bE%6e|s4AV>H3 z0}-13)A$I~j;9R?TG(%P$I#bE{%5FrsQCU0SlR){H|;6GZ~J_Hj%vq^xI%yDC%p08_P?_?GYX+WQA)bPtZguGp4zE2S zL7bQ4@7DNtwUWQ4#s=B@nG`8{(;g&m*__R`nU!*0eq3bnq>GzY{Slc?t+20nqsMJ( z*LS?@>77P;&Th9~k7&+5uo?)V?+ux8+A-YW_wdQiR17@~C5nxLd9!B-W}f(hDIHtO zQhW3-^~&Tu<~^9`k}`hN%b(OfNG3?bFRVA(K8x{E1-!*ds;-3sx0>Y`D|Ao%-R?^cnKkyJDOCf2{ef7G@x&U|%Y zrA*$Rp!j}_-y|wWXS;96ZV<2Mguv#S-QtBSHjHMB4O;at4^WH zb(kGQ#lIP)zS^|Tlr{O*k^oxCsW=(;v3N$ioVlZ?o57}L`hWUXg_Tuj7 z{R54_2EX;E>RJF7cxB_AzF)7QV0z?X z{?A}qx9kf(;8U!_f;5YokY5&H^e!e7)xTMCv9VA)?Iv=*@ZGZ{qYuJUd*^Od`Us;o zJs-H)%XL5HD4lREHq@r-7kz*UhXnkTV`+f^u5HD@nr2Q)Ha&)eK~&Dl$hd_aP~kUg zmhw|%8~Jv+Bj2sBu8L8pwiEKY#eKsZ(rt|4@$ACO7|g)HJ=2XgHMrs5!svHju~>P} zb49mE=<|4Uhu&6BB!7Oh^CwpYcx;I@c z!Zg~1BrZXv$HD5jzA+Z~H`%~R{*?$#5I_r*SUjt$fAw6%pBZ)5WnYA7$%gsHxk-`l z2=V2sto^spt)!GORpPBY=PmS3VD1Ba0Ch=zIZCe91uILe^(=B zB;x?HiEJ@Ic4!Zpn*oQTaK`w9!=lHAmTNC-`?8}3H#0DxPtE>DoA<|lU907aMsmWa zwtdv6l{b$N9w85rm%j~slYTg-wni!ZqNFOgz6?_y5zs1M?izGh*#%!QYP6vRVh`osN-VvX0)L6y27ej^q?8Z-i7`h+3!lQ`oprG7< zg2h}SXuhIq?pha9Auis~c${US8nNXTzIB}+?FUOS_!RtnpanW}@>zA{*CuE7yux)$ ze$ibzOB~u$dGt`<8`J@vd|h<$x_QtHGcT5Z3pLF#CDv&TqLThQc+ci>O5!ELUTrIW zGoD?(YN8Dr=~afemtqxddVvO)56(2TImWydGr}hVfH3fQ%P3K31pS{@R2Ji z0ydwF0uGQG+i`0!THU)kWUa$gM+y7ba{dRU9RHz9DSt2MpoN*)iLr?$8nBz6e*f{J zq+TERb!alDU)vFU3pGgg`1Wp>n0w-c?RmByXuuqDS2g4d8WzoBb{>~0`l)|u(+xxS zK5J*Hq$}<9S49V1gtfPM_-;61rWb>gkamxJw{qosV|l88_^~PDHo>l7w%rk5M4mv4 z9Y~zlU>3y*$6|}XtorX&Uo-KnrE}>$$^dtZ{Y#2jt2VIDpj1lSVu$mgWK2g5&hxp# zW@&Yih#s=AZLE5S?UmB3GN)f-hFPO+_COOEE{WO5ip@Z-1~T#S zC0Ia{&>#cuwp(MiTuli1&Ss(SOP=%nJHB|$(J^LGwvH^j=Vf*caWF7bv_bPV?oP9x ze2$G-R*u+V zvej-a0xAl%_{DXnnEG@>O0bhyZshJwu$9%L>nkRZoT9p&{VmIWY9$ue8|N!7T4t?A zCPD1GHez%qT5U&Ar$7<4nHp_mhWq}(?rf~Y`FF&;`@q+WPn~?{o90cMxpW3x<{ZK5 z^}9Z)!|VY2po=um*A{N}=rr;gs#{v+TVi*D(QR7A@+eJe@G4YkR%yv?$YVJ= z0I+vURaK98)($_*r7RN_Poo0^B$Z1hnHx2G|E8CZKEoKO00= ztRVEj%dt4r=pJ!DgkRg9@mS;0@KXfpYNedVjKA*dARrOT-Zz!ylAkRj=&1c61tbXG zWUkOuORlDteP)<_wY8|hdY$0DM&X#jT%{L}HVfrgHd=G-@>Zun;t22_g8l`K7}FWI z-tI$bfub8)<~2K{0gb`f)jl-5CVqaUYdWt1A{*L4(AxmSFd(+ ztc+@0STB}cowhpMS|p`i4-_&4f$0_Ov6Zs=(1VqhMI&}Xe>6-$rtUo$xcx85**CgT zH8^$l2@r_q*a)_VX@d5hu%;dLCaU*r=)07Lg6ioYOl-I~qFQG?*apK*X zn%OYPxfgmfoA0Pn(|a>HN0sLN+|Cijo5+;4I<8B1i9FM;+bLFM;lR3T*y`l~kYHQ* z$9k=TG7XKF;A(kMriw^~nSL6nl1iQFHh9fswj->k?(D%_rn)nQrdp$hdMGLM8|$@5 z+Ig+Os4xuwI{6rbiH)nIeu~NK8xBI7c`gN_%FYIswW&`tW0Fj(z(tCWNm?o!hys^h zl>)*;9r`i0~1k5SIfMN`Xwwb^S2*Z4_5>tn?F)8sVce#tbv`o4UgmTh-d6s-S0uAj?&t5F&wq2fZfc(>LJZ75yeG6R0H zM&v=hc!fE9*AlmCVRv`In#mu^mzJMEV&82&At_GDwsF`hyjR*gJ0*2TDtxuPh1K8# z7hk49JT!v{0^_){NNgf)V1n`a?TLWP&c-zV6YBCh7f zQ!7@N|Kw?1U!Rj+%!@Tc}#3P8Gtm8}Q@I$q>0i zqf+Z>*`B8zxfck)r=Aj z(ul;tS?3R+cq=F7<|U?LTNO(UR0h^eT<7qriA>G>6QJhKjT)G6Yrw{HOVXuNFFycU zag}0UZZ0Sge^53)>wILh?G%C!vrLifQxl`lMrV%CRN*R5D+Txcy$PZ>@Got>xy-4> zp3Qd&JLhi{kBr&#(Z?dph;ZJ!^f=F!WY^8-WQzO!qh8LI+^ecw$J&w|?jvn|Q$5|D zx5^Yv7jXW8(O9?{|BhIj5`QaEbbh5mTUoC*vndVG1+gwTzc`g-$~xVh-i8V>qTY~h zz=vo6{O=7CKHRR3F2fwGTcbi!S;5i1K0(gZG1rA6bseu;G^T;uL|fOp4@*nYsq)1w zG%)$kE|Suef5H3O){bJ?*CTmh7#kAEpsV5nDArBCf0(=%aVJF`*%Zr4?>&Q&pJNWL zJJQKr0k~nM3}l^M=$X>zn(Dk4twU!%%oEzYWtqv#uEU)D?YZ{)o2g*5GIJ`%Uv1+N z)CyNe`Xroi{Z;Wvx(cDx80FRecK41DP{b8}b<89eoB~DRp1fqc+#e1I6(PaL@$a|p zL3DKQpKZ?PYpR@WL+?5|xYS66iY$Dj-*kP8sHKwy;>A%ejy*gNn;`cEKwfj1J1cwY z)z^^O*r2=uZr_?2+#huohpn*Tf; zA5tja#|V^AU)}1j29pXWn1;r0)2;yA)WLB3Hj6ctdJiBmZp#zG|G=m0=ECMZ>4xAs z^7B&j5Vmqkef*%Hhpi~>H1I^Z0g&$QyD~gnj%`n|+myW>6s1f2Y)BO=aIT3?>6yt| zmA*wh>nX0Im&FI6Imxu~b0urwuXMkS3mN;e__$Xl%eDqJxmoCgn5>WhCw@6z*T?Za z_C*oYLC89U9NHc8j3WZOxp|z_!9WUeF5DLc4)C*GF{ApI%k=;NxLJm8zEsUyIA1oI zxEGo&tybcuzLSIx@7_HQLW0V z{1keYg;Zy%@X4XS85VszNYKyQs>T?Tnj z^0ij>_Iys{nWu?_bJV6m4x!h_nf2R}4P(C*f6NYa(GiQ=g+hgKBqpvhecHZaBWU!K zpa_L;@~UsY#N_ zdv7?6q^5x4*6#W_(tOy3&WFef_nu-KD@D@7J5CA$dAu6e_w0aHVMC>F5%m>CKRM9( z$E;KZb;Nnm+^-*H&x3hCr5C2xa#(3t$hOtjqYf0<0^dQ_MY^Bs_5I#3%)sWM=@0}a zH)#!oyC?wFrm5jA0scQ^jPuotS?HZXMW1BT*+RYW`e4HjFBr|VPnR^3CykA-Kl*EL z#6CSqolMq;!=+0-cx-bz{hk0mtRla;&n!FrU8RM4>w-TahoLwLf6h6)VgvSXozIyz zSZ0e?q5hUo7hX)kZ@n=O)%K7&*sN#Czoc8XCFgV$?r`TI_ego7-y(J`E{BS5|C?PG z>JnVSuKbCQUDry99mOA4>Ji=KPGKY|0v;cSkHr%HUa86j1phgDx`m|zP%WNa$>_9` za3TbA6uj-bNm;mIE+$3AS*do;+8uq8r92ICR2~_L_!Rx6X~0R_2)v9IPH4~vbg{Hl z66u#MklNT*$U!G6K3G~KLe=lf3*_j&fxoEgdxHZQEuu{?&0sTTS#uu;c+jWVn{76Z zJ)yEd%#JEH?x>fVmUOuSPA_O5c?Ar!QeURGV^2Y2it}a;p?SU@EomZ_VH-(X_HlVS zYM!PMN)3c=vgnGIT2Kb{K>+W2(`iNYC8h`H+9I?w`PU{dU2SABVI(A$IYqU7_}u%? z0gmL&IpJylmvbX=uA@6Q`Y^?<;Om+{YM#vJmei!H&xf3!1s#-$PcAkO@>eblxhtvt zS$%J;3(uy>@aA+-4M%a-HZv5$a-)!Y-W^#^d{jI7S58TDjDh9_#a{K0Dp>+kT!B`V zrdmzp11CNk*Ae{Y{EQC;x3C@Fv(tmN#J&54E%HZTq386_dzcEB!=oruuQnW4hrTz7 zDyG*-v0HJ0mlsC}^+P0d0y^5w$>8wP zIT^s8MohD?c-f6N#~byF5>3<9VYQ`ec8e+`-WC~uf=z=B*pAq4{yZkn54G(<#NQuO z(!-_Mu+XDA=>1Gg{Z8x(==0GhF}k&4?5#=xmoKy5N7Y&-FZ=t~xRfN07FKMua+Y?# z1Ws`3XQ_4C<#*&JnSu%pgeG!*tVvv#MmS?KYK8&2IhVB{a0j1l3?hR&T%m!XTl-s# z{YjIf!Af=6UWy&L9gUh#W-HXY)pH#wO_eB#X3juFPT(?E11QSy2S!QOw(jyqTQ(rDrpY6-Cf>jRuAEMylPZ zuH<#bd1GhR*5xQD{^ok(Gqh!O|>2XhJ%1H=903hF_Rp z?5RzR9hCwtprRk3-mtMM5gNo!_`9D6P!-WmDVI~$cmL%{iU;0gbiuZSk+fIRApT*2 z0slrG&LX_IDv?0j8{DwXafC{g}YlolrL}l&#Cunn5n&rzNJ<_WZ z`WO=+Oce{NwPOjtO!nz* zk&!S5-niB;W;1^~D#Rp4`UN`$wQ;MH@J~m%McL366z8md<(=(Z1IE(}cs z+d2}{9{-KLo)z!v7`l~POqXt>^E;jBUH0J|5q?)v;b(^Dt668QlgB@6R;QpGubuE| za|JFCGp8DYH%$+%5|i-U6wked+z;>Bc2_%u7S4ju4yWhMV6>S%qo44RV!OL)n)cn5 zVFjZ}8BM*_`!$L@Stl#AM!IRA`93u14{@_St2+-jB#Wu;y5$CAGC{_`tn*`;xj`Kg;4o+Z)uUVa&jQEKcAIKzP zLRwoq_$PkV&=^VNmZM|q*_;atI;r)$YI}z_$?%H5IpEk=b7gEvus9wq#MGJ8#DK7L z4UB!bsrp2a9bmwQyYP0{e7+g~z3O4_ptuFUkT+StCME!cQgT6_FE?;#!&#%NXx(31 z7fEO!^krm)7$tfT`%NpD5L}kZsj!e%#U$p#PKKN}hdX!lQL?|DtpS;Z@IQp#!dI-W zRkHoZxDuGG01w-a<*W}Q(9N3V*WMs|02~G%jRM1l`i0^4K|QPkQ{!VgzP_2rVi!HBItZXINJEi`6iVTL@odXDKAs?QbfZ`N-ejT)|j7J1CU zd-tumLx|Cc;62PUV8Ivld8sm14`0`?193q>anok0)bEFKJ`h$Gt>5YK%Bdr&vM+<1 zC%El%TwHkQ4v>{}QB!+Y>peyo&f;bCN|k;A7fNz%TNv-gk1jolo^7K1p7rl5W=MnR zk2b%*XTyJUS1|RGJ*+;&coonPT4YCs;s?^_E$p)E10#y`3P%3C;z;m?0R>di*8JQ$yqhmO)?+X+J0>_2&HkJG;hw3nxG;UZfRc7h_8Z!V+OZ_MY)%h!ANL!% zb4m|`jKLVx99{D!PB@8kMPyRWraP+ldr}&h`go!B(d7a52yeGG z^7AK|ZzOw4zpXI)-bmv(C}!_!oU1hP!pdpmxgXyCj^ckTt_VBQ zWbI{{#_kymn{rhXiJ9-|vSDw$Kx6opQ&(A|B5Y&|s}P9@Mu;v}a}@%2+=$42a!H74 z{stb*@>d?dFuw`<5u*r1WgFu+#y1>y2zJ?)C>Xaiv50TaWAKv%|0+B)5N62`d2@6J zu@U1&*FR&D)`h2N7v>40;aSTsz}-yBi2O`CyF%AFT0%8I-#(m54p3j;Ys870kd&Ea zIlXU;vIq-?5;1m9koTPu!dp%-=Thq${McfuyF*V&EOnc zFp#_R1`lP~*VG9QHNCV{o`1ro8qe}R>iuS&5}FXS~9FYMM~TI}2`3 z*u=be!FC--H0YHaT?-~6)ebsY)yIh!kkdsDHb%M#_xo!};p84Zs)dc)E6O&)(~GvJ zuQe571lJF{B7Gr@XB#DH8_iOC*W~K(WtDl&$1Li34+`zyEbG06((mddg58W@L9+cL z3lF>E{w_t2?0^FRkjJB|vE1^wruFY>@;`Z*ZXeTlV=$iHZo!#fSBE(AycBhOZ6%DE} zi(4|$E2yIS?~PsY`D;V)y^S$QBPUV#;=WNE}HPD{_B zfclhG*_DR--Mbm1^5purlzwe!y!SR1=xreG27)9w zUpDt$@0cB`6W&m^JCeytN_Xal99|CY1bO7Gw#d<<$F>4tJ5_r&dv_iJ_m|vPDraua z(2MkXnmyw7ZHxO$j|acH8q}n#a|3`rsDV|wl>Y}&8p~gH%A9s!F^|k~{WxsnYUC(z zml6Zo7myjgy*Iu!lRe&ed(&jHDKST;orK=JMCUedGVD;MLqze7Ych9oI!8rz0423@o$h2J=7sI$F-TVbf-~HxgpEKi$L+>AA?1l>PvQh5vGq}XHpO>r7v&G|{7e0~8f+OLJ(_43}HZO!23lIJ> z`>~LcJEft<$b~mZe3`NYdPw&2#!>ES51eW34!Bmww@{Zm;nw9>fc$6GAFvpR5hp!( z(?e|RLH_oJLV!}c)IsRJi!!ds{*n=USx2V$e1<>DopR1y%l*YmvENC8r|5FwUufqkG>A9=&kK#aZyK1zUJDZ+R zZ4GJv**WHO5CH`-(ztq4s}Lnl0YhprDLiI*TrYn>3N;F8lKQdM3a5-6Vvfn*O&0*K z??ib=UNE+bWZ(FHX>98wtjQ9$EmKh2shz>|<~_A0B9G;dHHW-j2Q7(ncTW*6u*Jt5 zpFh25<6lM_?aALLXflvAWmIy<;Ub)|0Kew)G3nngdJI8g&N{a7kb{7z)(bAl170%3 z&Nqwivf5&2-xyT;sVAeI5yFAqtJjYp(63+*eHhzITQbX1dK@x)_c(y%lGo7EV4$s) zz>}XQ+~E9yjkkasKRz(|Z5)!d=HVNH&)ol|=zTMhl;XvX#N(N|z4vIiwO)naf~^3? z)_2LJ^xR`?p=(VNEWO0?2mWY>L79o=MS4TZ8e;xDTcT%<0RaQ{6&B3RLDrT3^mo4F zp>1%O)$9tiJzA@__QPrg<=L2nn=2_gTqT`k=bE97CFjE5`hz-u7rJ;$UjhyI#NfIC zcs@O3#0)UqjFSBUkSyRT1Q;s5P9`4#w8whm35-eIk<%D zao|*^I^3M5PF%!tUBRX+uZIENdwDUWZL|U6!B70H-3O3G?tG0 z^E3MO5WZUc-Gcc+T*AuQeL}d3}I|rnddw}4@(aOaeGLAYdhdKczE#n4hgq628!TL zm$aBHGidA`dSq{4iX!g&SpI%c1!0FyhdVoCA&`ceMDFk_Ect@mmwYr@DT(Mo`HE%!D%S|rm z@N(e|j(8)=*)pSca5eP3S!7}>V+_750(a>S^J0q}7^MbpvWzInXiHIq3lY?~r9#{t za;~b=J-OdxKT5GW*Dx&LQ71TmS~+23b9J4!=IGyaiER-1A{lO-FI6iG+(fn^*Ivv4 zi^V)k1w>n1eEV!?O~t#Z;qd;z3G+fs)_@dyaX?Ypnfn>$0%558E(TB*@jOn9%qgy4 zw}Q6u;bn#M>}?M2x1{ksgZZS9FpXgR@Jf~g0*zTbb9!E?hTURxZ;BOm4snNAaqKHw zk!TrzeP_K$N@e^jR!O&3 z&t&M?d*)Ft!N~~g*9?_AcQwLkoeiVh`iHy~&=4)%lvku<@N&6SDC#B3$QvKbHr5Dai@W zoyJC--+kjb)TP4vZdRKSiZ<8v{n)Q*GCu5C6f;iJDJH{{N5|!D5Y!ARG38lP74glYtiCwew!Pk{^?ZawaPs)TnGNh!t;B$W^2Vr3 z%?QZgp0Tw-T39Vt%JU&P%OLwtSgT~_(nqwagq+hOb|8A>0teH>x^NF(C3S_?n~oG? zvY=TSVC0kGWCdBoHDVNbx8F2a&;mw%vdu(=d02K5e0wU3z|X=pQen7dfzTz1yj8jC z8iUk7wJnmYcbkfpKTm?TZW}0p1UCYwdd&=jG$tzHG7@M-dJzTCi8`D;Rjf?Hls^d+JWtK@OYlW)1wG^z}JYgW_7?R zcM-+7a?wKZ|C*2#g!=N-O^%?@sPpZ0h>P&VucUlA+9B@GX-SS_`=<=&u-5l#%x?0; zeJapq`!H{1v4T5eC*fo{_(?roS>0kRIq=`j7VdZ~wLcYCx}WNQY~dhAgx7n-(tGLk zmc5#Y%_wq5D*BJRn~HWo$=JQoV+-Ch_E-`vZo3QbBL(TiSq_0wfhFg6o>On$9 z@KR2Ea$&falrm&||jBmfl-cjoP^GOpXtxd>ga)Fgl5TdN)xx&EN=^TOOe!huAqX zddR)OhKv6XvV-|3Zvnd2d<5znQjHK}t8yH-_W{{&BBeC16LeEo&f8gU>tHwzY8Eq( zntKK*RAHO+ng2lg+y6PUV8=g4L_S-)zV!_iCKq~LKTzi|tE27G%z-qc)hnYjuaJjW z*Fnv$aDmkRcjYZxjY0-B=)~N5>HtD7FP4tP6oSiLBUS>JruO(|>B7P~T{^nczYDAX zP{^UQPI+#+AlL)Ho9-^Y%YnK&cj?TNbG;P#ylf%#Eo`$0jX z+g0jF>(-Rl_*eZLY-U~D#5TjA#IM}GPILhzQHhm-GIA~+2sKmCU2!XTw5Sl@0m zF8n|-I>F7 zyY&N=^LM#64S3X)?AP1F3_EI4cXb8ZiLHgk6j)kkL7$qwi+&S8!!PJ;=tKCI+fx=; z9`h7G-N=yO@kkXQS4CHDXev{+a?a5MI(k40a+5hJ;o>1NueUbUI@z)a2O?ECL(|Q; zmhjiria;GKcs2VZ$^OHy=}}%ox=4!Dn&r8DXJxcn!JIbdXrAYBL?H3P@86cMboW6+ z2(%9dYsOdaDwnB_gg>}YCxGy6%*Q#6vn3lHMb`1PLhLXRWFUh^;>u2{}5PC{b0VwXO#OG_N|FmY$Ya>xYqUwsL#h{BfMKOn4kt8H78 z6MhGx653QnJn+F6f036fbWJsS>y6}RnNnnR-Kby!S|oNG$dgwgyBN>bpQ|@tZ524n zI0PLTWTC?ChMk;9D)$VlTR-~@Ix>HgN5d&Sj>P9l(Jw*uyr)Bpa-1!DAp>nzW$3@Q z=ibM6j$kh>)a(n3f;cn!DxYemJ7t&qJo0E9n}GX#T8=>*P!p`_6^Wz6FI-Wp=9{39QWJ4N z38p_y5JU@aM(v0tvNZ+7+%Hxh4%-$H`oT=2F>)H`~q4(&BWSo~OHB0a!Na{*n z7`F=38RpyAh(W4OFyeW?uRCai?hWdeLltesd2a0R?J3i33q95dWTO&)|&R9^fU95!`l6#Q}Gb;s-kgk~Q z>>{T`gUz?lYGF6qA<`<@jT+vKy0!KnLrNR{A5;A_yZPIk6vY^D>pKxT+?+%;Sco5lEGc(cKtZrJ=lj2yI^Jc7Aekp zJz#K}YdwEmx7^VuIj+#omvA&AAjgegRMDU-Ie9e5gt}e zVS*9O`-(y&(O$+~+S^=Ydq!5e?^HX9vAZ}F!E0R@4e#r*JY9&<4wLUPkn=7~K#H(9 z=M?PCdZ~gzdIz)L!og|SapGa#@q${UU`cZGzWKR6e7Hg9Bd6pHB|<4~j?u7ICGssR zrWwL=C1Z>%IFd7(V9@u)@2#V-QzMLmwh;RW1LWiF>L0lZL#JjM8x z5k`JT6P1@7KdPX5(ug|WFcWT4oFU?zyQ+c01>jDblia$nrFiBvkluOEQ@0QQ%ohD1~^T5u&Xcj;B(l6IUbI$5+YA{y*LF#SQn6?u5+rf|}MJon{a_*H~ zU2T?}Khp-3Q-%xpJiRE68)$D9fU=F?QG^eVgP?->v42cpRZsHBp$#SSKkEZRGwvL- z_L;kZuKG$3{V{(X|Dnfsy1zqLk2RkXtSQ|(Prglmi{Tq?y`u!F1BUYgPtC{< zwVxI6K$&-cQM|dD?ICrdq=%MOlNmaI=DGCsH)&@Tbw^DO{I@zd9rMixud^58=+n8Q zP?#FBIgSl40uTFu|0Wm zbF}}j!xSujPnK5DHc!18dQrp~abhusaC0uDWvj&sGWA!rW32J<2lE#yITV_ey7dgX zQphuMd@x75xzZCF`u%cxkuvNS8J$(TS>gOrw4B*ZiXP$l9g<6a20p-tfvA|7W!#t8 z$>4HpC-1;BXF)k!5QGJAt>4>tAmc&s^5gn4%UW4h1qC4@#R ziiM~G3yWc2rj1t4`q}5fs$~CTXa%Fa-s%^8EEn%_Bdi9LdxKEXp{#NB>R>9V@YP(v z3-AF2)(|;AvhdOE1?t-G%cg-or%UkHgar*0zFEwQ%Ffwpskb3CgS4SvotFA~$@i!1 z=U9f2<^b1r(_9#!>P5=0N86+_SQjz+V6XNff@vN_KV5{+Xv z(AuLeD@BA%j1~S{(Nbv;_d09=NVqfsZzIZe3U-3@{}n1a@bc4bzY)%}=AqwToo-hg zBRsiRtB8zJ2Lk0PqQr-l;DuzH(Pi&544_rvF`JI4-FuUBJ_`jvUu&x6vg*J8$(T1b zk#?%=5};*}l5INkpppA3!rwP%J4py#DJtMDh}zb^rIpunmH+Tm9w9AjXc-V$!?`q*5>UYzA5J4z$Mz3QQi z6~D5?X&P$%B;9^SW#zfW)aWoNuWh=BTp(xQC2mB|pmbbuHn;M6yzt8S5<6(h{*v^r_brXbgZm_d zkwFxT)&w3|v(6*xH97IIK#BEJpS#oW#?%k{qbDN1V=uSh>99$+%2`~+rx&-bSjd6| zT~zVS6tdT3`&JN4h4J}u4jx}gDSKAss2nX@lNFLrBBEZsB!}*ER84W_j5QoC61Y0& z7yk^scjN75BH*Cgw?1{@oyfkl!Aczqon`LIO*`}3O|E+ zCCs-T{P(3|!_j~ChYX~s7%mjdRDbAQ>3+<9g0D|yVazBida0Pv=URi&a-flKqjZOI zqer0`$J{}2ji-1@kbk3{?!pVznTKDpIiKdk#y3WtGd5#>wX1pG3oYg9{g55k<(XPr{CnZ z@VTXtHEe%~in$XVIlhgzf&iA@TZTg`*?Dh}?t1%hEHtB}tosM1PWPfq6b9JKqJ@KWkUbvSO`NrS?^f&-}Wg9avQD2tPWk>K&;rM2#$@zJQ zzSQ0M>aommyA>c=Gwyl69BDGHFNOrMeeGnHJB|POlAXJd04g@Krg~oQJ_x4D_jD0C z?NRm9L?xL7K#kB3+Bh@6AB`YyRNZJ2e>5`ERL48I97p@-NsI%zD?Hn<3JS|g&(M28 z6x*!i8R^>+(mD!)Q1~t59@K-Hc55Hhjoz*ZG#>KbyKFoH20X4z3u#6DUSon$o}j&~ zk0cP(Yl!?Iv)FUNVS2@177VYNPTUhQ?^c@sA#W&#cC3{JPE`PBDf30hv-g9Oy<14^ zVNXLs)ql1Im~#38r=|L_)FoPDisI<8`?h};b*Y+(mjjOs8@^eT#EAyEA->rWen;Dn zR9$SdUCh);zrDRTQHdM2S_L94np`-#jHLRP2Hj+pe_s`!-W&l>g#f>pb_js1ox0g| z-^ydky~Zm1I8;7&AKTI5;bh2FQD0lky}|e2U}g9H($dm%)%NmKP-%f~2wzA!0CV;@ zob#s3A_-C(>DVWs@0U_h%5PgvAWMSF`I?mmD;}`6vhux!2k77bX`=ahu)o;NnX2~m zWYO@;_Zj!SlxEc=&Ri;l=|z?SO)P4q`qGuU?In8-HPSqkspkJ71plWd^t}TDPWSjP zIf(jRjT@pZ;6JD96k{yaJo~G}t$pfetUb@33`V|R&k;kH&Uzk@PX3>R?PYx#@<<#F zUdEQ$kh%=QDDI2+HPKcQ&!pN{6e8wanf1d;ZPw;t5<~K{|cW z?|J+r)2vwkQ~gzRW+Jl~UtisE_{k9Ai&6iEfb`_Y`uW25s?R2^IO@QvDXwkP1r4x@ zZATpD*x)%`T6K#<^PW|=aPC9-6o?+4gXU)}Z){-jS z3&2?S@7_ZBwgYuHdq=`;#OT3uoqICZuIfo(b;eweA6^p(q@TcFVqFcI)l=siRYBDi z7UQ4Vk&|Z6)Gr{!+*O{`)MHz$@c4pSpRy0MM`$TxM*{1$HR2ZF|eVepXRzFduOqd)3l2LuAd$B8wG(SIlE2Q#eVOO5QxFM zQ2ZaCPz_wY95K*wZfF3-zu}$UL9oyLF261>8ql*IuJtdb*tkc;v;?uAR^Tf?pXtAk zXz*A*g2Y6!agZ8JOoe;N9~x&+FT|>&xfk%Cluj3*&2044UFhozSjqaGll4%X$Hbi; zo6$Ya!O3NrnT=kzp`z#^O82KsR;|ExkVGmv?Nn?RM-1 z>{{hQR;bG^jkvN)+qSA7qYe0(-wBaX2T=E;w*^PLTO{|VjwRX9s;k{&<|Q_w zgD=DeBFe(*>ZjCL+Vu_xJ?p36`g4zMb2LyI7KQpJxS?`cD#cxO4Lic7U^*vH{k`t_ zP~Bo%uF~SAJi21x!L@q3U5$uV{oO-@fub-##l1i!rWM`#)bWpZatX>lVfOg#*wG;O zJ%0UCSdn~?R7kLUjI!u#qek#`zkPsDU*VxV2~M^(_(91vqLfbW!@-WKS)Tw{F*-T( zbDZO>($nT>f4y*;bX-vPvz~&(nc7CgDoCD($Qf+9D!o4QyT*g#1cko&yXAV?q~(i` zpi9-oqg_{;G_uMHnU7qqi550VL-GH0)O@~Ldz1t_-N@Ed!2M|%e|#q)B8t~-PBmrM zRia`HEzH%gBUgWV*A?@iNG5Y4GW*mXZHn{ zeFs9^6(;IZJN;pU$H1FL@73F;W!6BK=Y<*c;n{%4S}nDy^^u#d2}*Y*_}I8T{jP0` zw_l+~ADr-a-8>1p>ue9>K118tepLxb~whjMEQJj+`P<;Tv7b>=@R*H*#&+4A@t&+}%J zk5&1u87Ye;Q>YpFmarSRf-d|*m>=lo=gX@tnRO7H=_4*~I2&!4V!D(PZ{0pN^lR*h zc?53qKbgu1A+bdYSabAjD&=#cj9iJeyq{0fB<RsYJCu z_=2~&~ zEo+c;#PyxK9nR7dCcx8jnsp^-OwUIFWu?cOo_f#+5?NIZxU&+YQUI0$f%VHUZ141n zOpH9hDO7|W)i4Q04LF=Tq3OR}P_U3gIph^8sW&&k?DBo-oYVC;dt^MDYSyagrtitN zO+hg|*m)?^K~twk(__M^ZXB0lHkr2awwr7{d#2gHiy3ythzP(yh8y?)y&Z>RXqt$9 za8xzl1XucHYRoxRb=MA_DyED%lja@s<9SkD_^0S@SVp}QOE|o8_d@7CS2kxZY)EyI zqtoeoW|rc6`4}+Eq%1UczaQRSv!NT?vCUByTdZ9>nI#y6XPK839Oa9j2e@R`ZgNKE zg1R`daXwFy&1xAztz+_5$Li>+k zX}&|ncL=L3zAw~OC$9#*LK*0L_txiEgX159E*8#9Ke~G?La(BZUO+YjBF=f3FaUd#=yK9oOiSXXz? ziSP^4D7uCC6wl*oht-zVdsUBpDO*ej_nf6VIMsej4Fq_mub^^Yy4QXH^de2 zZSlNMxv(G%G0W!I55p{7b?3_|>Cjy-&#~fK9m+z~tinEUSyVX>TTYDVRN*-cyIboO zN@jqx@`CyBKzp~w*ST<)Ia(Ebm5Uhw-(&Z zmghf>9K8-0FHq=SDXN;SPj|~__XRyM=33X)Wv;$pL)8<1vO{-%D_-sjr zK)t$PO==c^aHBUlU*hhUy#C?lweO!9*b{iP57}u0E3IN0uOgsCr0il*X}NcU(QT`o z7<}_Ys-NnK6%&6I2BmCC=RjeRNSlLs#z4E`?37x0E#=ANi4{B$peWA%OYE#Eh3-41 z|GxWwk?8$zu&Vt(goOY2>i#d0`2GL$<^3-P{{NQcAFJN~B__oGO=`OTV@jmh%oyz6 bzl+ypu4@=pQ~%ZvX35r1Q~ literal 0 HcmV?d00001 diff --git a/doc/doxygen/stylesheet.css b/doc/doxygen/stylesheet.css index 1562747f8c..36c1fa4fed 100644 --- a/doc/doxygen/stylesheet.css +++ b/doc/doxygen/stylesheet.css @@ -26,3 +26,11 @@ div.doxygen-generated-exception-message { padding: 5px; font-size: 95%; } + +div.image img[src="fe_enriched_1d.png"] { + width:500px; +} + +div.image img[src="fe_enriched_h-refinement.png"] { + width:500px; +} diff --git a/doc/news/changes.h b/doc/news/changes.h index de2cc39093..9859b2b5d2 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -206,6 +206,12 @@ inconvenience this causes.

General

    +
  1. New: FE_Enriched finite element class implements the partition of unitity + method which allows to enrich the finite element space based on a priori + knowledge about solution. +
    + (Denis Davydov, 2016/09/28) +
  2. Improved: The doxygen documentation now contains nicely formatted boxes containing the text message of each exception. Several messages diff --git a/include/deal.II/fe/fe_enriched.h b/include/deal.II/fe/fe_enriched.h new file mode 100644 index 0000000000..b3131e46da --- /dev/null +++ b/include/deal.II/fe/fe_enriched.h @@ -0,0 +1,613 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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_h +#define dealii__fe_enriched_h + +#include + +#ifdef DEAL_II_WITH_CXX11 +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + +/** + * Implementation of a partition of unity finite element method (PUM) by + * Babuska and Melenk which enriches a standard + * finite element with an enrichment function multiplied with another (usually + * linear) finite element: + * \f[ + * U(\mathbf x) = \sum_i N_i(\mathbf x) U_i + \sum_j N_j(\mathbf x) \sum_k F_k(\mathbf x) U_{jk} + * \f] + * where $ N_i(\mathbf x) $ and $ N_j(\mathbf x) $ are the underlying finite elements (including + * the mapping from the isoparametric element to the real element); $ F_k(\mathbf x) $ + * are the scalar enrichment functions in real space (e.g. $ 1/r $, $ \exp(-r) $, etc); + * $ U_i $ and $ U_{jk} $ are the standard and enriched DoFs. This allows to + * include in the finite element space a priori knowledge about the partial + * differential equation being solved which in turn improves the local + * approximation properties of the spaces. This can be useful for highly oscillatory + * solutions, problems with domain corners or on unbounded domains or sudden + * changes of boundary conditions. PUM method uses finite element spaces which + * satisfy the partition of unity property (e.g. FE_Q). Among other properties + * this makes the resulting space to reproduce enrichment functions exactly. + * + * The simplest constructor of this class takes two finite element objects and an + * enrichment function to be used. For example + * + * @code + * FE_Enriched fe(FE_Q(2), + * FE_Q(1), + * function) + * @endcode + * + * In this case, standard DoFs are distributed by FE_Q(2), + * whereas enriched DoFs are coming from a single finite element + * FE_Q(1) used with a single enrichment function + * function. In this case, the total number of DoFs on the + * enriched element is the sum of DoFs from FE_Q(2) and + * FE_Q(1). + * + * As an example of an enrichment function, consider $ \exp(-x) $, which + * leads to the following shape functions on the unit element: + * + * + * + * + * + * + * + * + * + *
    + * @image html fe_enriched_1d.png + * + * @image html fe_enriched_h-refinement.png + *
    + * 1d element, base and enriched shape functions. + * + * enriched shape function corresponding to the central vertex. + *
    + * + * Note that evaluation of gradients (hessians) of the enriched shape functions + * or the finite element field requires evaluation of gradients (gradients and hessians) + * of the enrichment functions: + * @f{align*} + * U(\mathbf x) + * &= \sum_i N_i(\mathbf x) U_i + * + \sum_{j,k} N_j(\mathbf x) F_k(\mathbf x) U_{jk} \\ + * \mathbf \nabla U(\mathbf x) + * &= \sum_i \mathbf \nabla N_i(\mathbf x) U_i + * + \sum_{j,k} \left[\mathbf \nabla N_j(\mathbf x) F_k(\mathbf x) + + * N_j(\mathbf x) \mathbf \nabla F_k(\mathbf x) \right] U_{jk} \\ + * \mathbf \nabla \mathbf \nabla U(\mathbf x) + * &= \sum_i \mathbf \nabla \mathbf \nabla N_i(\mathbf x) U_i + * + \sum_{j,k} \left[\mathbf \nabla \mathbf \nabla N_j(\mathbf x) F_k(\mathbf x) + + * \mathbf \nabla F_k(\mathbf x) \mathbf \nabla N_j(\mathbf x) + + * \mathbf \nabla N_j(\mathbf x) \mathbf \nabla F_k(\mathbf x) + + * N_j(\mathbf x) \mathbf \nabla \mathbf \nabla F_k(\mathbf x) \right] U_{jk} + * @f} + * + *

    Using enriched and non-enriched FEs together

    + * + * In most applications it is beneficial to introduce enrichments only in + * some part of the domain (e.g. around a crack tip) and use standard FE (e.g. FE_Q) + * elsewhere. + * This can be achieved by using the hp finite element framework in deal.II + * that allows for the use of different elements on different cells. To make + * the resulting space $C^0$ continuous, it is then necessary for the DoFHandler + * class and DoFTools::make_hanging_node_constraints() function to be able to + * figure out what to do at the interface between enriched and non-enriched + * cells. Specifically, we want the degrees of freedom corresponding to + * enriched shape functions to be zero at these interfaces. These classes and + * functions can not to do this automatically, but the effect can be achieved + * by using not just a regular FE_Q on cells without enrichment, but to wrap + * the FE_Q into an FE_Enriched object without actually enriching it. + * This can be done as follows: + * @code + * FE_Enriched fe_non_enriched(FE_Q(1)); + * @endcode + * This constructor is equivalent to calling + * @code + * FE_Enriched fe_non_enriched(FE_Q(1), + * FE_Nothing(1,true), + * NULL); + * @endcode + * and will result in the correct constraints for enriched DoFs attributed to + * support points on the interface between the two regions. + * + *

    References

    + * + * When using this class, please cite + * @code{.bib} + @Article{Davydov2016, + Title = {On the h-adaptive PUM and hp-adaptive FEM approaches applied to PDEs in quantum mechanics.}, + Author = {Davydov, D and Gerasimov, T and Pelteret, J.-P. and Steinmann, P.}, + Journal = {Journal of Computational Physics. Submitted}, + Year = {2016}, + } + * @endcode + * The PUM was introduced in + * @code{.bib} + @Article{Melenk1996, + Title = {The partition of unity finite element method: Basic theory and applications }, + Author = {Melenk, J.M. and Babu\v{s}ka, I.}, + Journal = {Computer Methods in Applied Mechanics and Engineering}, + Year = {1996}, + Number = {1--4}, + Pages = {289 -- 314}, + Volume = {139}, +} +@Article{Babuska1997, + Title = {The partition of unity method}, + Author = {Babu\v{s}ka, I. and Melenk, J. M.}, + Journal = {International Journal for Numerical Methods in Engineering}, + Year = {1997}, + Number = {4}, + Pages = {727--758}, + Volume = {40}, +} + * @endcode + * + *

    Implementation

    + * + * The implementation of the class is based on FESystem which is aggregated as + * a private member. The simplest constructor FE_Enriched fe(FE_Q(2), FE_Q(1),function) + * will internally initialize FESystem as + * + * @code + * FESystem fe_system(FE_Q(2),1, + * FE_Q(1),1); + * @endcode + * + * Note that it would not be wise to have this class derived + * from FESystem as the latter concatenates the given elements into different + * components of a vector element, whereas the current class combines the given + * elements into the same components. For instance, if two scalar elements are + * given, the resulting element will be scalar rather than have two components + * when doing the same with an FESystem. + * + * The ordering of the shape function, @p interface_constrains, the @p prolongation (embedding) + * and the @p restriction matrices are taken from the FESystem class. + * + * @note This class is only available when deal.II is compiled with C++11. + * + * @ingroup fe + * + * @author Denis Davydov, 2016. + */ +template +class FE_Enriched : public FiniteElement +{ +public: + + /** + * Constructor which takes base FiniteElement @p fe_base and the enrichment + * FiniteElement @p fe_enriched which will be multiplied by the @p enrichment_function. + * + * In case @p fe_enriched is other than FE_Nothing, the lifetime of the + * @p enrichment_function must be at least as long as the FE_Enriched object. + */ + FE_Enriched (const FiniteElement &fe_base, + const FiniteElement &fe_enriched, + const Function *enrichment_function); + + /** + * Constructor which only wraps the base FE @p fe_base. + * As for the enriched finite element space, FE_Nothing is used. + * Continuity constraints will be automatically generated when + * this non-enriched element is used in conjunction with enriched finite element + * within the hp::DoFHandler. + * + * See the discussion in the class documentation on how to use this element + * in the context of hp finite element methods. + */ + FE_Enriched (const FiniteElement &fe_base); + + /** + * Constructor which takes pointer to the base FiniteElement @p fe_base and + * a vector of enriched FiniteElement's @p fe_enriched . @p fe_enriched[i] + * finite element will be enriched with functions in @p functions[i]. + * + * This is the most general public constructor which also allows to have + * different enrichment functions in different disjoint parts of the domain. + * To that end the last argument provides an association of cell iterator + * to a Function. This is done to simplify the usage of this class when the + * number of disjoint domains with different functions is more than a few. + * Otherwise one would have to use different instance of this class for each + * disjoint enriched domain. + * + * @note When using the same finite element for enrichment with N + * different functions, it is advised to have the second argument of size 1 + * and the last argument of size 1 x N. The same can be achieved by providing + * N equivalent enrichment elements while keeping the last argument of size + * N x 1. However this will be much more computationally expensive. + * + * @note When using different enrichment functions on disjoint domains, no + * checks are done by this class that the domains are actually disjoint. + */ + FE_Enriched (const FiniteElement *fe_base, + const std::vector * > &fe_enriched, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions); + +private: + /** + * The most general private constructor. The first two input parameters are + * consistent with those in FESystem. It is used internally only with + * multiplicities[0]=1, which is a logical requirement for this finite element. + */ + FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions); +public: + + /** + * @p clone function instead of a copy constructor. + * + * This function is needed by the constructors of @p FESystem. + */ + virtual FiniteElement *clone() const; + + virtual + UpdateFlags + requires_update_flags (const UpdateFlags update_flags) const; + + /** + * Return a string that identifies a finite element. + */ + virtual std::string get_name () const; + + /** + * Access to a composing element. The index needs to be smaller than the + * number of base elements. In the context of this class, the number of + * base elements is always more than one: a non-enriched element plus an + * element to be enriched, which could be FE_Nothing. + */ + virtual const FiniteElement & + base_element (const unsigned int index) const; + + /** + * Return the value of the @p ith shape function at the point @p p. @p p is a + * point on the reference element. + * + * This function returns meaningful values only for non-enriched element as + * real-space enrichment requires evaluation of the function at the point in + * real-space. + */ + virtual double shape_value(const unsigned int i, + const Point< dim > &p) const; + + + /** + * @name Functions to support hp + * @{ + */ + + /** + * Return whether this element implements hp constraints. + * + * This function returns @p true if and only if all its base elements return @p true + * for this function. + */ + virtual bool hp_constraints_are_implemented () const; + + /** + * Return the matrix interpolating from a face of of one element to the face + * of the neighboring element. The size of the matrix is then + * source.dofs_per_face times this->dofs_per_face. + * + * Base elements of this element will have to implement this function. They + * may only provide interpolation matrices for certain source finite + * elements, for example those from the same family. If they don't implement + * interpolation from a given element, then they must throw an exception of + * type FiniteElement::ExcInterpolationNotImplemented, which + * will get propagated out from this element. + */ + virtual void + get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + + /** + * Return the matrix interpolating from a face of of one element to the + * subface of the neighboring element. The size of the matrix is then + * source.dofs_per_face times this->dofs_per_face. + * + * Base elements of this element will have to implement this function. They + * may only provide interpolation matrices for certain source finite + * elements, for example those from the same family. If they don't implement + * interpolation from a given element, then they must throw an exception of + * type FiniteElement::ExcInterpolationNotImplemented, which + * will get propagated out from this element. + */ + virtual void + get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const; + + /** + * If, on a vertex, several finite elements are active, the hp code first + * assigns the degrees of freedom of each of these FEs different global + * indices. It then calls this function to find out which of them should get + * identical values, and consequently can receive the same global DoF index. + * This function therefore returns a list of identities between DoFs of the + * present finite element object with the DoFs of @p fe_other, which is a + * reference to a finite element object representing one of the other finite + * elements active on this particular vertex. The function computes which of + * the degrees of freedom of the two finite element objects are equivalent, + * both numbered between zero and the corresponding value of dofs_per_vertex + * of the two finite elements. The first index of each pair denotes one of + * the vertex dofs of the present element, whereas the second is the + * corresponding index of the other finite element. + */ + virtual + std::vector > + hp_vertex_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as hp_vertex_dof_indices(), except that the function treats degrees + * of freedom on lines. + */ + virtual + std::vector > + hp_line_dof_identities (const FiniteElement &fe_other) const; + + /** + * Same as hp_vertex_dof_indices(), except that the function treats degrees + * of freedom on quads. + */ + virtual + std::vector > + hp_quad_dof_identities (const FiniteElement &fe_other) const; + + /** + * Return whether this element dominates the one given as argument when they + * meet at a common face, whether it is the other way around, whether + * neither dominates, or if either could dominate. + * + * For a definition of domination, see FiniteElementBase::Domination and in + * particular the + * @ref hp_paper "hp paper". + */ + virtual + FiniteElementDomination::Domination + compare_for_face_domination (const FiniteElement &fe_other) const; + //@} + + + /** + * Return enrichment functions + */ + const std::vector *(const typename Triangulation::cell_iterator &) > > > + get_enrichments() const; + + /** + * Return the underlying FESystem object. + */ + const FESystem & + get_fe_system() const; + +protected: + + /** + * A class to hold internal data needed for evaluation of this FE at quadrature points. + */ + class InternalData : public FiniteElement::InternalDataBase + { + public: + /** + * For each Finite Element (base number) and each enrichment function (base_index) + * this struct will contain values, gradients and hessians of the enrichment functions. + */ + struct EnrichmentValues + { + std::vector values; + std::vector > gradients; + std::vector > hessians; + }; + + /** + * Constructor. Is used inside setup_data to wrap FESystem's internal + * data object. The former is called from get_data, get_subface_data and + * get_face_data which FE_Enriched has to implement. + * + * Since FESystem::get_data(), FESystem::get_face_data() and FESystem::get_subface_data() + * just create an object and return a pointer to it (i.e. they don't retain + * ownership), we store the cast result in a std::unique_ptr to indicate + * that InternalData owns the object. + */ + InternalData ( std::unique_ptr::InternalData> fesystem_data); + + /** + * Gives read-access to the pointer to a @p InternalData of the @p + * base_noth base element of FESystem's data. + */ + typename FiniteElement::InternalDataBase & + get_fe_data (const unsigned int base_no) const; + + /** + * Gives read-access to the pointer to an object into which the + * base_noth base element will write its output when calling + * FiniteElement::fill_fe_values() and similar functions. + */ + internal::FEValues::FiniteElementRelatedData & + get_fe_output_object (const unsigned int base_no) const; + + /** + * Aggregate FESystem's internal data. It is used every time + * we call FESystem's fill_fe_values() and alike. + */ + std::unique_ptr::InternalData> fesystem_data; + + /** + * For each FE used in enrichment (base number i) and each enrichment function + * (base multiplicity j), enrichment_values[i][j] will be used to store + * possibly requested values, gradients and hessians of enrichment function j. + * + * The variable is made mutable as InternalData's provided to fill_fe_values and alike + * are const. + * + * @note We do not want to store this information in the finite element object itself, + * because this would mean that (i) only one FEValues object could use a finite element object at a time, + * and (ii) that these objects could not be used in a multithreaded context. + */ + mutable std::vector > enrichment; + }; + + /** + * For each finite element @p i used in enrichment and each enrichment function + * @p j associated with it (essentially its multiplicity), + * @p base_no_mult_local_enriched_dofs[i][j] contains the associated local DoFs + * on the FE_Enriched finite element. + */ + std::vector > > base_no_mult_local_enriched_dofs; + + /** + * Enrichment functions. + * The size of the first vector is the same as the number of FiniteElement spaces used + * with enrichment. Whereas the size of the inner vector corresponds to the number + * of enrichment functions associated with a single FiniteElement. + */ + const std::vector *(const typename Triangulation::cell_iterator &) > > > enrichments; + + /** + * Auxiliary variable used to distinguish between the case when we do enrichment + * and when the class simply wraps another FiniteElement. + * + * This variable is initialized in the constructor by looping over a vector of + * enrichment elements and checking if all of them are FE_Nothing. If this is + * the case, then the value is set to false, otherwise it is + * true. + */ + const bool is_enriched; + + /** + * Auxiliary function called from get_data, get_face_data and + * get_subface_data. It take internal data of FESystem object in @p fes_data + * and the quadrature rule @p qudrature. + * + * This function essentially take the internal data from an instance of + * FESystem class and wraps it into our own InternalData class which + * additionally has objects to hold values/gradients/hessians of + * enrichment functions at each quadrature point depending on @p flags. + */ + template + typename FiniteElement::InternalDataBase * + setup_data (std::unique_ptr::InternalDataBase> fes_data, + const UpdateFlags flags, + const Quadrature &quadrature) const; + + /** + * Prepare internal data structures and fill in values independent of the + * cell. Returns a pointer to an object of which the caller of this function + * (FEValues) then has to assume ownership (which includes destruction when it is no + * more needed). + */ + virtual typename FiniteElement::InternalDataBase * + get_data (const UpdateFlags flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const; + + virtual typename FiniteElement::InternalDataBase * + get_face_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const; + + virtual typename FiniteElement::InternalDataBase * + get_subface_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData &output_data) const; + + virtual + void fill_fe_values (const typename Triangulation::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &output_data + ) const; + + virtual + void + fill_fe_face_values ( const typename Triangulation::cell_iterator &cell, + const unsigned int face_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &output_data + ) const; + + virtual + void + fill_fe_subface_values (const typename Triangulation< dim, spacedim >::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature &quadrature, + const Mapping &mapping, + const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + dealii::internal::FEValues::FiniteElementRelatedData &output_data + ) const; + +private: + /** + * This function sets up the index table for the system as well as @p + * restriction and @p prolongation matrices. + */ + void initialize (const std::vector*> &fes, + const std::vector &multiplicities); + + /** + * The underlying FESystem object. + */ + FESystem fe_system; + + /** + * After calling fill_fe_(face/subface_)values this function + * implements the chain rule to multiply stored shape values/gradient/hessians + * by those of enrichment function evaluated at quadrature points. + */ + template + void + multiply_by_enrichment (const Quadrature &quadrature, + const InternalData &fe_data, + const internal::FEValues::MappingRelatedData &mapping_data, + const typename Triangulation< dim, spacedim >::cell_iterator &cell, + internal::FEValues::FiniteElementRelatedData &output_data) const; +}; + +//} +DEAL_II_NAMESPACE_CLOSE + +#endif // CXX11 + +#endif // dealii__fe_enriched_h + diff --git a/include/deal.II/fe/fe_system.h b/include/deal.II/fe/fe_system.h index 1f87be95e3..89b017bd1e 100644 --- a/include/deal.II/fe/fe_system.h +++ b/include/deal.II/fe/fe_system.h @@ -28,6 +28,8 @@ DEAL_II_NAMESPACE_OPEN +template class FE_Enriched; + /** * This class provides an interface to group several elements together into @@ -959,7 +961,6 @@ private: unsigned int> > base_elements; - /** * Initialize the @p unit_support_points field of the FiniteElement class. * Called from the constructor. @@ -1071,6 +1072,8 @@ private: * Mutex for protecting initialization of restriction and embedding matrix. */ mutable Threads::Mutex mutex; + + friend class FE_Enriched; }; diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index 0e8f549797..f0a674c08c 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -44,6 +44,7 @@ SET(_src fe_raviart_thomas_nodal.cc fe_series.cc fe_system.cc + fe_enriched.cc fe_tools.cc fe_tools_interpolate.cc fe_trace.cc @@ -86,6 +87,7 @@ SET(_inst fe_raviart_thomas.inst.in fe_raviart_thomas_nodal.inst.in fe_system.inst.in + fe_enriched.inst.in fe_tools.inst.in fe_tools_interpolate.inst.in fe_trace.inst.in diff --git a/source/fe/fe_enriched.cc b/source/fe/fe_enriched.cc new file mode 100644 index 0000000000..7486d6bedd --- /dev/null +++ b/source/fe/fe_enriched.cc @@ -0,0 +1,926 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- +#include +#ifdef DEAL_II_WITH_CXX11 + +#include + +DEAL_II_NAMESPACE_OPEN + +namespace +{ + /** + * Auxiliary function to create multiplicity vector from input enrichment functions. + */ + template + std::vector + build_multiplicities(const std::vector > &functions ) + { + std::vector multiplicities; + multiplicities.push_back(1); // the first one is non-enriched FE + for (unsigned int i = 0; i < functions.size(); i++) + multiplicities.push_back(functions[i].size()); + + return multiplicities; + } + + + /** + * Auxiliary function to build FiniteElement's vector + */ + template + std::vector< const FiniteElement< dim, spacedim > * > + build_fes(const FiniteElement *fe_base, + const std::vector * > &fe_enriched) + { + std::vector< const FiniteElement< dim, spacedim > * > fes; + fes.push_back(fe_base); + for (unsigned int i = 0; i < fe_enriched.size(); i++) + fes.push_back(fe_enriched[i]); + + return fes; + } + + + /** + * Auxiliary function which check consistency of the input parameters. + * Returns true if everything is ok. + */ + template + bool + consistency_check (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) + { + AssertThrow(fes.size() > 0, + ExcMessage("FEs size should be >=1")); + AssertThrow(fes.size() == multiplicities.size(), + ExcMessage("FEs and multiplicities should have the same size")); + + AssertThrow (functions.size() == fes.size() - 1, + ExcDimensionMismatch(functions.size(), fes.size()-1)); + + AssertThrow(multiplicities[0] == 1, + ExcMessage("First multiplicity should be 1")); + + const unsigned int n_comp_base = fes[0]->n_components(); + + // start from fe=1 as 0th is always non-enriched FE. + for (unsigned int fe=1; fe < fes.size(); fe++) + { + const FE_Nothing *fe_nothing = dynamic_cast*>(fes[fe]); + if (fe_nothing) + AssertThrow (fe_nothing->is_dominating(), + ExcMessage("Only dominating FE_Nothing can be used in FE_Enriched")); + + AssertThrow (fes[fe]->n_components() == n_comp_base, + ExcMessage("All elements must have the same number of components")); + } + return true; + } + + + /** + * Auxiliary function which determines whether the FiniteElement will be enriched. + */ + template + bool + check_if_enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes) + { + // start from fe=1 as 0th is always non-enriched FE. + for (unsigned int fe=1; fe < fes.size(); fe++) + if (dynamic_cast*>(fes[fe]) == 0) + // this is not FE_Nothing => there will be enrichment + return true; + + return false; + } +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement &fe_base) + : + FE_Enriched(fe_base, + FE_Nothing(fe_base.n_components(),true), + NULL) +{ +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement &fe_base, + const FiniteElement &fe_enriched, + const Function *enrichment_function) + : + FE_Enriched + (&fe_base, { &fe_enriched }, +{ + { + [=] (const typename Triangulation::cell_iterator &) -> const Function * + { + return enrichment_function; + } + } +}) +{ +} + + +template +FE_Enriched::FE_Enriched (const FiniteElement *fe_base, + const std::vector * > &fe_enriched, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) + : + FE_Enriched (build_fes(fe_base,fe_enriched), + build_multiplicities(functions), + functions) +{ +} + + +template +FE_Enriched::FE_Enriched (const std::vector< const FiniteElement< dim, spacedim > * > &fes, + const std::vector< unsigned int > &multiplicities, + const std::vector *(const typename Triangulation::cell_iterator &) > > > &functions) + : + FiniteElement (FETools::Compositing::multiply_dof_numbers(fes,multiplicities,false), + FETools::Compositing::compute_restriction_is_additive_flags(fes,multiplicities), + FETools::Compositing::compute_nonzero_components(fes,multiplicities,false)), + enrichments(functions), + is_enriched(check_if_enriched(fes)), + fe_system(fes,multiplicities) +{ + // descriptive error are thrown within the function. + Assert(consistency_check(fes,multiplicities,functions), + ExcInternalError()); + + initialize(fes, multiplicities); + + // resize to be consistent with all FEs used to construct the FE_Enriched, + // even though we will never use the 0th element. + base_no_mult_local_enriched_dofs.resize(fes.size()); + for (unsigned int fe=1; fe < fes.size(); fe++) + base_no_mult_local_enriched_dofs[fe].resize(multiplicities[fe]); + + Assert (base_no_mult_local_enriched_dofs.size() == this->n_base_elements(), + ExcDimensionMismatch(base_no_mult_local_enriched_dofs.size(), + this->n_base_elements())); + + // build the map: (base_no, base_m) -> vector of local element DoFs + for (unsigned int system_index=0; system_indexdofs_per_cell; + ++system_index) + { + const unsigned int base_no = this->system_to_base_table[system_index].first.first; + if (base_no == 0) // 0th is always non-enriched FE + continue; + + const unsigned int base_m = this->system_to_base_table[system_index].first.second; + const unsigned int + base_index = this->system_to_base_table[system_index].second; + + Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(), + ExcMessage("Size mismatch for base_no_mult_local_enriched_dofs: " + "base_index = " + std::to_string(base_index) + + "; base_no = " + std::to_string(base_no) + + "; base_m = " + std::to_string(base_m) + + "; system_index = " + std::to_string(system_index))); + + Assert (base_m < base_no_mult_local_enriched_dofs[base_no].size(), + ExcDimensionMismatch(base_m, + base_no_mult_local_enriched_dofs[base_no].size())); + + base_no_mult_local_enriched_dofs[base_no][base_m].push_back(system_index); + } + + // make sure that local_enriched_dofs.size() is correct, that is equals to DoFs + // per cell of the corresponding FE. + for (unsigned int base_no = 1; base_no < base_no_mult_local_enriched_dofs.size(); base_no++) + { + for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++) + Assert ( base_no_mult_local_enriched_dofs[base_no][m].size() == fes[base_no]->dofs_per_cell, + ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no][m].size(), + fes[base_no]->dofs_per_cell)); + } +} + + +template +const std::vector *(const typename Triangulation::cell_iterator &) > > > +FE_Enriched::get_enrichments() const +{ + return enrichments; +} + + +template +double +FE_Enriched::shape_value(const unsigned int i, + const Point< dim > &p) const +{ + Assert(!is_enriched, ExcMessage("For enriched finite elements shape_value() can not be defined on the reference element.")); + return fe_system.shape_value(i,p); +} + + +template +FiniteElement * +FE_Enriched::clone() const +{ + std::vector< const FiniteElement< dim, spacedim > * > fes; + std::vector< unsigned int > multiplicities; + + for (unsigned int i=0; in_base_elements(); i++) + { + fes.push_back( & base_element(i) ); + multiplicities.push_back(this->element_multiplicity(i) ); + } + + return new FE_Enriched(fes, multiplicities, get_enrichments()); +} + + +template +UpdateFlags +FE_Enriched::requires_update_flags (const UpdateFlags flags) const +{ + UpdateFlags out = fe_system.requires_update_flags(flags); + + if (is_enriched) + { + // if we ask for values or gradients, then we would need quadrature points + if (flags & (update_values | update_gradients)) + out |= update_quadrature_points; + + // if need gradients, add update_values due to product rule + if (out & update_gradients) + out |= update_values; + } + + Assert (!(flags & update_3rd_derivatives), + ExcNotImplemented()); + + return out; +} + + +template +template +typename FiniteElement::InternalDataBase * +FE_Enriched::setup_data (std::unique_ptr::InternalDataBase> fes_data, + const UpdateFlags flags, + const Quadrature &quadrature) const +{ + Assert ((dynamic_cast::InternalData *> (fes_data.get()) != NULL), + ExcInternalError()); + typename FESystem::InternalData *data_fesystem = + static_cast::InternalData *> (fes_data.get()); + + // FESystem::InternalData will be aggregated (owned) by + // our InternalData. + fes_data.release(); + InternalData *data = new InternalData(std::unique_ptr::InternalData>(data_fesystem)); + + // copy update_each from FESystem data: + data->update_each = data_fesystem->update_each; + + // resize cache array according to requested flags + data->enrichment.resize(this->n_base_elements()); + + const unsigned int n_q_points = quadrature.size(); + + for (unsigned int base=0; base < this->n_base_elements(); ++base) + { + data->enrichment[base].resize(this->element_multiplicity(base)); + for (unsigned int m = 0; m < this->element_multiplicity(base); ++m) + { + if (flags & update_values) + data->enrichment[base][m].values.resize(n_q_points); + + if (flags & update_gradients) + data->enrichment[base][m].gradients.resize(n_q_points); + + if (flags & update_hessians) + data->enrichment[base][m].hessians.resize(n_q_points); + } + } + + return data; +} + + +template +typename FiniteElement::InternalDataBase * +FE_Enriched::get_face_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_face_data(update_flags,mapping,quadrature,output_data)), + update_flags, + quadrature); +} + + +template +typename FiniteElement::InternalDataBase * +FE_Enriched::get_subface_data (const UpdateFlags update_flags, + const Mapping &mapping, + const Quadrature &quadrature, + dealii::internal::FEValues::FiniteElementRelatedData &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_subface_data(update_flags,mapping,quadrature,output_data)), + update_flags, + quadrature); +} + + +template +typename FiniteElement::InternalDataBase * +FE_Enriched::get_data (const UpdateFlags flags, + const Mapping &mapping, + const Quadrature &quadrature, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data) const +{ + return setup_data(std::unique_ptr::InternalDataBase>(fe_system.get_data(flags,mapping,quadrature,output_data)), + flags, + quadrature); +} + + +template +void FE_Enriched::initialize (const std::vector*> &fes, + const std::vector &multiplicities) +{ + Assert (fes.size() == multiplicities.size(), + ExcDimensionMismatch (fes.size(), multiplicities.size()) ); + + // Note that we need to skip every fe with multiplicity 0 in the following block of code + this->base_to_block_indices.reinit(0, 0); + + for (unsigned int i=0; i0) + this->base_to_block_indices.push_back( multiplicities[i] ); + + { + // If the system is not primitive, these have not been initialized by + // FiniteElement + this->system_to_component_table.resize(this->dofs_per_cell); + this->face_system_to_component_table.resize(this->dofs_per_face); + + FETools::Compositing::build_cell_tables(this->system_to_base_table, + this->system_to_component_table, + this->component_to_base_table, + *this, + false); + + FETools::Compositing::build_face_tables(this->face_system_to_base_table, + this->face_system_to_component_table, + *this, + false); + } + + // restriction and prolongation matrices are built on demand + + // now set up the interface constraints for h-refinement. + // take them from fe_system: + this->interface_constraints = fe_system.interface_constraints; + + // if we just wrap another FE (i.e. use FE_Nothing as a second FE) + // then it makes sense to have support points. + // However, functions like interpolate_boundary_values() need all FEs inside + // FECollection to be able to provide support points irrespectively whether + // this FE sits on the boundary or not. Thus for moment just copy support + // points from fe system: + { + this->unit_support_points = fe_system.unit_support_points; + this->unit_face_support_points = fe_system.unit_face_support_points; + } + + // take adjust_quad_dof_index_for_face_orientation_table from FESystem: + { + this->adjust_line_dof_index_for_line_orientation_table = fe_system.adjust_line_dof_index_for_line_orientation_table; + } +} + + +template +std::string +FE_Enriched::get_name () const +{ + std::ostringstream namebuf; + + namebuf << "FE_Enriched<" + << Utilities::dim_string(dim,spacedim) + << ">["; + for (unsigned int i=0; i< this->n_base_elements(); ++i) + { + namebuf << base_element(i).get_name(); + if (this->element_multiplicity(i) != 1) + namebuf << '^' << this->element_multiplicity(i); + if (i != this->n_base_elements()-1) + namebuf << '-'; + } + namebuf << ']'; + + return namebuf.str(); +} + + +template +const FiniteElement & +FE_Enriched::base_element (const unsigned int index) const +{ + return fe_system.base_element(index); +} + + +template +void +FE_Enriched::fill_fe_values (const typename Triangulation< dim, spacedim >::cell_iterator &cell, + const CellSimilarity::Similarity cell_similarity, + const Quadrature< dim > &quadrature, + const Mapping< dim, spacedim > &mapping, + const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data + ) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + // call FESystem's method to fill everything without enrichment function + fe_system.fill_fe_values(cell, + cell_similarity, + quadrature, + mapping, + mapping_internal, + mapping_data, + *fe_data.fesystem_data, + output_data); + + if (is_enriched) + multiply_by_enrichment(quadrature, + fe_data, + mapping_data, + cell, + output_data); +} + + +template +void +FE_Enriched::fill_fe_face_values +(const typename Triangulation< dim, spacedim >::cell_iterator &cell, + const unsigned int face_no, + const Quadrature< dim-1 > &quadrature, + const Mapping< dim, spacedim > &mapping, + const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data +) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + // call FESystem's method to fill everything without enrichment function + fe_system.fill_fe_face_values(cell, + face_no, + quadrature, + mapping, + mapping_internal, + mapping_data, + *fe_data.fesystem_data, + output_data); + + if (is_enriched) + multiply_by_enrichment(quadrature, + fe_data, + mapping_data, + cell, + output_data); +} + + +template +void +FE_Enriched::fill_fe_subface_values +(const typename Triangulation< dim, spacedim >::cell_iterator &cell, + const unsigned int face_no, + const unsigned int sub_no, + const Quadrature< dim-1 > &quadrature, + const Mapping< dim, spacedim > &mapping, + const typename Mapping< dim, spacedim >::InternalDataBase &mapping_internal, + const dealii::internal::FEValues::MappingRelatedData< dim, spacedim > &mapping_data, + const typename FiniteElement::InternalDataBase &fe_internal, + internal::FEValues::FiniteElementRelatedData< dim, spacedim > &output_data +) const +{ + Assert (dynamic_cast (&fe_internal) != NULL, + ExcInternalError()); + const InternalData &fe_data = static_cast (fe_internal); + + // call FESystem's method to fill everything without enrichment function + fe_system.fill_fe_subface_values(cell, + face_no, + sub_no, + quadrature, + mapping, + mapping_internal, + mapping_data, + *fe_data.fesystem_data, + output_data); + + if (is_enriched) + multiply_by_enrichment(quadrature, + fe_data, + mapping_data, + cell, + output_data); + +} + + +template +template +void +FE_Enriched::multiply_by_enrichment +(const Quadrature &quadrature, + const InternalData &fe_data, + const internal::FEValues::MappingRelatedData &mapping_data, + const typename Triangulation< dim, spacedim >::cell_iterator &cell, + internal::FEValues::FiniteElementRelatedData &output_data) const +{ + // mapping_data will contain quadrature points on the real element. + // fe_internal is needed to get update flags + // finally, output_data should store all the results we need. + + // Either dim_1==dim + // (fill_fe_values) or dim_1==dim-1 + // (fill_fe_(sub)face_values) + Assert(dim_1==dim || dim_1==dim-1, ExcInternalError()); + const UpdateFlags flags = fe_data.update_each; + + const unsigned int n_q_points = quadrature.size(); + + // First, populate output_data object (that shall hold everything requested such + // as shape value, gradients, hessians, etc) from each base element. + // That is almost identical to FESystem::compute_fill_one_base(), + // the difference being that we do it irrespectively of cell_similarity + // and use base_fe_data.update_flags + + // TODO: do we need it only for dim_1 == dim (i.e. fill_fe_values)? + if (dim_1 == dim) + for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++) + { + const FiniteElement & + base_fe = base_element(base_no); + typename FiniteElement::InternalDataBase & + base_fe_data = fe_data.get_fe_data(base_no); + internal::FEValues::FiniteElementRelatedData & + base_data = fe_data.get_fe_output_object(base_no); + + const UpdateFlags base_flags = base_fe_data.update_each; + + for (unsigned int system_index=0; system_indexdofs_per_cell; + ++system_index) + if (this->system_to_base_table[system_index].first.first == base_no) + { + const unsigned int + base_index = this->system_to_base_table[system_index].second; + Assert (base_indexn_nonzero_components(i); + unsigned int in_index = 0; + for (unsigned int i=0; in_nonzero_components(system_index) == + base_fe.n_nonzero_components(base_index), + ExcInternalError()); + for (unsigned int s=0; sn_nonzero_components(system_index); ++s) + { + if (base_flags & update_values) + for (unsigned int q=0; qn_base_elements(); base_no++) + { + Assert (base_no_mult_local_enriched_dofs[base_no].size() == fe_data.enrichment[base_no].size(), + ExcDimensionMismatch(base_no_mult_local_enriched_dofs[base_no].size(), + fe_data.enrichment[base_no].size())); + for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++) + { + Assert (enrichments[base_no-1][m](cell) !=NULL, + ExcMessage("The pointer to the enrichment function is NULL")); + + Assert (enrichments[base_no-1][m](cell)->n_components == 1, + ExcMessage("Only scalar-valued enrichment functions are allowed")); + + if (flags & update_hessians) + { + Assert (fe_data.enrichment[base_no][m].hessians.size() == n_q_points, + ExcDimensionMismatch(fe_data.enrichment[base_no][m].hessians.size(), + n_q_points)); + for (unsigned int q=0; qhessian (mapping_data.quadrature_points[q]); + } + + if (flags & update_gradients) + { + Assert (fe_data.enrichment[base_no][m].gradients.size() == n_q_points, + ExcDimensionMismatch(fe_data.enrichment[base_no][m].gradients.size(), + n_q_points)); + for (unsigned int q=0; qgradient(mapping_data.quadrature_points[q]); + } + + if (flags & update_values) + { + Assert (fe_data.enrichment[base_no][m].values.size() == n_q_points, + ExcDimensionMismatch(fe_data.enrichment[base_no][m].values.size(), + n_q_points)); + for (unsigned int q=0; qvalue (mapping_data.quadrature_points[q]); + } + } + } + + // Finally, update the standard data stored in output_data + // by expanding the product rule for enrichment function. + // note that the order if important, namely + // output_data.shape_XYZ contains values of standard FEM and we overwrite + // it with the updated one in the following order: hessians -> gradients -> values + if (flags & update_hessians) + { + for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++) + { + for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++) + for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++) + { + const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i]; + for (unsigned int q=0; q grad_grad = outer_product(output_data.shape_gradients[enriched_dof][q], + fe_data.enrichment[base_no][m].gradients[q]); + const Tensor<2,spacedim,double> sym_grad_grad = symmetrize (grad_grad) * 2.0; // symmetrize does [s+s^T]/2 + + output_data.shape_hessians [enriched_dof][q]*= fe_data.enrichment[base_no][m].values[q]; + output_data.shape_hessians [enriched_dof][q]+= + sym_grad_grad + + output_data.shape_values [enriched_dof][q] * fe_data.enrichment[base_no][m].hessians[q]; + } + } + } + } + + if (flags & update_gradients) + for (unsigned int base_no = 1; base_no < this->n_base_elements(); base_no++) + { + for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++) + for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++) + { + const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i]; + for (unsigned int q=0; qn_base_elements(); base_no++) + { + for (unsigned int m=0; m < base_no_mult_local_enriched_dofs[base_no].size(); m++) + for (unsigned int i=0; i < base_no_mult_local_enriched_dofs[base_no][m].size(); i++) + { + const unsigned int enriched_dof = base_no_mult_local_enriched_dofs[base_no][m][i]; + for (unsigned int q=0; q +const FESystem & +FE_Enriched:: +get_fe_system() const +{ + return fe_system; +} + + +template +bool +FE_Enriched:: +hp_constraints_are_implemented () const +{ + return true; +} + + +template +void +FE_Enriched:: +get_face_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&source)) + { + fe_system.get_face_interpolation_matrix(fe_enr_other->get_fe_system(), + matrix); + } + else + { + typedef FiniteElement FEL; + AssertThrow(false,typename FEL::ExcInterpolationNotImplemented()); + } +} + + +template +void +FE_Enriched:: +get_subface_interpolation_matrix (const FiniteElement &source, + const unsigned int subface, + FullMatrix &matrix) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&source)) + { + fe_system.get_subface_interpolation_matrix(fe_enr_other->get_fe_system(), + subface, + matrix); + } + else + { + typedef FiniteElement FEL; + AssertThrow(false,typename FEL::ExcInterpolationNotImplemented()); + } +} + + +template +std::vector > +FE_Enriched:: +hp_vertex_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_vertex_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +std::vector > +FE_Enriched:: +hp_line_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_line_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +std::vector > +FE_Enriched:: +hp_quad_dof_identities (const FiniteElement &fe_other) const +{ + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.hp_quad_dof_identities(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return std::vector >(); + } +} + + +template +FiniteElementDomination::Domination +FE_Enriched:: +compare_for_face_domination (const FiniteElement &fe_other) const +{ + // need to decide which element constrain another. + // for example Q(2) dominate Q(4) and thus some DoFs of Q(4) will be constrained. + // If we have Q(2) and Q(4)+POU, then it's clear that Q(2) dominates, + // namely our DoFs will be constrained to make field continuous. + // However, we need to check for situations like Q(4) vs Q(2)+POU. + // In that case the domination for the underlying FEs should be the other way, + // but this implies that we can't constrain POU dofs to make the field continuous. + // In that case, through an error + + // if it's also enriched, do domination based on each one's FESystem + if (const FE_Enriched *fe_enr_other + = dynamic_cast*>(&fe_other)) + { + return fe_system.compare_for_face_domination(fe_enr_other->get_fe_system()); + } + else + { + Assert (false, ExcNotImplemented()); + return FiniteElementDomination::neither_element_dominates; + } +} + +/* ----------------------- FESystem::InternalData ------------------- */ + + +template +FE_Enriched::InternalData::InternalData(std::unique_ptr::InternalData> fesystem_data) + : + fesystem_data(std::move(fesystem_data)) +{} + + +template +typename FiniteElement::InternalDataBase & +FE_Enriched:: +InternalData::get_fe_data (const unsigned int base_no) const +{ + return fesystem_data->get_fe_data(base_no); +} + + +template +internal::FEValues::FiniteElementRelatedData & +FE_Enriched:: +InternalData::get_fe_output_object (const unsigned int base_no) const +{ + return fesystem_data->get_fe_output_object(base_no); +} + + +// explicit instantiations +#include "fe_enriched.inst" + +DEAL_II_NAMESPACE_CLOSE + +#endif // CXX11 diff --git a/source/fe/fe_enriched.inst.in b/source/fe/fe_enriched.inst.in new file mode 100644 index 0000000000..853f36a622 --- /dev/null +++ b/source/fe/fe_enriched.inst.in @@ -0,0 +1,23 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) +{ +#if deal_II_dimension <= deal_II_space_dimension + template class FE_Enriched; +#endif +} + diff --git a/tests/fe/fe_enriched_01.cc b/tests/fe/fe_enriched_01.cc new file mode 100644 index 0000000000..9413d6fe4a --- /dev/null +++ b/tests/fe/fe_enriched_01.cc @@ -0,0 +1,127 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 basic output of the class (name, n_blocks, n_compoennts, n_dofs_per_cell) + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + +template +void test_base() +{ + deallog << "Test basic functions:"< function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + deallog << fe.get_name()<(); + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_01.with_cxx11=on.output b/tests/fe/fe_enriched_01.with_cxx11=on.output new file mode 100644 index 0000000000..d5e7fca629 --- /dev/null +++ b/tests/fe/fe_enriched_01.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::Test basic functions: +DEAL::FE_Enriched<3>[FE_Q<3>(1)-FE_Q<3>(1)] +DEAL::2 +DEAL::1 +DEAL::16 diff --git a/tests/fe/fe_enriched_02.cc b/tests/fe/fe_enriched_02.cc new file mode 100644 index 0000000000..fb6ac2069e --- /dev/null +++ b/tests/fe/fe_enriched_02.cc @@ -0,0 +1,163 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check FEValues::shape_value() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + +template +void test1() +{ + deallog << "FEValues.shape_value()"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(1); + FEValues fe_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + const unsigned int n_q_points = quadrature.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const std::vector > &q_points = fe_values.get_quadrature_points(); + + for (unsigned int i=0; i(); + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_02.with_cxx11=on.output b/tests/fe/fe_enriched_02.with_cxx11=on.output new file mode 100644 index 0000000000..00053ee909 --- /dev/null +++ b/tests/fe/fe_enriched_02.with_cxx11=on.output @@ -0,0 +1,20 @@ + +DEAL::FEValues.shape_value() +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL::dof=0 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=1 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=2 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=3 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=4 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=5 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=6 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=7 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=8 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=9 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=10 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=11 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=12 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=13 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL::dof=14 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.1250 +DEAL::dof=15 qp=0.5000 0.5000 0.5000 f(qp)=0.4206 N(qp)=0.0526 +DEAL:: diff --git a/tests/fe/fe_enriched_03.cc b/tests/fe/fe_enriched_03.cc new file mode 100644 index 0000000000..73ec999994 --- /dev/null +++ b/tests/fe/fe_enriched_03.cc @@ -0,0 +1,157 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// Check FEFaceValues::shape_value() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test2() +{ + deallog << "FEFaceValues"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(1); + FEFaceValues fe_face_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; ++face) + { + fe_face_values.reinit (cell,face); + const unsigned int n_q_points = quadrature.size(); + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const std::vector > &q_points = fe_face_values.get_quadrature_points(); + + for (unsigned int i=0; i(); + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_03.with_cxx11=on.output b/tests/fe/fe_enriched_03.with_cxx11=on.output new file mode 100644 index 0000000000..ef9efb93fa --- /dev/null +++ b/tests/fe/fe_enriched_03.with_cxx11=on.output @@ -0,0 +1,99 @@ + +DEAL::FEFaceValues +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL::dof=0 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=3 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=4 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=5 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=6 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=7 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=8 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=9 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=10 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=11 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=12 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=13 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=14 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.0000 0.5000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=3 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=4 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=5 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=6 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=7 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=8 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=9 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=10 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=11 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=12 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=13 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=14 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=1.0000 0.5000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=0 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=3 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=4 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=5 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=6 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=7 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=8 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=9 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=10 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=11 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=12 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=13 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=14 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.5000 0.0000 0.5000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=3 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=4 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=5 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=6 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=7 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=8 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=9 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=10 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=11 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=12 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=13 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=14 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=0.5000 1.0000 0.5000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=0 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=1 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=2 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=3 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=4 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=5 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=6 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.2500 +DEAL::dof=7 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0.1233 +DEAL::dof=8 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=9 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=10 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=11 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=12 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=13 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=14 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=15 qp=0.5000 0.5000 0.0000 f(qp)=0.4931 N(qp)=0 +DEAL::dof=0 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=1 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=2 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=3 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=4 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=5 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=6 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=7 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0 +DEAL::dof=8 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=9 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=10 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=11 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=12 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=13 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 +DEAL::dof=14 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.2500 +DEAL::dof=15 qp=0.5000 0.5000 1.0000 f(qp)=0.2938 N(qp)=0.0735 diff --git a/tests/fe/fe_enriched_04.cc b/tests/fe/fe_enriched_04.cc new file mode 100644 index 0000000000..2d3d921658 --- /dev/null +++ b/tests/fe/fe_enriched_04.cc @@ -0,0 +1,244 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 FEValues::get_function_values() + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test3() +{ + deallog << "FEValues.get_function_values()"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + dof_handler.distribute_dofs(fe); + + QGauss quadrature(2); + FEValues fe_values(fe, quadrature, + update_values | update_gradients | + update_JxW_values); + + Vector solution_fe(dof_handler.n_dofs()), + solution_pou(dof_handler.n_dofs()), + solution(dof_handler.n_dofs()); + solution_fe = 0; + solution_pou = 0; + + // groups are one after another + // thus set to 1.0 the same shape function in the underlying FE + solution_fe[0] = 2.0; + solution_pou[1] = 2.0; + solution = solution_fe; + solution+= solution_pou; + + const unsigned int n_q_points = quadrature.size(); + + std::vector solution_values_fe(n_q_points), + solution_values_pou(n_q_points), + solution_values(n_q_points); + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + const std::vector > &q_points = fe_values.get_quadrature_points(); + fe_values.get_function_values (solution_fe, + solution_values_fe); + fe_values.get_function_values (solution_pou, + solution_values_pou); + fe_values.get_function_values (solution, + solution_values); + + for (unsigned int q_point=0; q_point data_out; + data_out.attach_dof_handler (dof_handler); + + 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_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + try + { + test3<3>(); +#ifdef DATA_OUT_FE_ENRICHED + plot_shape_function<1>(); + plot_shape_function<2>(); + plot_shape_function<3>(); +#endif + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_04.with_cxx11=on.output b/tests/fe/fe_enriched_04.with_cxx11=on.output new file mode 100644 index 0000000000..785ca6ae19 --- /dev/null +++ b/tests/fe/fe_enriched_04.with_cxx11=on.output @@ -0,0 +1,11 @@ + +DEAL::FEValues.get_function_values() +DEAL::for same underlying FEs: f(qp) * N_{fe}(qp) == N_{pou}(qp) +DEAL:: qp=0.2113 0.2113 0.2113 f(qp)=0.6935 U_fe(qp)=0.9811 U_pou(qp)=0.6804 U_[fe+pou](qp)=1.6615 +DEAL:: qp=0.7887 0.2113 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.2113 0.7887 0.2113 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.7887 0.7887 0.2113 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.2113 0.2113 0.7887 f(qp)=0.4302 U_fe(qp)=0.2629 U_pou(qp)=0.1131 U_[fe+pou](qp)=0.3760 +DEAL:: qp=0.7887 0.2113 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.2113 0.7887 0.7887 f(qp)=0.3214 U_fe(qp)=0.0704 U_pou(qp)=0.0226 U_[fe+pou](qp)=0.0931 +DEAL:: qp=0.7887 0.7887 0.7887 f(qp)=0.2551 U_fe(qp)=0.0189 U_pou(qp)=0.0048 U_[fe+pou](qp)=0.0237 diff --git a/tests/fe/fe_enriched_05.cc b/tests/fe/fe_enriched_05.cc new file mode 100644 index 0000000000..6a2a63ce30 --- /dev/null +++ b/tests/fe/fe_enriched_05.cc @@ -0,0 +1,186 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 h-refinement with DoFHandler and print constraints. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test4() +{ + deallog << "h-refinement:"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(2), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + + // one global refinement and 1 refinement of 1 active cell (local) + triangulation.refine_global(); + triangulation.begin_active()->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + + dof_handler.distribute_dofs(fe); + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + constraints.print(deallog.get_file_stream()); + +#ifdef DATA_OUT_FE_ENRICHED + std::vector> shape_functions; + std::vector names; + for (unsigned int s=0; s < dof_handler.n_dofs(); s++) + if (!constraints.is_constrained(s)) + { + names.push_back(std::string("N_") + + dealii::Utilities::int_to_string(s)); + + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[s] = 1.0; + + // make continuous: + constraints.distribute(shape_function); + shape_functions.push_back(shape_function); + } + + DataOut data_out; + data_out.attach_dof_handler (dof_handler); + + 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 = "h-shape_functions_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test4<2>(); + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_05.with_cxx11=on.output b/tests/fe/fe_enriched_05.with_cxx11=on.output new file mode 100644 index 0000000000..764434aac4 --- /dev/null +++ b/tests/fe/fe_enriched_05.with_cxx11=on.output @@ -0,0 +1,20 @@ + +DEAL::h-refinement: + 42 8: 1.00000 + 43 1: 0.500000 + 43 5: 0.500000 + 44 0: 0.375000 + 44 4: -0.125000 + 44 8: 0.750000 + 48 21: 1.00000 + 49 5: 0.500000 + 49 14: 0.500000 + 52 4: -0.125000 + 52 13: 0.375000 + 52 21: 0.750000 + 54 0: -0.125000 + 54 4: 0.375000 + 54 8: 0.750000 + 55 4: 0.375000 + 55 13: -0.125000 + 55 21: 0.750000 diff --git a/tests/fe/fe_enriched_06.cc b/tests/fe/fe_enriched_06.cc new file mode 100644 index 0000000000..01b33e104e --- /dev/null +++ b/tests/fe/fe_enriched_06.cc @@ -0,0 +1,223 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// no refinement, just 4 cells. +// Make sure fe_values and shape functions are continuous + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test5() +{ + deallog << "4 cells:"< triangulation; + DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + FE_Enriched fe(FE_Q(1), + FE_Q(1), + &function); + + GridGenerator::hyper_cube(triangulation); + + triangulation.refine_global(); + + dof_handler.distribute_dofs(fe); + + std::vector> shape_functions; + std::vector names; + for (unsigned int s=0; s < dof_handler.n_dofs(); s++) + { + names.push_back(std::string("N_") + + dealii::Utilities::int_to_string(s)); + + Vector shape_function; + shape_function.reinit(dof_handler.n_dofs()); + shape_function[s] = 1.0; + + shape_functions.push_back(shape_function); + } + + // output 11th: + { + const unsigned int global_dof = 11; + const Vector &solution = shape_functions[global_dof]; + QTrapez quadrature; + FEValues fe_values(fe, quadrature, + update_values); + + const unsigned int n_q_points = quadrature.size(); + std::vector solution_values(n_q_points); + + const unsigned int dofs_per_cell = fe.dofs_per_cell; + std::vector local_dof_indices; + local_dof_indices.resize (dofs_per_cell); + + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + { + fe_values.reinit (cell); + + // find out which + unsigned int local_dof = 0; + cell->get_dof_indices (local_dof_indices); + for ( ; local_dof < dofs_per_cell; local_dof++) + if (local_dof_indices[local_dof]==global_dof) + break; + + const std::vector > &q_points = fe_values.get_quadrature_points(); + fe_values.get_function_values (solution, + solution_values); + + deallog <<" cell="<center()<(); + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_06.with_cxx11=on.output b/tests/fe/fe_enriched_06.with_cxx11=on.output new file mode 100644 index 0000000000..a45dcea52b --- /dev/null +++ b/tests/fe/fe_enriched_06.with_cxx11=on.output @@ -0,0 +1,22 @@ + +DEAL::4 cells: +DEAL:: cell=0.2500 0.2500 +DEAL:: qp=0.0000 0.0000 f(qp)=1.0000 U(qp)=0 +DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0 +DEAL:: cell=0.7500 0.2500 +DEAL:: qp=0.5000 0.0000 f(qp)=0.6065 N(7,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.0000 f(qp)=0.3679 N(7,qp)=0 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(7,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(7,qp)=0.3269 U(qp)=0.3269 +DEAL:: cell=0.2500 0.7500 +DEAL:: qp=0.0000 0.5000 f(qp)=0.6065 U(qp)=0 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 U(qp)=0 +DEAL:: qp=0.0000 1.0000 f(qp)=0.3679 U(qp)=0 +DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 U(qp)=0 +DEAL:: cell=0.7500 0.7500 +DEAL:: qp=0.5000 0.5000 f(qp)=0.4931 N(3,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 0.5000 f(qp)=0.3269 N(3,qp)=0.3269 U(qp)=0.3269 +DEAL:: qp=0.5000 1.0000 f(qp)=0.3269 N(3,qp)=0 U(qp)=0 +DEAL:: qp=1.0000 1.0000 f(qp)=0.2431 N(3,qp)=0 U(qp)=0 diff --git a/tests/fe/fe_enriched_07.cc b/tests/fe/fe_enriched_07.cc new file mode 100644 index 0000000000..d0512022fe --- /dev/null +++ b/tests/fe/fe_enriched_07.cc @@ -0,0 +1,237 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// hp-refinement. Output constraints. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test6(const bool do_href, + const unsigned int p_feq = 1, + const unsigned int p_feen = 2) +{ + deallog << "hp: "< triangulation; + hp::DoFHandler dof_handler(triangulation); + + EnrichmentFunction function; + + hp::FECollection fe_collection; + fe_collection.push_back(FE_Enriched(FE_Q(p_feq))); + fe_collection.push_back(FE_Enriched(FE_Q(p_feen), + FE_Q(1), + &function)); + // push back to be able to resolve hp constrains: + fe_collection.push_back(FE_Enriched(FE_Q(p_feen))); + + GridGenerator::hyper_cube(triangulation); + + // one global refinement and 1 refinement of 1 active cell + // a sketch with material ids: + // --------------- + // | | | + // | 1 | 0 | + // |------|------- + // |1 |1 | | + // |------| 0 | + // |1 |1 | | + // |------|------| + triangulation.refine_global(); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell->set_active_fe_index(1);//POU + cell++; + cell->set_active_fe_index(1);//POU + } + dof_handler.distribute_dofs(fe_collection); + + if (do_href) + { + triangulation.begin_active()->set_refine_flag (); + triangulation.execute_coarsening_and_refinement (); + dof_handler.distribute_dofs(fe_collection); + } + + deallog << "n_cells: "<> 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(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_ref="+std::to_string(do_href) + +"_p_feq="+std::to_string(p_feq) + +"_p_feenr="+std::to_string(p_feen) + +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test6<2>(false,1,2);// 1 vs 2+1 + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_07.with_cxx11=on.output b/tests/fe/fe_enriched_07.with_cxx11=on.output new file mode 100644 index 0000000000..e4f68ec21f --- /dev/null +++ b/tests/fe/fe_enriched_07.with_cxx11=on.output @@ -0,0 +1,10 @@ + +DEAL::hp: 0 1 2 +DEAL::n_cells: 4 + 4 = 0 + 5 = 0 + 9 18: 0.500000 + 9 19: 0.500000 + 13 = 0 + 16 19: 0.500000 + 16 22: 0.500000 diff --git a/tests/fe/fe_enriched_08.cc b/tests/fe/fe_enriched_08.cc new file mode 100644 index 0000000000..c59363d6bc --- /dev/null +++ b/tests/fe/fe_enriched_08.cc @@ -0,0 +1,219 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// check constraints on hp dofhandler with 2 neighboring cells where one +// cell is enriched and nother is not. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// argument for build_patches() +const unsigned int patches = 10; + +using namespace dealii; + +// uncomment when debugging +// #define DATA_OUT_FE_ENRICHED + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } +}; + + +template +void test2cells(const unsigned int p_feq=2, + const unsigned int p_feen=1) +{ + deallog << "2cells: "< triangulation; + { + Triangulation triangulationL; + Triangulation triangulationR; + GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain + GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain + Point shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + } + + hp::DoFHandler dof_handler(triangulation); + EnrichmentFunction function; + + hp::FECollection fe_collection; + fe_collection.push_back(FE_Enriched(FE_Q(p_feq))); + fe_collection.push_back(FE_Enriched(FE_Q(p_feen), + FE_Q(1), + &function)); + // push back to be able to resolve hp constrains: + fe_collection.push_back(FE_Enriched(FE_Q(p_feen))); + + dof_handler.begin_active()->set_active_fe_index(1);//POU + + + dof_handler.distribute_dofs(fe_collection); + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + constraints.print(deallog.get_file_stream()); + +#ifdef DATA_OUT_FE_ENRICHED + // 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(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 = "2cell-shape_functions_p_feq="+std::to_string(p_feq) + +"_p_feenr="+std::to_string(p_feen) + +"_"+dealii::Utilities::int_to_string(dim)+"D.vtu"; + std::ofstream output (filename.c_str ()); + data_out.write_vtu (output); +#endif + + dof_handler.clear(); +} + + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test2cells<2>(1,2); // 1 vs 2+1 + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_08.with_cxx11=on.output b/tests/fe/fe_enriched_08.with_cxx11=on.output new file mode 100644 index 0000000000..3a075cc441 --- /dev/null +++ b/tests/fe/fe_enriched_08.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::2cells: 2 1 2 + 2 = 0 + 5 = 0 + 7 11: 0.500000 + 7 13: 0.500000 diff --git a/tests/fe/fe_enriched_08a.cc b/tests/fe/fe_enriched_08a.cc new file mode 100644 index 0000000000..23fd6dbdae --- /dev/null +++ b/tests/fe/fe_enriched_08a.cc @@ -0,0 +1,133 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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. +// +// --------------------------------------------------------------------- + + +// Same as fe_enriched_08, but checks the constraints from FESystem directly. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +using namespace dealii; + +template +void test2cellsFESystem(const unsigned int p_feq=2, + const unsigned int p_feen=1) +{ + deallog << "2cells: "< triangulation; + { + Triangulation triangulationL; + Triangulation triangulationR; + GridGenerator::hyper_cube (triangulationL, -1,0); //create a square [-1,0]^d domain + GridGenerator::hyper_cube (triangulationR, -1,0); //create a square [-1,0]^d domain + Point shift_vector; + shift_vector[0] = 1.0; + GridTools::shift(shift_vector,triangulationR); + GridGenerator::merge_triangulations (triangulationL, triangulationR, triangulation); + } + + hp::DoFHandler dof_handler(triangulation); + + hp::FECollection fe_collection; + fe_collection.push_back(FESystem(FE_Q(p_feq),1, + FE_Nothing(),1)); + fe_collection.push_back(FESystem(FE_Q(p_feen),1, + FE_Q(1),1)); + + // push back to be able to resolve hp constrains: + fe_collection.push_back(FESystem(FE_Q(p_feen),1, + FE_Nothing(),1)); + + dof_handler.begin_active()->set_active_fe_index(1); + + dof_handler.distribute_dofs(fe_collection); + + ConstraintMatrix constraints; + constraints.clear(); + dealii::DoFTools::make_hanging_node_constraints (dof_handler, constraints); + constraints.close (); + + constraints.print(deallog.get_file_stream()); + + dof_handler.clear(); +} + +int main (int argc,char **argv) +{ + std::ofstream logfile ("output"); + deallog << std::setprecision(4); + deallog << std::fixed; + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + + try + { + test2cellsFESystem<2>(1,2); // 1+0 vs 2+1 + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_08a.with_cxx11=on.output b/tests/fe/fe_enriched_08a.with_cxx11=on.output new file mode 100644 index 0000000000..3a075cc441 --- /dev/null +++ b/tests/fe/fe_enriched_08a.with_cxx11=on.output @@ -0,0 +1,6 @@ + +DEAL::2cells: 2 1 2 + 2 = 0 + 5 = 0 + 7 11: 0.500000 + 7 13: 0.500000 diff --git a/tests/fe/fe_enriched_compare_to_fe_system.cc b/tests/fe/fe_enriched_compare_to_fe_system.cc new file mode 100644 index 0000000000..f19cb7c120 --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system.cc @@ -0,0 +1,320 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 FE_Enriched class by comparing with FE_System for a single enrichment +// function. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// create a mesh without hanging nodes, +// move nodes randomly +// and compare FEEnriched to FESystem with explicit treatment of +// the product rule for +// valus, gradients, hessians on +// elements and faces. +// The comparison is straight forward because local dofs are enumerated +// in the same way for FE_System and FEEnriched. +const unsigned int patches = 10; + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction() + : Function(1) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + return std::exp(-point.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> res = point; + Assert (point.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + res *= -value(point)/point.norm(); + return res; + } + + virtual SymmetricTensor<2,dim> hessian (const Point &p, + const unsigned int component=0) const + { + Tensor<1,dim> dir = p; + const double r = dir.norm(); + Assert (r>0.0, ExcMessage("r is not positive")); + dir/=r; + SymmetricTensor<2,dim> dir_x_dir; + for (unsigned int i=0; i()/r ); + } + +}; + +/** + * + * @param p point + * @param func function + * @param v_e value of enriched FE + * @param g_e gradient of enriched FE + * @param h_e hessian of enriched FE + * @param v_s0 value of 0th component of FE_System + * @param g_s0 gradient of 0th component of FE_System + * @param h_s0 hessian of 0th component of FE_System + * @param v_s1 value of 1st component of FE_System + * @param g_s1 gradient of 1st component of FE_System + * @param h_s1 hessian of 1st component of FE_System + */ +template +void check_consistency(const Point &p, + const Function &func, + const double &v_e, + const Tensor<1,dim> &g_e, + const Tensor<2,dim> &h_e, + const double &v_s0, + const Tensor<1,dim> &g_s0, + const Tensor<2,dim> &h_s0, + const double &v_s1, + const Tensor<1,dim> &g_s1, + const Tensor<2,dim> &h_s1) +{ + const double v_f = func.value(p); + const Tensor<1,dim> g_f = func.gradient(p); + const SymmetricTensor<2,dim> h_f = func.hessian(p); + + // product rule: + const double v_s = v_s0 + v_s1 * v_f; + const Tensor<1,dim> g_s = g_s0 + g_s1 * v_f + v_s1 * g_f; + Tensor<2,dim> op = outer_product(g_s1 , g_f ); + + const SymmetricTensor<2,dim> sp = symmetrize (op) * 2.0; // symmetrize does [s+s^T]/2 + // Hessians should be symmetric, but due to round off errors and the fact + // that deal.II stores hessians in tensors, we may end up with failing check + // for symmetry in SymmetricTensor class. + // For a moment stick with usual Tensor + const Tensor<2,dim> h_s = h_s0 + h_s1 * v_f + sp + v_s1 * h_f; + + const double v_d = v_s - v_e; + AssertThrow ( std::abs(v_d) < eps, + ExcInternalError()); + const Tensor<1,dim> g_d = g_s - g_e; + AssertThrow ( g_d.norm() < eps, + ExcInternalError()); + + // see note above. + const Tensor<2,dim> h_d = h_s - h_e; + AssertThrow ( h_d.norm() < eps, + ExcInternalError()); +} + +template +void test(const FiniteElement &fe1, + const FiniteElement &fe2, + const Quadrature &volume_quad, + const Quadrature &face_quad, + const bool distort) +{ + Triangulation triangulation; + DoFHandler dof_handler_enriched(triangulation); + DoFHandler dof_handler_system(triangulation); + + EnrichmentFunction function; + FE_Enriched fe_enriched(fe1, + fe2, + &function); + FESystem fe_system(fe1,1, + fe2,1); + + { + Point p1,p2; + std::vector repetitions(dim); + for (unsigned int d=0; d< dim; d++) + { + p1[d]=-1.0; + p2[d]=2.0; + repetitions[d] = 3; + } + GridGenerator::subdivided_hyper_rectangle(triangulation, + repetitions, + p1,p2); + + if (distort) + GridTools::distort_random (0.1,triangulation); + } + + dof_handler_enriched.distribute_dofs(fe_enriched); + dof_handler_system.distribute_dofs(fe_system); + + FEValues fe_values_enriched(fe_enriched, volume_quad, + update_values | update_gradients | update_hessians ); + + FEValues fe_values_system(fe_system, volume_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + FEFaceValues fe_face_values_enriched(fe_enriched, face_quad, + update_values | update_gradients | update_hessians); + + FEFaceValues fe_face_values_system (fe_system, face_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + const unsigned int dofs_per_cell = fe_enriched.dofs_per_cell; + Assert (fe_enriched.dofs_per_cell == fe_system.dofs_per_cell, + ExcInternalError()); + + typename DoFHandler::active_cell_iterator + cell_enriched = dof_handler_enriched.begin_active (), + endc_enriched = dof_handler_enriched.end (), + cell_system = dof_handler_system.begin_active(), + endc_system = dof_handler_system.end(); + for (; cell_enriched!=endc_enriched; ++cell_enriched, ++cell_system) + { + fe_values_enriched.reinit (cell_enriched); + fe_values_system.reinit (cell_system); + const unsigned int n_q_points = volume_quad.size(); + const std::vector > &q_points = fe_values_system.get_quadrature_points(); + + // check shape functions + for (unsigned int i=0; i::faces_per_cell; ++face) + { + fe_face_values_enriched.reinit (cell_enriched,face); + fe_face_values_system.reinit (cell_system,face); + const unsigned int n_q_points_face = face_quad.size(); + const std::vector > &q_points = fe_face_values_system.get_quadrature_points(); + + for (unsigned int i=0; i(3),FE_Q(2),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),true); + } + + { + const unsigned int dim = 3; + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),QGauss(2),QGauss(2),true); + } + + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output new file mode 100644 index 0000000000..f2b98e8f5b --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system.with_cxx11=on.output @@ -0,0 +1,5 @@ + +DEAL::Ok +DEAL::Ok +DEAL::Ok +DEAL::Ok diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.cc b/tests/fe/fe_enriched_compare_to_fe_system_2.cc new file mode 100644 index 0000000000..ca0e3c8468 --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system_2.cc @@ -0,0 +1,403 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 FE_Enriched class by comparing with FE_System for multiple enrichment +// functions. + +#include "../tests.h" + +#include +#include + +#include + +#include +#include + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include +#include + +const double eps = 1e-10; + +// create a mesh without hanging nodes, +// move nodes randomly +// and compare FEEnriched to FESystem with explicit treatment of +// the product rule for +// valus, gradients, hessians on +// elements and faces. +// The comparison is straight forward because local dofs are enumerated +// in the same way for FE_System and FEEnriched. +// Same as fe_enriched_compare_to_fe_system.cc but now multiple enrichment functions. +const unsigned int patches = 10; + +using namespace dealii; + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin) + : Function(1), + origin(origin) + {} + + virtual double value(const Point &point, + const unsigned int component = 0 ) const + { + Tensor<1,dim> d = point-origin; + return std::exp(-d.norm()); + } + + virtual Tensor<1,dim> gradient(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> d = point -origin; + Assert (d.norm() > 0, + dealii::ExcMessage("gradient is not defined at zero")); + d *= -std::exp(-d.norm())/d.norm(); + return d; + } + + virtual SymmetricTensor<2,dim> hessian (const Point &p, + const unsigned int component=0) const + { + Tensor<1,dim> dir = p - origin; + const double r = dir.norm(); + Assert (r>0.0, ExcMessage("r is not positive")); + dir/=r; + SymmetricTensor<2,dim> dir_x_dir; + for (unsigned int i=0; i()/r ); + } + +private: + Point origin; + +}; + +/** + * + * @param p point + * @param func function + * @param v_e value of enriched FE + * @param g_e gradient of enriched FE + * @param h_e hessian of enriched FE + * @param v_s0 value of 0th component of FE_System + * @param g_s0 gradient of 0th component of FE_System + * @param h_s0 hessian of 0th component of FE_System + * @param v_s1 value of 1st component of FE_System + * @param g_s1 gradient of 1st component of FE_System + * @param h_s1 hessian of 1st component of FE_System + */ +template +void check_consistency(const Point &p, + const Function &func1, + const Function &func2, + const Function &func3, + const double &v_e, + const Tensor<1,dim> &g_e, + const Tensor<2,dim> &h_e, + const double &v_s0, + const Tensor<1,dim> &g_s0, + const Tensor<2,dim> &h_s0, + const double &v_s1, + const Tensor<1,dim> &g_s1, + const Tensor<2,dim> &h_s1, + const double &v_s2, + const Tensor<1,dim> &g_s2, + const Tensor<2,dim> &h_s2, + const double &v_s3, + const Tensor<1,dim> &g_s3, + const Tensor<2,dim> &h_s3) +{ + const double v_f1 = func1.value(p); + const Tensor<1,dim> g_f1 = func1.gradient(p); + const SymmetricTensor<2,dim> h_f1 = func1.hessian(p); + + const double v_f2 = func2.value(p); + const Tensor<1,dim> g_f2 = func2.gradient(p); + const SymmetricTensor<2,dim> h_f2 = func2.hessian(p); + + const double v_f3 = func3.value(p); + const Tensor<1,dim> g_f3 = func3.gradient(p); + const SymmetricTensor<2,dim> h_f3 = func3.hessian(p); + + + + // product rule: + const double v_s = v_s0 + v_s1 * v_f1 + + v_s2 * v_f2 + + v_s3 * v_f3; + const Tensor<1,dim> g_s = g_s0 + g_s1 * v_f1 + v_s1 * g_f1 + + g_s2 * v_f2 + v_s2 * g_f2 + + g_s3 * v_f3 + v_s3 * g_f3; + Tensor<2,dim> op1 = outer_product(g_s1 , g_f1 ), + op2 = outer_product(g_s2 , g_f2 ), + op3 = outer_product(g_s3 , g_f3 ); + + const SymmetricTensor<2,dim> sp1 = symmetrize (op1) * 2.0; // symmetrize does [s+s^T]/2 + const SymmetricTensor<2,dim> sp2 = symmetrize (op2) * 2.0; // symmetrize does [s+s^T]/2 + const SymmetricTensor<2,dim> sp3 = symmetrize (op3) * 2.0; // symmetrize does [s+s^T]/2 + // Hessians should be symmetric, but due to round off errors and the fact + // that deal.II stores hessians in tensors, we may end up with failing check + // for symmetry in SymmetricTensor class. + // For a moment stick with usual Tensor + const Tensor<2,dim> h_s = h_s0 + h_s1 * v_f1 + sp1 + v_s1 * h_f1 + + h_s2 * v_f2 + sp2 + v_s2 * h_f2 + + h_s3 * v_f3 + sp3 + v_s3 * h_f3; + + const double v_d = v_s - v_e; + AssertThrow ( std::abs(v_d) < eps, + ExcInternalError()); + const Tensor<1,dim> g_d = g_s - g_e; + AssertThrow ( g_d.norm() < eps, + ExcInternalError()); + + // see note above. + const Tensor<2,dim> h_d = h_s - h_e; + AssertThrow ( h_d.norm() < eps, + ExcInternalError()); +} + +/** + * {fe_en1,1}, {fe_en2,2} => 3 functions + */ +template +void test(const FiniteElement &fe_base, + const FiniteElement &fe_en1, + const FiniteElement &fe_en2, + const Quadrature &volume_quad, + const Quadrature &face_quad, + const bool distort) +{ + Triangulation triangulation; + DoFHandler dof_handler_enriched(triangulation); + DoFHandler dof_handler_system(triangulation); + + Point p1, p2, p3; + for (unsigned int d=0; d< dim; d++) + { + p1[d] = 0.0; + p2[d] = 0.0; + p3[d] = 0.0; + } + p2[0] = 10.0; + p3[0] = 5.0; + EnrichmentFunction fun1(p1), fun2(p2),fun3(p3); + std::vector *> fe_enrichements(2); + fe_enrichements[0] = &fe_en1; + fe_enrichements[1] = &fe_en2; + std::vector *(const typename Triangulation::cell_iterator &) > > > functions(2); + functions[0].resize(1); + functions[0][0] = [&](const typename Triangulation::cell_iterator &) + { + return &fun1; + }; + functions[1].resize(2); + functions[1][0] = [&](const typename Triangulation::cell_iterator &) + { + return &fun2; + }; + functions[1][1] = [&](const typename Triangulation::cell_iterator &) + { + return &fun3; + }; + + FE_Enriched fe_enriched(&fe_base, + fe_enrichements, + functions); + FESystem fe_system(fe_base,1, + fe_en1, 1, + fe_en2, 2); + + { + Point p1,p2; + std::vector repetitions(dim); + for (unsigned int d=0; d< dim; d++) + { + p1[d]=-1.0; + p2[d]=2.0; + repetitions[d] = 3; + } + GridGenerator::subdivided_hyper_rectangle(triangulation, + repetitions, + p1,p2); + + if (distort) + GridTools::distort_random (0.1,triangulation); + } + + dof_handler_enriched.distribute_dofs(fe_enriched); + dof_handler_system.distribute_dofs(fe_system); + + FEValues fe_values_enriched(fe_enriched, volume_quad, + update_values | update_gradients | update_hessians ); + + FEValues fe_values_system(fe_system, volume_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + FEFaceValues fe_face_values_enriched(fe_enriched, face_quad, + update_values | update_gradients | update_hessians); + + FEFaceValues fe_face_values_system (fe_system, face_quad, + update_values | update_gradients | update_hessians | update_quadrature_points); + + AssertThrow(dof_handler_enriched.n_dofs() == dof_handler_system.n_dofs(), + ExcInternalError()); + Vector solution(dof_handler_enriched.n_dofs()); + + const unsigned int dofs_per_cell = fe_enriched.dofs_per_cell; + Assert (fe_enriched.dofs_per_cell == fe_system.dofs_per_cell, + ExcInternalError()); + + typename DoFHandler::active_cell_iterator + cell_enriched = dof_handler_enriched.begin_active (), + endc_enriched = dof_handler_enriched.end (), + cell_system = dof_handler_system.begin_active(), + endc_system = dof_handler_system.end(); + for (; cell_enriched!=endc_enriched; ++cell_enriched, ++cell_system) + { + fe_values_enriched.reinit (cell_enriched); + fe_values_system.reinit (cell_system); + const unsigned int n_q_points = volume_quad.size(); + const std::vector > &q_points = fe_values_system.get_quadrature_points(); + + // check shape functions + for (unsigned int i=0; i::faces_per_cell; ++face) + { + fe_face_values_enriched.reinit (cell_enriched,face); + fe_face_values_system.reinit (cell_system,face); + const unsigned int n_q_points_face = face_quad.size(); + const std::vector > &q_points = fe_face_values_system.get_quadrature_points(); + + for (unsigned int i=0; i(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),true); + } + + { + const unsigned int dim = 3; + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),false); + test(FE_Q(3),FE_Q(2),FE_Q(1),QGauss(2),QGauss(2),true); + } + + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output b/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output new file mode 100644 index 0000000000..f2b98e8f5b --- /dev/null +++ b/tests/fe/fe_enriched_compare_to_fe_system_2.with_cxx11=on.output @@ -0,0 +1,5 @@ + +DEAL::Ok +DEAL::Ok +DEAL::Ok +DEAL::Ok diff --git a/tests/fe/fe_enriched_step-36.cc b/tests/fe/fe_enriched_step-36.cc new file mode 100644 index 0000000000..a22183a967 --- /dev/null +++ b/tests/fe/fe_enriched_step-36.cc @@ -0,0 +1,779 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 FE_Enriched in real-life application on eigenvalue problem similar +// to Step-36. That involves assembly (shape values and gradients) and +// error estimator (Kelly - > face gradients) and MPI run. + +#include "../tests.h" + +#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 + +const unsigned int dim = 3; + +using namespace dealii; + +/** + * Coulomb potential + */ +template +class PotentialFunction : public dealii::Function +{ +public: + PotentialFunction() + : dealii::Function(1) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +double PotentialFunction::value(const dealii::Point &p, + const unsigned int ) const +{ + return -1.0 / std::sqrt(p.square()); +} + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin, + const double &Z, + const double &radius) + : Function(1), + origin(origin), + Z(Z), + radius(radius) + {} + + virtual double value(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> dist = point-origin; + const double r = dist.norm(); + return std::exp(-Z*r); + } + + bool is_enriched(const Point &point) const + { + if (origin.distance(point) < radius) + return true; + else + return false; + } + + 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(); + dist/=r; + return -Z*std::exp(-Z*r)*dist; + } + +private: + /** + * origin + */ + const Point origin; + + /** + * charge + */ + const double Z; + + /** + * enrichment radius + */ + const double radius; +}; + +namespace Step36 +{ + + /** + * Main class + */ + template + class EigenvalueProblem + { + public: + EigenvalueProblem (); + virtual ~EigenvalueProblem(); + void run (); + + private: + bool cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const; + + std::pair setup_system (); + void assemble_system (); + std::pair solve (); + void constrain_pou_dofs(); + void estimate_error (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection q_collection; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector eigenfunctions_locally_relevant; + std::vector eigenvalues; + PETScWrappers::MPI::SparseMatrix stiffness_matrix, mass_matrix; + + ConstraintMatrix constraints; + + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + + ConditionalOStream pcout; + + const unsigned int number_of_eigenvalues; + + PotentialFunction potential; + + EnrichmentFunction enrichment; + + const FEValuesExtractors::Scalar fe_extractor; + const unsigned int fe_group; + const unsigned int fe_fe_index; + const unsigned int fe_material_id; + const FEValuesExtractors::Scalar pou_extractor; + const unsigned int pou_group; + const unsigned int pou_fe_index; + const unsigned int pou_material_id; + + std::vector > vec_estimated_error_per_cell; + Vector estimated_error_per_cell; + }; + + template + EigenvalueProblem::EigenvalueProblem () + : + dof_handler (triangulation), + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process(dealii::Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout (std::cout, + (this_mpi_process == 0)), + number_of_eigenvalues(1), + enrichment(Point(), + /*Z*/1.0, + /*radius*/2.5), // radius is set such that 8 cells are marked as enriched + fe_extractor(/*dofs start at...*/0), + fe_group(/*in FE*/0), + fe_fe_index(0), + fe_material_id(0), + pou_extractor(/*dofs start at (scalar fields!)*/1), + pou_group(/*FE*/1), + pou_fe_index(1), + pou_material_id(1) + { + dealii::GridGenerator::hyper_cube (triangulation, -10, 10); + triangulation.refine_global (2); // 64 cells + + for (auto cell= triangulation.begin_active(); cell != triangulation.end(); ++cell) + if (std::sqrt(cell->center().square()) < 5.0) + cell->set_refine_flag(); + + triangulation.execute_coarsening_and_refinement(); // 120 cells + + q_collection.push_back(QGauss(4)); + q_collection.push_back(QGauss(10)); + + // usual elements (active_fe_index ==0): + fe_collection.push_back (FE_Enriched (FE_Q(2))); + + + // enriched elements (active_fe_index==1): + fe_collection.push_back (FE_Enriched (FE_Q(2), + FE_Q(1), + &enrichment)); + + // assign fe index in the constructor so that FE/FE+POU is determined + // on the coarsest mesh to avoid issues like + // +---------+----+----+ + // | | fe | | + // | fe +----+----+ + // | | pou| | + // +---------+----+----+ + // see discussion in Step46. + for (typename hp::DoFHandler::cell_iterator cell= dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (enrichment.is_enriched(cell->center())) + cell->set_material_id(pou_material_id); + else + cell->set_material_id(fe_material_id); + } + + template + std::pair + EigenvalueProblem::setup_system () + { + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + { + if (cell->material_id() == fe_material_id) + cell->set_active_fe_index (fe_fe_index); + else if (cell->material_id() == pou_material_id) + cell->set_active_fe_index (pou_fe_index); + else + Assert (false, ExcNotImplemented()); + } + + GridTools::partition_triangulation (n_mpi_processes, triangulation); + dof_handler.distribute_dofs (fe_collection); + DoFRenumbering::subdomain_wise (dof_handler); + std::vector locally_owned_dofs_per_processor + = DoFTools::locally_owned_dofs_per_subdomain (dof_handler); + locally_owned_dofs = locally_owned_dofs_per_processor[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); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (1), + constraints); + constraints.close (); + + // Initialise the stiffness and mass matrices + DynamicSparsityPattern csp (locally_relevant_dofs); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, + /* keep constrained dofs */ 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_processor[i].n_elements(); + + SparsityTools::distribute_sparsity_pattern + (csp, + n_locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + + stiffness_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + mass_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + // reinit vectors + eigenfunctions.resize (number_of_eigenvalues); + eigenfunctions_locally_relevant.resize(number_of_eigenvalues); + vec_estimated_error_per_cell.resize(number_of_eigenvalues); + for (unsigned int i=0; i::cell_iterator cell= dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + n_pou_cells++; + else + n_fem_cells++; + + return std::make_pair (n_fem_cells, + n_pou_cells); + } + + template + bool + EigenvalueProblem::cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const + { + return cell->material_id() == pou_material_id; + } + + template + void + EigenvalueProblem::constrain_pou_dofs() + { + std::vector local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face); + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (!cell->at_boundary(f)) + { + bool face_is_on_interface = false; + if ((cell->neighbor(f)->has_children() == false /* => it is active */) + && + (! cell_is_pou(cell->neighbor(f)))) + face_is_on_interface = true; + else if (cell->neighbor(f)->has_children() == true) + { + // The neighbor does + // have + // children. See if + // any of the cells + // on the other + // side are vanilla FEM + for (unsigned int sf=0; sfface(f)->n_children(); ++sf) + if (!cell_is_pou (cell->neighbor_child_on_subface(f, sf))) + { + face_is_on_interface = true; + break; + } + } + // add constraints + if (face_is_on_interface) + { + cell->face(f)->get_dof_indices (local_face_dof_indices, pou_fe_index); + for (unsigned int i=0; i 0) + constraints.add_line (local_face_dof_indices[i]); + } + } + } + + template + void + EigenvalueProblem::assemble_system() + { + stiffness_matrix = 0; + mass_matrix = 0; + + dealii::FullMatrix cell_stiffness_matrix; + dealii::FullMatrix cell_mass_matrix; + std::vector local_dof_indices; + std::vector potential_values; + std::vector enrichment_values; + std::vector> enrichment_gradients; + + hp::FEValues fe_values_hp(fe_collection, q_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + 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; + + potential_values.resize(n_q_points); + enrichment_values.resize(n_q_points); + enrichment_gradients.resize(n_q_points); + + local_dof_indices.resize (dofs_per_cell); + cell_stiffness_matrix.reinit (dofs_per_cell,dofs_per_cell); + cell_mass_matrix.reinit (dofs_per_cell,dofs_per_cell); + + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + potential.value_list (fe_values.get_quadrature_points(), + potential_values); + + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + + stiffness_matrix.compress (dealii::VectorOperation::add); + mass_matrix.compress (dealii::VectorOperation::add); + } + + template + std::pair + EigenvalueProblem::solve() + { + dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false); + dealii::SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control, + mpi_communicator); + + eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); + eigensolver.set_problem_type (EPS_GHEP); + + eigensolver.solve (stiffness_matrix, mass_matrix, + eigenvalues, eigenfunctions, + eigenfunctions.size()); + + for (unsigned int i = 0; i < eigenfunctions.size(); i++) + { + constraints.distribute(eigenfunctions[i]); + eigenfunctions_locally_relevant[i] = eigenfunctions[i]; + } + + return std::make_pair (solver_control.last_step (), + solver_control.last_value ()); + + } + + template + EigenvalueProblem::~EigenvalueProblem() + { + dof_handler.clear (); + } + + template + void EigenvalueProblem::estimate_error () + { + { + std::vector sol (number_of_eigenvalues); + std::vector *> error (number_of_eigenvalues); + + for (unsigned int i = 0; i < number_of_eigenvalues; i++) + { + sol[i] = &eigenfunctions_locally_relevant[i]; + error[i] = &vec_estimated_error_per_cell[i]; + } + + hp::QCollection face_quadrature_formula; + face_quadrature_formula.push_back (dealii::QGauss(3)); + face_quadrature_formula.push_back (dealii::QGauss(3)); + + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_formula, + typename FunctionMap::type(), + sol, + error); + } + + // sum up for a global: + for (unsigned int c=0; c < estimated_error_per_cell.size(); c++) + { + double er = 0.0; + for (unsigned int i = 0; i < number_of_eigenvalues; i++) + er += vec_estimated_error_per_cell[i][c] * vec_estimated_error_per_cell[i][c]; + + estimated_error_per_cell[c] = sqrt(er); + } + } + + template + void EigenvalueProblem::refine_grid () + { + const double threshold = 0.9 * estimated_error_per_cell.linfty_norm(); + GridRefinement::refine (triangulation, + estimated_error_per_cell, + threshold); + + triangulation.prepare_coarsening_and_refinement (); + triangulation.execute_coarsening_and_refinement (); + } + + template + class Postprocessor : public DataPostprocessorScalar + { + public: + Postprocessor(const EnrichmentFunction &enrichment); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + private: + const EnrichmentFunction &enrichment; + }; + + template + Postprocessor::Postprocessor(const EnrichmentFunction &enrichment) + : + DataPostprocessorScalar ("total_solution", + update_values | update_q_points), + enrichment(enrichment) + {} + + template + void + Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const + { + const unsigned int n_quadrature_points = uh.size(); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + for (unsigned int q=0; q + void EigenvalueProblem::output_results (const unsigned int cycle) const + { + dealii::Vector fe_index(triangulation.n_active_cells()); + { + typename dealii::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(); + } + + Assert (cycle < 10, ExcNotImplemented()); + if (this_mpi_process==0) + { + std::string filename = "solution-"; + filename += ('0' + cycle); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + Postprocessor postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (eigenfunctions_locally_relevant[0], "solution"); + data_out.build_patches (6); + data_out.write_vtk (output); + } + + // second output without making mesh finer + if (this_mpi_process==0) + { + std::string filename = "mesh-"; + filename += ('0' + cycle); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (fe_index, "fe_index"); + data_out.add_data_vector (estimated_error_per_cell, "estimated_error"); + data_out.build_patches (); + data_out.write_vtk (output); + } + + // scalar data for plotting + //output scalar data (eigenvalues, energies, ndofs, etc) + if (this_mpi_process == 0) + { + const std::string scalar_fname = "scalar-data.txt"; + + std::ofstream output (scalar_fname.c_str(), + std::ios::out | + (cycle==0 + ? + std::ios::trunc : std::ios::app)); + + std::streamsize max_precision = std::numeric_limits::digits10; + + std::string sep(" "); + output << cycle << sep + << std::setprecision(max_precision) + << triangulation.n_active_cells() << sep + << dof_handler.n_dofs () << sep + << std::scientific; + + for (unsigned int i = 0; i < eigenvalues.size(); i++) + output << eigenvalues[i] << sep; + + output< + void + EigenvalueProblem::run() + { + for (unsigned int cycle = 0; cycle < 5; cycle++) + { + pcout << "Cycle "< n_cells = setup_system (); + + pcout << " Number of active cells: " + << triangulation.n_active_cells () + << std::endl + << " FE / POUFE : " + << n_cells.first << " / " << n_cells.second + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs () + << std::endl; + + assemble_system (); + + const std::pair res = solve (); + AssertThrow(res.second < 5e-8, + ExcInternalError()); + + estimate_error (); + // output_results(cycle); + refine_grid (); + + pcout << std::endl; + for (unsigned int i=0; i step36; + dealii::PETScWrappers::set_option_value("-eps_target","-1.0"); + dealii::PETScWrappers::set_option_value("-st_type","sinvert"); + dealii::PETScWrappers::set_option_value("-st_ksp_type","cg"); + dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg"); + dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); + step36.run(); + } + + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output new file mode 100644 index 0000000000..fa33870cc1 --- /dev/null +++ b/tests/fe/fe_enriched_step-36.with_petsc=on.with_slepc=on.with_cxx11=on.output @@ -0,0 +1,30 @@ +Cycle 0 + Number of active cells: 120 + FE / POUFE : 120 / 8 + Number of degrees of freedom: 1432 + + Eigenvalue 0 : -0.483398 +Cycle 1 + Number of active cells: 176 + FE / POUFE : 120 / 72 + Number of degrees of freedom: 2206 + + Eigenvalue 0 : -0.498786 +Cycle 2 + Number of active cells: 512 + FE / POUFE : 504 / 72 + Number of degrees of freedom: 5734 + + Eigenvalue 0 : -0.499788 +Cycle 3 + Number of active cells: 568 + FE / POUFE : 504 / 136 + Number of degrees of freedom: 6508 + + Eigenvalue 0 : -0.499818 +Cycle 4 + Number of active cells: 736 + FE / POUFE : 504 / 328 + Number of degrees of freedom: 8542 + + Eigenvalue 0 : -0.499947 diff --git a/tests/fe/fe_enriched_step-36b.cc b/tests/fe/fe_enriched_step-36b.cc new file mode 100644 index 0000000000..4a92b81d87 --- /dev/null +++ b/tests/fe/fe_enriched_step-36b.cc @@ -0,0 +1,903 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2016 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 FE_Enriched in real-life application on eigenvalue problem similar +// to Step-36. That involves assembly (shape values and gradients) and +// error estimator (Kelly - > face gradients) and MPI run. +// +// same as fe_enriched_step-36, but using FESystem only. Note that the results +// of KellyErrorEstimator won't be the same, thus this tests does only a single +// refinement cycle where the refinement path is the same. + +#include "../tests.h" + +#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 + +const unsigned int dim = 3; + +using namespace dealii; + +/** + * Coulomb potential + */ +template +class PotentialFunction : public dealii::Function +{ +public: + PotentialFunction() + : dealii::Function(1) + {} + + virtual double value(const dealii::Point &point, + const unsigned int component = 0 ) const; +}; + +template +double PotentialFunction::value(const dealii::Point &p, + const unsigned int ) const +{ + return -1.0 / std::sqrt(p.square()); +} + +template +class EnrichmentFunction : public Function +{ +public: + EnrichmentFunction(const Point &origin, + const double &Z, + const double &radius) + : Function(1), + origin(origin), + Z(Z), + radius(radius) + {} + + virtual double value(const Point &point, + const unsigned int component = 0) const + { + Tensor<1,dim> dist = point-origin; + const double r = dist.norm(); + return std::exp(-Z*r); + } + + bool is_enriched(const Point &point) const + { + if (origin.distance(point) < radius) + return true; + else + return false; + } + + 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(); + dist/=r; + return -Z*std::exp(-Z*r)*dist; + } + +private: + /** + * origin + */ + const Point origin; + + /** + * charge + */ + const double Z; + + /** + * enrichment radius + */ + const double radius; +}; + +namespace Step36 +{ + + /** + * Main class + */ + template + class EigenvalueProblem + { + public: + EigenvalueProblem (); + virtual ~EigenvalueProblem(); + void run (); + + private: + bool cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const; + + std::pair setup_system (); + void assemble_system (); + std::pair solve (); + void constrain_pou_dofs(); + void estimate_error (); + void refine_grid (); + void output_results (const unsigned int cycle) const; + + Triangulation triangulation; + hp::DoFHandler dof_handler; + hp::FECollection fe_collection; + hp::QCollection q_collection; + + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + + std::vector eigenfunctions; + std::vector eigenfunctions_locally_relevant; + std::vector eigenvalues; + PETScWrappers::MPI::SparseMatrix stiffness_matrix, mass_matrix; + + ConstraintMatrix constraints; + + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + + ConditionalOStream pcout; + + const unsigned int number_of_eigenvalues; + + PotentialFunction potential; + + EnrichmentFunction enrichment; + + const FEValuesExtractors::Scalar fe_extractor; + const unsigned int fe_group; + const unsigned int fe_fe_index; + const unsigned int fe_material_id; + const FEValuesExtractors::Scalar pou_extractor; + const unsigned int pou_group; + const unsigned int pou_fe_index; + const unsigned int pou_material_id; + + std::vector > vec_estimated_error_per_cell; + Vector estimated_error_per_cell; + }; + + template + EigenvalueProblem::EigenvalueProblem () + : + dof_handler (triangulation), + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes(dealii::Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process(dealii::Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout (std::cout, + (this_mpi_process == 0)), + number_of_eigenvalues(1), + enrichment(Point(), + /*Z*/1.0, + /*radius*/2.5), // radius is set such that 8 cells are marked as enriched + fe_extractor(/*dofs start at...*/0), + fe_group(/*in FE*/0), + fe_fe_index(0), + fe_material_id(0), + pou_extractor(/*dofs start at (scalar fields!)*/1), + pou_group(/*FE*/1), + pou_fe_index(1), + pou_material_id(1) + { + dealii::GridGenerator::hyper_cube (triangulation, -10, 10); + triangulation.refine_global (2); // 64 cells + + for (auto cell= triangulation.begin_active(); cell != triangulation.end(); ++cell) + if (std::sqrt(cell->center().square()) < 5.0) + cell->set_refine_flag(); + + triangulation.execute_coarsening_and_refinement(); // 120 cells + + q_collection.push_back(QGauss(4)); + q_collection.push_back(QGauss(10)); + + // usual elements (active_fe_index ==0): + fe_collection.push_back (FESystem (FE_Q(2), 1, + FE_Nothing(), 1)); + + // enriched elements (active_fe_index==1): + fe_collection.push_back (FESystem (FE_Q(2), 1, + FE_Q(1), 1)); + + // assign fe index in the constructor so that FE/FE+POU is determined + // on the coarsest mesh to avoid issues like + // +---------+----+----+ + // | | fe | | + // | fe +----+----+ + // | | pou| | + // +---------+----+----+ + // see discussion in Step46. + for (typename hp::DoFHandler::cell_iterator cell= dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (enrichment.is_enriched(cell->center())) + cell->set_material_id(pou_material_id); + else + cell->set_material_id(fe_material_id); + } + + template + std::pair + EigenvalueProblem::setup_system () + { + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + { + if (cell->material_id() == fe_material_id) + cell->set_active_fe_index (fe_fe_index); + else if (cell->material_id() == pou_material_id) + cell->set_active_fe_index (pou_fe_index); + else + Assert (false, ExcNotImplemented()); + } + + GridTools::partition_triangulation (n_mpi_processes, triangulation); + dof_handler.distribute_dofs (fe_collection); + DoFRenumbering::subdomain_wise (dof_handler); + std::vector locally_owned_dofs_per_processor + = DoFTools::locally_owned_dofs_per_subdomain (dof_handler); + locally_owned_dofs = locally_owned_dofs_per_processor[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); + VectorTools::interpolate_boundary_values (dof_handler, + 0, + ZeroFunction (2), + constraints); + constrain_pou_dofs(); + constraints.close (); + + // Initialise the stiffness and mass matrices + DynamicSparsityPattern csp (locally_relevant_dofs); + DoFTools::make_sparsity_pattern (dof_handler, csp, + constraints, + /* keep constrained dofs */ 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_processor[i].n_elements(); + + SparsityTools::distribute_sparsity_pattern + (csp, + n_locally_owned_dofs, + mpi_communicator, + locally_relevant_dofs); + + stiffness_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + mass_matrix.reinit (locally_owned_dofs, + locally_owned_dofs, + csp, + mpi_communicator); + + // reinit vectors + eigenfunctions.resize (number_of_eigenvalues); + eigenfunctions_locally_relevant.resize(number_of_eigenvalues); + vec_estimated_error_per_cell.resize(number_of_eigenvalues); + for (unsigned int i=0; i::cell_iterator cell= dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + n_pou_cells++; + else + n_fem_cells++; + + return std::make_pair (n_fem_cells, + n_pou_cells); + } + + template + bool + EigenvalueProblem::cell_is_pou(const typename hp::DoFHandler::cell_iterator &cell) const + { + return cell->material_id() == pou_material_id; + } + + template + void + EigenvalueProblem::constrain_pou_dofs() + { + std::vector local_face_dof_indices(fe_collection[pou_fe_index].dofs_per_face); + for (typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(); + cell != dof_handler.end(); ++cell) + if (cell_is_pou(cell)) + for (unsigned int f=0; f::faces_per_cell; ++f) + if (!cell->at_boundary(f)) + { + bool face_is_on_interface = false; + if ((cell->neighbor(f)->has_children() == false /* => it is active */) + && + (! cell_is_pou(cell->neighbor(f)))) + face_is_on_interface = true; + else if (cell->neighbor(f)->has_children() == true) + { + // The neighbor does + // have + // children. See if + // any of the cells + // on the other + // side are vanilla FEM + for (unsigned int sf=0; sfface(f)->n_children(); ++sf) + if (!cell_is_pou (cell->neighbor_child_on_subface(f, sf))) + { + face_is_on_interface = true; + break; + } + } + // add constraints + if (face_is_on_interface) + { + cell->face(f)->get_dof_indices (local_face_dof_indices, pou_fe_index); + for (unsigned int i=0; i 0) + constraints.add_line (local_face_dof_indices[i]); + } + } + } + + template + void + EigenvalueProblem::assemble_system() + { + stiffness_matrix = 0; + mass_matrix = 0; + + dealii::FullMatrix cell_stiffness_matrix; + dealii::FullMatrix cell_mass_matrix; + std::vector local_dof_indices; + std::vector potential_values; + std::vector enrichment_values; + std::vector> enrichment_gradients; + + hp::FEValues fe_values_hp(fe_collection, q_collection, + update_values | update_gradients | + update_quadrature_points | update_JxW_values); + + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (; cell!=endc; ++cell) + 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; + + potential_values.resize(n_q_points); + enrichment_values.resize(n_q_points); + enrichment_gradients.resize(n_q_points); + + local_dof_indices.resize (dofs_per_cell); + cell_stiffness_matrix.reinit (dofs_per_cell,dofs_per_cell); + cell_mass_matrix.reinit (dofs_per_cell,dofs_per_cell); + + cell_stiffness_matrix = 0; + cell_mass_matrix = 0; + + potential.value_list (fe_values.get_quadrature_points(), + potential_values); + + if (cell->active_fe_index() == 0) // plain FE + { + for (unsigned int q_point=0; q_pointactive_fe_index() == 1, ExcInternalError()); + + enrichment.value_list (fe_values.get_quadrature_points(), + enrichment_values); + enrichment.gradient_list (fe_values.get_quadrature_points(), + enrichment_gradients); + + for (unsigned int i=0; iget_fe().system_to_base_index(i).first.first; + for (unsigned int j=i; jget_fe().system_to_base_index(j).first.first; + + if ( i_group == fe_group && j_group == fe_group ) // fe - fe + { + for (unsigned int q_point=0; q_pointget_dof_indices (local_dof_indices); + + constraints + .distribute_local_to_global (cell_stiffness_matrix, + local_dof_indices, + stiffness_matrix); + constraints + .distribute_local_to_global (cell_mass_matrix, + local_dof_indices, + mass_matrix); + } + + stiffness_matrix.compress (dealii::VectorOperation::add); + mass_matrix.compress (dealii::VectorOperation::add); + } + + template + std::pair + EigenvalueProblem::solve() + { + dealii::SolverControl solver_control (dof_handler.n_dofs(), 1e-9,false,false); + dealii::SLEPcWrappers::SolverKrylovSchur eigensolver (solver_control, + mpi_communicator); + + eigensolver.set_which_eigenpairs (EPS_SMALLEST_REAL); + eigensolver.set_problem_type (EPS_GHEP); + + eigensolver.solve (stiffness_matrix, mass_matrix, + eigenvalues, eigenfunctions, + eigenfunctions.size()); + + for (unsigned int i = 0; i < eigenfunctions.size(); i++) + { + constraints.distribute(eigenfunctions[i]); + eigenfunctions_locally_relevant[i] = eigenfunctions[i]; + } + + return std::make_pair (solver_control.last_step (), + solver_control.last_value ()); + + } + + template + EigenvalueProblem::~EigenvalueProblem() + { + dof_handler.clear (); + } + + template + void EigenvalueProblem::estimate_error () + { + { + std::vector sol (number_of_eigenvalues); + std::vector *> error (number_of_eigenvalues); + + for (unsigned int i = 0; i < number_of_eigenvalues; i++) + { + sol[i] = &eigenfunctions_locally_relevant[i]; + error[i] = &vec_estimated_error_per_cell[i]; + } + + hp::QCollection face_quadrature_formula; + face_quadrature_formula.push_back (dealii::QGauss(3)); + face_quadrature_formula.push_back (dealii::QGauss(3)); + + KellyErrorEstimator::estimate (dof_handler, + face_quadrature_formula, + typename FunctionMap::type(), + sol, + error); + } + + // sum up for a global: + for (unsigned int c=0; c < estimated_error_per_cell.size(); c++) + { + double er = 0.0; + for (unsigned int i = 0; i < number_of_eigenvalues; i++) + er += vec_estimated_error_per_cell[i][c] * vec_estimated_error_per_cell[i][c]; + + estimated_error_per_cell[c] = sqrt(er); + } + } + + template + void EigenvalueProblem::refine_grid () + { + const double threshold = 0.9 * estimated_error_per_cell.linfty_norm(); + GridRefinement::refine (triangulation, + estimated_error_per_cell, + threshold); + + triangulation.prepare_coarsening_and_refinement (); + triangulation.execute_coarsening_and_refinement (); + } + + template + class Postprocessor : public DataPostprocessorScalar + { + public: + Postprocessor(const EnrichmentFunction &enrichment); + + virtual + void + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &duh, + const std::vector > > &dduh, + const std::vector > &normals, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const; + + private: + const EnrichmentFunction &enrichment; + }; + + template + Postprocessor::Postprocessor(const EnrichmentFunction &enrichment) + : + DataPostprocessorScalar ("total_solution", + update_values | update_q_points), + enrichment(enrichment) + {} + + template + void + Postprocessor:: + compute_derived_quantities_vector (const std::vector > &uh, + const std::vector > > &/*duh*/, + const std::vector > > &/*dduh*/, + const std::vector > &/*normals*/, + const std::vector > &evaluation_points, + std::vector > &computed_quantities) const + { + const unsigned int n_quadrature_points = uh.size(); + Assert (computed_quantities.size() == n_quadrature_points, ExcInternalError()); + for (unsigned int q=0; q + void EigenvalueProblem::output_results (const unsigned int cycle) const + { + dealii::Vector fe_index(triangulation.n_active_cells()); + { + typename dealii::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(); + } + + Assert (cycle < 10, ExcNotImplemented()); + if (this_mpi_process==0) + { + std::string filename = "solution-"; + filename += ('0' + cycle); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + Postprocessor postprocessor(enrichment); // has to live until the DataOut object is destroyed; objects are destroyed in reverse order of declaration + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (eigenfunctions_locally_relevant[0], "solution"); + data_out.add_data_vector (eigenfunctions_locally_relevant[0], postprocessor); + data_out.build_patches (6); + data_out.write_vtk (output); + } + + // second output without making mesh finer + if (this_mpi_process==0) + { + std::string filename = "mesh-"; + filename += ('0' + cycle); + filename += ".vtk"; + std::ofstream output (filename.c_str()); + + DataOut > data_out; + data_out.attach_dof_handler (dof_handler); + data_out.add_data_vector (fe_index, "fe_index"); + data_out.add_data_vector (estimated_error_per_cell, "estimated_error"); + data_out.build_patches (); + data_out.write_vtk (output); + } + + // scalar data for plotting + //output scalar data (eigenvalues, energies, ndofs, etc) + if (this_mpi_process == 0) + { + const std::string scalar_fname = "scalar-data.txt"; + + std::ofstream output (scalar_fname.c_str(), + std::ios::out | + (cycle==0 + ? + std::ios::trunc : std::ios::app)); + + std::streamsize max_precision = std::numeric_limits::digits10; + + std::string sep(" "); + output << cycle << sep + << std::setprecision(max_precision) + << triangulation.n_active_cells() << sep + << dof_handler.n_dofs () << sep + << std::scientific; + + for (unsigned int i = 0; i < eigenvalues.size(); i++) + output << eigenvalues[i] << sep; + + output< + void + EigenvalueProblem::run() + { + for (unsigned int cycle = 0; cycle < 2; cycle++) + { + pcout << "Cycle "< n_cells = setup_system (); + + pcout << " Number of active cells: " + << triangulation.n_active_cells () + << std::endl + << " FE / POUFE : " + << n_cells.first << " / " << n_cells.second + << std::endl + << " Number of degrees of freedom: " + << dof_handler.n_dofs () + << std::endl; + + assemble_system (); + + const std::pair res = solve (); + AssertThrow(res.second < 5e-8, + ExcInternalError()); + + estimate_error (); + // output_results(cycle); + refine_grid (); + + pcout << std::endl; + for (unsigned int i=0; i step36; + dealii::PETScWrappers::set_option_value("-eps_target","-1.0"); + dealii::PETScWrappers::set_option_value("-st_type","sinvert"); + dealii::PETScWrappers::set_option_value("-st_ksp_type","cg"); + dealii::PETScWrappers::set_option_value("-st_pc_type", "gamg"); + dealii::PETScWrappers::set_option_value("-st_ksp_tol", "1e-11"); + step36.run(); + } + + } + 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; + }; +} diff --git a/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output b/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output new file mode 100644 index 0000000000..516893b783 --- /dev/null +++ b/tests/fe/fe_enriched_step-36b.with_petsc=on.with_slepc=on.with_cxx11=on.output @@ -0,0 +1,12 @@ +Cycle 0 + Number of active cells: 120 + FE / POUFE : 120 / 8 + Number of degrees of freedom: 1432 + + Eigenvalue 0 : -0.483398 +Cycle 1 + Number of active cells: 176 + FE / POUFE : 120 / 72 + Number of degrees of freedom: 2206 + + Eigenvalue 0 : -0.498786 -- 2.39.5