From eafc440daf71643c862971f302624aaab1b88cbc Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Wed, 12 Aug 2015 04:50:50 -0500 Subject: [PATCH] Implement FE_Q_Bubbles that is based on TensorProductPolynomialsBubbles --- .../images/fe_q_bubbles_conditioning.png | Bin 0 -> 8062 bytes doc/news/changes.h | 5 + include/deal.II/fe/fe_poly.templates.h | 1 + include/deal.II/fe/fe_q_base.h | 11 +- include/deal.II/fe/fe_q_bubbles.h | 205 +++++++ source/fe/CMakeLists.txt | 2 + source/fe/fe_poly.cc | 1 + source/fe/fe_poly.inst.in | 3 +- source/fe/fe_q_base.cc | 88 +-- source/fe/fe_q_base.inst.in | 1 + source/fe/fe_q_bubbles.cc | 527 ++++++++++++++++++ source/fe/fe_q_bubbles.inst.in | 24 + 12 files changed, 829 insertions(+), 39 deletions(-) create mode 100644 doc/doxygen/images/fe_q_bubbles_conditioning.png create mode 100644 include/deal.II/fe/fe_q_bubbles.h create mode 100644 source/fe/fe_q_bubbles.cc create mode 100644 source/fe/fe_q_bubbles.inst.in diff --git a/doc/doxygen/images/fe_q_bubbles_conditioning.png b/doc/doxygen/images/fe_q_bubbles_conditioning.png new file mode 100644 index 0000000000000000000000000000000000000000..a0618721d80da0bc729e61442af89242265496dc GIT binary patch literal 8062 zcmaiZ2T)U8@NW`AuOd~zAXq>V5Kur$0!Wo69Vv;GPi%lRDM6B5Wl(E z1v?N3f(3!VlU!`f6^NkHFms}2Yw2LjU@({~d3kxP9|8!3We`9N5SdIsfWCbL8Qhv; z5M)7E1OdxXMPXw>AOse~z%nq5p~90VnYSV)`gLUx7z`#OBcrCK zhCm>6b#=|o&N4-;sHr89lPZ49tvIesG8p9X2*!uzN-~)-0%A%7U1xq{gi{Qq%ttJD_AVU7(~PRkzRzX%*8d;-L8+~-Pjn*!^$!W zLF5Q>8fb7lZ#>Bs;domGD+0obh+q-SfyJT-SQG=x!jE9&F$ao~=ffZn7)*T_gzmbW z6~;>5hwf%Zb2p(XjDg5w$l`D~B9VwhB6D+dVK5jTkB7tIF?O!*Or5!X%v>LUKmuP_ zAMo>20#YE5%$oTH1Bc+Ob5Iv0TZ`>eYj?K zZeGl&Ljl^N%J=^omX_B&GIVl_y!GK^)i^Kx?l!I;=PKww8(UOyC55Vx;ghZ5i{mZ^N0x7C|;O~2Z^8zhG1c=(DE8vW`?ow<{mam1l-xfDeLbiI=Tt$*6De4)B zCd&Hjt(}!d`-6o;vEhFIYX_jL_b(5YVeBv8&FoM*M;UvcyP*E;%VvgU9&|U%A!nWF zD7+K;9;^6gokt|UUoy9n2*~gTZM`M-`nxG4{;NKzk$gGZJ{;$GR-;nV|7W5<+~4RC z?T0q$MsTMo4oGjk`u6t1sdEK~W4=Y$#KJ!%f>j<^5yp#0kd8YmX9w?kpw&Y!Fl6E1 zU(dniP=lCJLG97*JuRm@VXsIZ5Qha*TvyJPZ{DFEEo|KBP||#>0j1$ zg{!I)&2iw0Mvqv$RKWVO`z7IFlS8vzE8`CbIyxm9E;j`Lfw1?{o->zBz~e3Py9|Gs zBC^K!3!bIO-W>Pt5B;VUIPdQyb-^0^H7BF-v;A@s1Bc_VgDX}6&j!X`at^tsPSw=} z2eIGzO}+EF&2`C7fP>l&rd}dlmgQDmydZGCELT zzTjO)gL3eJF;4>0>MjQ#fZuEmu`WAa7gDOD@5<==`OBpSRWqR0^Jy2D zJm!pZCsi5EWHE?!xJ=cW@wXUXR9iUGJ)OtVV4{{-!DE&uxEwKS9mq z@b>~*Pkw07ZhOi2h?Rc}lJU8S&;+>i_{^$DhPxbF{jINt9EF~DyvJy$Qa^HvHmjT$ zEI8^Xy%wMN7$Fm|`&=|oIT(m@y9a3qUu|4m@6_^Z-}+0~ zY>_!{@SCHU&9u2)PaUctn9XrR=LR6<2-H-gC!c2l)Dw2%K3k-kgZJ0@3&S^;-eMn> z(z2`R8>Q&ta#+CW-C}eS#4hc(%DM6)b@a@~YEl9kjnTJC@R}XVPI1`C%tSf@_gvv~ zqo3sOZ-PXU-W(;4+n#%YHSz`gO|zXcF7usIfM^HLojmMB>c zj@m)T^OsBaMyy*ePk`01aD~IyWZt5k@ zE&or}u*HF$5|Poyu>rx$gs-&kPgQNQa*Z1N(Q?L{ykyVm!?pxm)c8eD;0%Y~+AGI{a*liQRVo^cX*?!9BFoSk^tvnoE?W*a2Zr~KhnIKLw;)_Z3( z_J?Ry@|@1^!gXKX^5Ic4Oz#xK2*mTxNL3?pE`e6-aXPoXa6S6OwCHGMbp}0$hL7Ce zt>j&B>h$s>?LG81zLmkhX@cpUzrsLy&gmw+j|z7r7}tK<0JaNjCb+y~DlOr1^+w^Eq2-ypYeOeX;f#?2OBy z=HYx=!_}JI4cJrP*GS0oun*tR&6Tz3zrC7=u*~U(1Z;iRy!2hAsx zgMA|8pzf7};_lZm37n5GH@LoSH{;Qj;ubATPu%fajNpt8ZlcpmnY{r%w26b4tEZj! zu=Lb1x_?Uzw#%l+b7c~wl8nVJUiqA|cUR0y|AU;;o6=m=+UnvSn|KnAdALPz{`BMY z4SP~$z7zw@YlcjL6O&Z3@d5_4C4j9P+ibBWE_FPB6t?T{iI;j-iWUaG?s01K3OF%C z{^fAt`f@eB_4cMyFQOQkyKQ?j_xl>z5J5b)(-Za7<2-S~Z}pX)qrZm*1`}0*#7A{^ zX(GUp&wn@Cm@f_N7^yyX1V#pkxdDyW%`hz#Fgu?UQ@(2L#ra=#og;sA!=9AoaZzb|))H+BVlOzJdU2itwBO6-KH@^l9QG7@`W^vR|L9Ne$^Fyw5%pOeR>s!C zbd*Xva{0B*!HQDMyiH*Kp(qD@ZDEES^tyndOHfsoZU97wyaEP1=ft)J^-9lY>o{YF^b`uw!V<6~N1Xv^ZdY&>u62VB}$ z-w#f?JfH*Tx`Xx#$#-Xi0uOgigy7tOu~kJv`fDEOEz0gXKMV`;<^Z?J1$(swZRDU! zo|#@4e^z4Agr@v)(hCaXS$nc=h)AIe#MAlmARP_8A=eJb^5@Hv8ROg{oazRJ#2btScK_k1vYnYD(vez&v~^G&mw zj8TOC&dUF=;T9}i&F!TY8?g9P6XNLqsBF(KsXJXdbof|;bZ}e^s(0HhECa#X9*DlT zHU$T`v!dEE1Pb{x&~#G_PZkT9?*dhR)n-yd*FxId=;y9AtYJ|!EsD8ug=u;~&=AkJ zk4yPFZ6hCEn6tGPJ;*$01w2=1nffAqFAucV!+dF2|~1(bT1B>i z0KJGg0x)w+mZ-UHwIv48ogmDqFh%8EWcGJVTdHywB7gifyAXHnrm+5B)`E7xg)vtC zVe;1cXCsyjvt+syY(k)lx4~bsg!5bmrh{0XWnOb%1+FUNEcc~r32oto{;I0*_)qvX z8Yj!jSVS4@>q75-&xt@})3L0KmsK;+e`Hwj?XIzRV4BvnQ8~=+ul`q!a_FkeT(+P_hbD^((zV!fAjbqK(1G65Vg1^w zUCESv(A&q;(6?h%lqR`tr&lh%ofz>5GxrkaVhqTLjK7vGX^OXg^%UepDi&hRtN&eamX`5rY;1>6}CG9xfk zDV!g!?U_6nQiLzysT|hM0vFPcYqQa1U!cv4;J%-xx)&+|v|Jn5S_>HP49n0OTUo9v zIjP226H+ZexA*FYEJ%3{184NDzK7imiuh_^tqY`BVdM}km{*C~Ccgrv7jHKl)DZB7 zEY>ZAAKXmDP2C=aE?A=)8mI%O45|B(d_i4`)e(-9}pjmC^*M!x za%m>^o4wcMI9mSiYx#HR_T8YO{!FP^mtm!C54WHfzmA&N0@YH0q@l(q?Wmzn%s|7z z^8&q(tuQ;x8u^L^A|QnflN7*IP5NTR-CFYIn~MD9hN_E59IOkZfh#CSDtrzTD4ztdgevS#MC z_fHqYp2B|&9BZ_{e|!HBOFf3|iYC!D+Ia0Y>geMctJ5)tB7lV%^LB$$Ng$jj45^Am ze=LGGFJa+NUn89;=V{ID>PPNL=Hj(bMMmn=3(l4_o&G<$_a7}`C1^MP5bZVh*84sz zc5rpaS1-UoJQBw6x~*<)fP2<6pbzdn360A~n&-p2Yd)_CS`l1m&(&#zy1Cg{8!qG&TZj%xIJ?)qS_?|o3g z)S2m+BIH~?o_|Y2h+jo^I;I$j%g4{>gdW01ja*2SFG9N}(9xw{F%v6#zB615Yt7ib z0iwZ{h7xM%_X)psN?hn6G1Y4Hs6QSb9a=`1iR{zw!FZ4^^ZUE)lkAqxpFJ6b426srPK4e+?SuT%p(vs8y(WV0i<<~F`$-(@Cc*;6T8w9{e?*?GdznWkemPG2xu6hD#R@wtA) zJ^>4&P2owJw$n8V#t1GTZx9*0rE!BRCJ*LEi=D|Sf`ZRb*{m?h!bGzVkc9NXCw+oy zKx?2)Y6TWxZ#AGst_wEnS>YhwyNHe0Mc9RDsdi@qXAC-aLh~xE#m^?S5StjiVmyZ#DayU!~ zZeE2jRW_VtCsn+BSEe-GNQNRODI$+;IPoR`&+}&S29KSX(1@P#Cu({Ux}7IxIY-Wp zw$zna6_VaaC%SD=0X?kbs-|3YSVr@Utjxx2v&uVIX zBMQ;!d#5hLS&hMr>TMFG0kbt@zTf`LD4+*&om8@esxIx9U+ET8nv^iokb=qw6*4P@ z))#>4>|d$Fwm%^(Q`3RHUL0Cn%DJjTra-0Hot&CZg$ z8d*yX8p2bTb#hA2g}d$Z|KpdxYem$5e`4|qQEZ=4sJj(rv|J<6Q^v*zzF$@X>tE3Q zl6v#+^X^x8E9C(Et)F!ysaAhc7K$GzM?8$)TAO(FGXe5-dT;t7bWE4@7jM{>NA7T z2Pjxfq2M3ueY2M<(XDDG!hvVGHtrCn5sPH+ZRrf&pFMf#&iHBDS2DA(&p7m;qsEW( z_fZIeiPcH!5Y3kb9J~A5^^l0Ch<%B3~}j3Jh?moA9F_3u%G$V`%n4S4Im<=2J1Fn zsGe~HxN4Vxb&FNFO*6)LxoKZlnM|zJ(E+Sz2cRwGZ1Q;&Qs3*d*M?QpBzq%tHqgo$ zRp`}r2SoM`npX`$d{{-`BX{t>b~mUt zOKsuzHBeCM8JY9$3TFfHP3SGbZg;sd*Oe(Xy>@-CjKcI}_|TwX2(g$(cNFAvt~blNW~dNhx4o1<+94|G&Nlq z3preK|Jvg=+A)qStSR+Zj>jkHZn~HYgia0(D*oV28uemgWCmLmWu&?zb`Tjc$yHxQ zd~~cyTzssU>oL}&*X`~mj&mec8z|)3{x~@^uU5ZnCF%mU$u^xKYr4up|G1Ya5v_B- z%u>^pB_NAhX%C!vxPn)Y-AICda`^N)AOG*811ncaLl>(s=FE_&cWCOWNgOi2j3A>2 zIgJQAO4|xqbTyD6!dFqi8l#<5jvWK$4BkL+^}F4VHQ0Z*rr$BMLs(;WlJc;+w9cng z=x?&-sg#3Bjl@sb11m1Kl!EXN!JY?}KgX)Q}Q}^qxP#S{u5X% z2{;9I0e{>*&E^2s21Frmu~~ErupPoPOpy}-Mvehu&>XA@e6v!hkIC9u!UclU0-kfE zVdt(nx#Ck|j^vf6n>r>|Vc3NY2sZ1kEYNaKWl^V{6E6GGPNiyHBWx)D2%_Z&-#v<{ zbbwgh$SqE!4J8Is!c%ub3x)qVxvPsA^4f}F^dr^&5{M_z%DHaB+V`65xHX=K=jbo4 z@07nZG!U{qB_S_RWc%V8;Y1}t^^r@_#QX~U|9;_@ROM8FNe{9=d%9U4p0icrqIdvr zT#-(b!s|RImX-iywK)ilb)b*DT$PHq8WbIJ}f-ZD+;iS zyVYYPpSR<@Bo^+(TWkW zE_jT6P0?&)HW!Ea$cl=r!M);_u2Q=Au#$VaUZ()c=mV!#(j$O(A^e@}!Z0{hg015o zreaH+Xx6Lq3ECPFjJF{L6rXJ#6^*Cxxg5-@zEp#vsOZyOg`%!V{Qt(wK?i8MD(Yp+ zLMMOKBTJ=|vxhroFTGsryD#^z(pD)#fR_Zx^1-rP%X55K`wg`G^p-m;NNc qTq=fHI=CQ?pH!0l-|mYfV15u4OHRN!1emXcLFUHR7v347BK{9L$P6U_ literal 0 HcmV?d00001 diff --git a/doc/news/changes.h b/doc/news/changes.h index 2db0318d82..2253787358 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -66,6 +66,11 @@ inconvenience this causes.
    +
  1. New: FE_Q_Bubbles describes a FiniteElement based on FE_Q + enriched by bubble functions. +
    + (Daniel Arndt, 2015/08/12) +
  2. New: MultithreadInfo::set_thread_limit() can now be called more than once and the environment variable DEAL_II_NUM_THREADS will be respected diff --git a/include/deal.II/fe/fe_poly.templates.h b/include/deal.II/fe/fe_poly.templates.h index d93c690e58..ce2a7bf4d6 100644 --- a/include/deal.II/fe/fe_poly.templates.h +++ b/include/deal.II/fe/fe_poly.templates.h @@ -18,6 +18,7 @@ #include #include #include +#include #include #include diff --git a/include/deal.II/fe/fe_q_base.h b/include/deal.II/fe/fe_q_base.h index b89de91721..a5a9f83377 100644 --- a/include/deal.II/fe/fe_q_base.h +++ b/include/deal.II/fe/fe_q_base.h @@ -27,8 +27,8 @@ DEAL_II_NAMESPACE_OPEN /*@{*/ /** - * This class collects the basic methods used in FE_Q and FE_Q_DG0. There is - * no public constructor for this class as it is not functional as a stand- + * This class collects the basic methods used in FE_Q, FE_Q_DG0 and FE_Q_Bubbles. + * There is no public constructor for this class as it is not functional as a stand- * alone. The completion of definitions is left to the derived classes. * * @author Wolfgang Bangerth, 1998, 2003; Guido Kanschat, 2001; Ralf Hartmann, @@ -324,6 +324,13 @@ private: * Mutex for protecting initialization of restriction and embedding matrix. */ mutable Threads::Mutex mutex; + + /* + * The highest polynomial degree of the underlying tensor product space + * without any enrichment. For FE_Q*(p) this is p. Note that enrichments + * may lead to a difference to degree. + */ + const unsigned int q_degree; }; diff --git a/include/deal.II/fe/fe_q_bubbles.h b/include/deal.II/fe/fe_q_bubbles.h new file mode 100644 index 0000000000..1e6b9cd015 --- /dev/null +++ b/include/deal.II/fe/fe_q_bubbles.h @@ -0,0 +1,205 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2015 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#ifndef __deal2__fe_q_bubbles_h +#define __deal2__fe_q_bubbles_h + +#include +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +/*!@addtogroup fe */ +/*@{*/ + +/** + * Implementation of a scalar Lagrange finite element @p Q_p^+ that yields the + * finite element space of continuous, piecewise polynomials of degree @p p in + * each coordinate direction plus some bubble enrichment space spanned by + * $(2x_j-1)^{p-1}\prod_{i=0}^{dim-1}(x_i(1-x_i))$. Therefore the highest polynomial + * degree is $p+1$. + * This class is realized using tensor product polynomials based on equidistant + * or given support points. + * + * The standard constructor of this class takes the degree @p p of this finite + * element. Alternatively, it can take a quadrature formula @p points defining + * the support points of the Lagrange interpolation in one coordinate direction. + * + * For more information about the spacedim template parameter check + * the documentation of FiniteElement or the one of Triangulation. + * + * Due to the fact that the enrichments are small almost everywhere + * for large p, the condition number for the mass and stiffness matrix fastly + * increaseses with increasing p. + * Below you see a comparison with FE_Q(QGaussLobatto(p+1)) for dim=1. + * + *

    + * @image html fe_q_bubbles_conditioning.png + *

    + * + * Therefore, this element should be used with care for $p>3$. + * + *

    Implementation

    + * + * The constructor creates a TensorProductPolynomials object that includes the + * tensor product of @p LagrangeEquidistant polynomials of degree @p p plus the + * bubble enrichments. This @p TensorProductPolynomialsBubbles object + * provides all values and derivatives of the shape functions. In case a + * quadrature rule is given, the constructor creates a + * TensorProductPolynomialsBubbles object that includes the tensor product of @p + * Lagrange polynomials with the support points from @p points and the bubble enrichments + * as defined above. + * + * Furthermore the constructor fills the @p interface_constrains, the @p + * prolongation (embedding) and the @p restriction matrices. + * + *

    Numbering of the degrees of freedom (DoFs)

    + * + * The original ordering of the shape functions represented by the + * TensorProductPolynomialsBubbles is a tensor product + * numbering. However, the shape functions on a cell are renumbered + * beginning with the shape functions whose support points are at the + * vertices, then on the line, on the quads, and finally (for 3d) on + * the hexes. Finally, there are support points for the bubble enrichments + * in the middle of the cell. + * + */ +template +class FE_Q_Bubbles : public FE_Q_Base,dim,spacedim> +{ +public: + /** + * Constructor for tensor product polynomials of degree @p p plus bubble enrichments + * + */ + FE_Q_Bubbles (const unsigned int p); + + /** + * Constructor for tensor product polynomials with support points @p points + * plus bubble enrichments based on a one-dimensional quadrature + * formula. The degree of the finite element is points.size(). + * Note that the first point has to be 0 and the last one 1. + */ + FE_Q_Bubbles (const Quadrature<1> &points); + + /** + * Return a string that uniquely identifies a finite element. This class + * returns FE_Q_Bubbles(degree), with @p dim and @p degree + * replaced by appropriate values. + */ + virtual std::string get_name () const; + + /** + * Interpolate a set of scalar values, computed in the generalized support + * points. + */ + virtual void interpolate(std::vector &local_dofs, + const std::vector &values) const; + + /** + * Interpolate a set of vector values, computed in the generalized support + * points. + * + * Since a finite element often only interpolates part of a vector, + * offset is used to determine the first component of the vector to + * be interpolated. Maybe consider changing your data structures to use the + * next function. + */ + virtual void interpolate(std::vector &local_dofs, + const std::vector > &values, + unsigned int offset = 0) const; + + /** + * Interpolate a set of vector values, computed in the generalized support + * points. + */ + virtual void interpolate( + std::vector &local_dofs, + const VectorSlice > > &values) const; + + /** + * Return the matrix interpolating from the given finite element to the + * present one. The size of the matrix is then @p dofs_per_cell times + * source.dofs_per_cell. + * + * These matrices are only available if the source element is also a @p + * FE_Q_Bubbles element. Otherwise, an exception of type + * FiniteElement::ExcInterpolationNotImplemented is thrown. + */ + virtual void + get_interpolation_matrix (const FiniteElement &source, + FullMatrix &matrix) const; + + virtual const FullMatrix & + get_prolongation_matrix (const unsigned int child, + const RefinementCase &refinement_case) const; + + virtual const FullMatrix & + get_restriction_matrix (const unsigned int child, + const RefinementCase &refinement_case) const; + + /** + * Check for non-zero values on a face. + * + * This function returns @p true, if the shape function @p shape_index has + * non-zero values on the face @p face_index. + * + * Implementation of the interface in FiniteElement + */ + virtual bool has_support_on_face (const unsigned int shape_index, + const unsigned int face_index) const; + +protected: + /** + * @p clone function instead of a copy constructor. + * + * This function is needed by the constructors of @p FESystem. + */ + virtual FiniteElement *clone() const; + +private: + + /** + * Returns the restriction_is_additive flags. + * Only the last components for the bubble enrichments are true. + */ + static std::vector get_riaf_vector(const unsigned int degree); + + /** + * Only for internal use. Its full name is @p get_dofs_per_object_vector + * function and it creates the @p dofs_per_object vector that is needed + * within the constructor to be passed to the constructor of @p + * FiniteElementData. + */ + static std::vector get_dpo_vector(const unsigned int degree); + + /** + * Number of additional bubble functions + */ + const unsigned int n_bubbles; +}; + + + +/*@}*/ + + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/fe/CMakeLists.txt b/source/fe/CMakeLists.txt index b4130f5407..73fee37d1e 100644 --- a/source/fe/CMakeLists.txt +++ b/source/fe/CMakeLists.txt @@ -35,6 +35,7 @@ SET(_src fe_poly_tensor.cc fe_q_base.cc fe_q.cc + fe_q_bubbles.cc fe_q_dg0.cc fe_q_hierarchical.cc fe_q_iso_q1.cc @@ -72,6 +73,7 @@ SET(_inst fe_poly.inst.in fe_poly_tensor.inst.in fe_q_base.inst.in + fe_q_bubbles.inst.in fe_q_dg0.inst.in fe_q_hierarchical.inst.in fe_q.inst.in diff --git a/source/fe/fe_poly.cc b/source/fe/fe_poly.cc index 147b4eb984..ebd25e7906 100644 --- a/source/fe/fe_poly.cc +++ b/source/fe/fe_poly.cc @@ -17,6 +17,7 @@ #include #include #include +#include #include #include #include diff --git a/source/fe/fe_poly.inst.in b/source/fe/fe_poly.inst.in index 6331952dd9..dcd33ed098 100644 --- a/source/fe/fe_poly.inst.in +++ b/source/fe/fe_poly.inst.in @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 1998 - 2014 by the deal.II authors +// Copyright (C) 1998 - 2015 by the deal.II authors // // This file is part of the deal.II library. // @@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension <= deal_II_space_dimension template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; + template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly >, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; template class FE_Poly, deal_II_dimension, deal_II_space_dimension>; diff --git a/source/fe/fe_q_base.cc b/source/fe/fe_q_base.cc index 866290fac4..f6eb9fdbdc 100644 --- a/source/fe/fe_q_base.cc +++ b/source/fe/fe_q_base.cc @@ -19,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -124,6 +125,10 @@ struct FE_Q_Base::Implementation { const unsigned int dim = 2; + unsigned int q_deg = fe.degree; + if (types_are_equal >::value) + q_deg = fe.degree-1; + // restricted to each face, the traces of the shape functions is an // element of P_{k} (in 2d), or Q_{k} (in 3d), where k is the degree of // the element. from this, we interpolate between mother and cell face. @@ -172,10 +177,10 @@ struct FE_Q_Base::Implementation // Add midpoint constraint_points.push_back (Point (0.5)); - if (fe.degree>1) + if (q_deg>1) { - const unsigned int n=fe.degree-1; - const double step=1./fe.degree; + const unsigned int n=q_deg-1; + const double step=1./q_deg; // subface 0 for (unsigned int i=1; i<=n; ++i) constraint_points.push_back ( @@ -197,13 +202,13 @@ struct FE_Q_Base::Implementation const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(fe.degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); for (unsigned int i=0; i p; p[0] = constraint_points[i](0); @@ -227,6 +232,10 @@ struct FE_Q_Base::Implementation { const unsigned int dim = 3; + unsigned int q_deg = fe.degree; + if (types_are_equal >::value) + q_deg = fe.degree-1; + // For a detailed documentation of the interpolation see the // FE_Q_Base<2>::initialize_constraints function. @@ -253,10 +262,10 @@ struct FE_Q_Base::Implementation constraint_points.push_back (Point (0.5, 0)); constraint_points.push_back (Point (0.5, 1)); - if (fe.degree>1) + if (q_deg>1) { - const unsigned int n=fe.degree-1; - const double step=1./fe.degree; + const unsigned int n=q_deg-1; + const double step=1./q_deg; std::vector > line_support_points(n); for (unsigned int i=0; i::Implementation // Now construct relation between destination (child) and source (mother) // dofs. - const unsigned int pnts=(fe.degree+1)*(fe.degree+1); + const unsigned int pnts=(q_deg+1)*(q_deg+1); // use that the element evaluates to 1 at index 0 and along the line at // zero const std::vector &index_map_inverse = fe.poly_space.get_numbering_inverse(); const std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(fe.degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_deg); Assert(std::abs(fe.poly_space.compute_value(index_map_inverse[0],Point()) - 1.) < 1e-14, ExcInternalError()); @@ -327,7 +336,7 @@ struct FE_Q_Base::Implementation for (unsigned int i=0; i constraint_point; @@ -373,14 +382,14 @@ struct FE_Q_Base::Implementation for (unsigned int j=0; j::FE_Q_Base (const POLY &poly_space, const std::vector &restriction_is_additive_flags) : FE_Poly(poly_space, fe_data, restriction_is_additive_flags, - std::vector(1, std::vector(1,true))) + std::vector(1, std::vector(1,true))), + q_degree (types_are_equal >::value + ?this->degree-1 + :this->degree) {} @@ -422,17 +434,18 @@ FE_Q_Base::initialize (const std::vector > &points) // distinguish q/q_dg0 case: need to be flexible enough to allow more // degrees of freedom than there are FE_Q degrees of freedom for derived // class FE_Q_DG0 that otherwise shares 95% of the code. - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); Assert(q_dofs_per_cell == this->dofs_per_cell || - q_dofs_per_cell+1 == this->dofs_per_cell, ExcInternalError()); + q_dofs_per_cell+1 == this->dofs_per_cell || + q_dofs_per_cell+dim == this->dofs_per_cell , ExcInternalError()); { std::vector renumber(q_dofs_per_cell); - const FiniteElementData fe(get_dpo_vector(this->degree),1, - this->degree); + const FiniteElementData fe(get_dpo_vector(q_degree),1, + q_degree); FETools::hierarchic_to_lexicographic_numbering (fe, renumber); - if (this->dofs_per_cell > q_dofs_per_cell) - renumber.push_back(q_dofs_per_cell); + for (unsigned int i= q_dofs_per_cell; idofs_per_cell; ++i) + renumber.push_back(i); this->poly_space.set_numbering(renumber); } @@ -471,7 +484,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, x_source_fe.dofs_per_cell)); // only evaluate Q dofs - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); const unsigned int source_q_dofs_per_cell = Utilities::fixed_power(source_fe->degree+1); // evaluation is simply done by evaluating the other FE's basis functions on @@ -504,7 +517,7 @@ get_interpolation_matrix (const FiniteElement &x_source_fe, } // cut off very small values - const double eps = 2e-13*this->degree*dim; + const double eps = 2e-13*q_degree*dim; for (unsigned int i=0; idofs_per_cell; ++i) for (unsigned int j=0; jdofs_per_cell; ++j) if (std::fabs(interpolation_matrix(i,j)) < eps) @@ -596,7 +609,7 @@ get_subface_interpolation_matrix (const FiniteElement &x_source_fe // Rule of thumb for FP accuracy, that can be expected for a given // polynomial degree. This value is used to cut off values close to // zero. - double eps = 2e-13*this->degree*(dim-1); + double eps = 2e-13*q_degree*(dim-1); // compute the interpolation matrix by simply taking the value at the // support points. @@ -888,11 +901,11 @@ void FE_Q_Base::initialize_unit_face_support_points (const st return; const unsigned int codim = dim-1; - this->unit_face_support_points.resize(Utilities::fixed_power(this->degree+1)); + this->unit_face_support_points.resize(Utilities::fixed_power(q_degree+1)); // find renumbering of faces and assign from values of quadrature std::vector face_index_map = - FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(this->degree); + FE_Q_Helper::face_lexicographic_to_hierarchic_numbering(q_degree); Quadrature<1> support_1d(points); Quadrature support_quadrature(support_1d); this->unit_face_support_points.resize(support_quadrature.size()); @@ -914,7 +927,7 @@ FE_Q_Base::initialize_quad_dof_index_permutation () Assert (this->adjust_quad_dof_index_for_face_orientation_table.n_elements()==8*this->dofs_per_quad, ExcInternalError()); - const unsigned int n=this->degree-1; + const unsigned int n=q_degree-1; Assert(n*n==this->dofs_per_quad, ExcInternalError()); // alias for the table to fill @@ -1137,7 +1150,7 @@ FE_Q_Base return this->prolongation[refinement_case-1][child]; // distinguish q/q_dg0 case: only treat Q dofs first - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); // compute the interpolation matrices in much the same way as we do for // the constraints. it's actually simpler here, since we don't have this @@ -1145,7 +1158,7 @@ FE_Q_Base // interpolation matrix is formed by a permutation of the indices of the // cell matrix. The value eps is used a threshold to decide when certain // evaluations of the Lagrange polynomials are zero or one. - const double eps = 1e-15*this->degree*dim; + const double eps = 1e-15*q_degree*dim; #ifdef DEBUG // in DEBUG mode, check that the evaluation of support points in the @@ -1167,7 +1180,7 @@ FE_Q_Base // the tensor product structure of this element and only evaluate 1D // information from the polynomial. This makes the cost of this function // almost negligible also for high order elements - const unsigned int dofs1d = this->degree+1; + const unsigned int dofs1d = q_degree+1; std::vector > subcell_evaluations (dim, Table<2,double>(dofs1d, dofs1d)); const std::vector &index_map_inverse = @@ -1315,7 +1328,7 @@ FE_Q_Base FullMatrix restriction(this->dofs_per_cell, this->dofs_per_cell); // distinguish q/q_dg0 case - const unsigned int q_dofs_per_cell = Utilities::fixed_power(this->degree+1); + const unsigned int q_dofs_per_cell = Utilities::fixed_power(q_degree+1); // for these Lagrange interpolation polynomials, construction of the // restriction matrices is relatively simple. the reason is that the @@ -1339,11 +1352,11 @@ FE_Q_Base // (compute on one child) by the same value (compute on a later child), // so we don't have to care about this - const double eps = 1e-15*this->degree*dim; + const double eps = 1e-15*q_degree*dim; const std::vector &index_map_inverse = this->poly_space.get_numbering_inverse(); - const unsigned int dofs1d = this->degree+1; + const unsigned int dofs1d = q_degree+1; std::vector > evaluations1d (dofs1d); restriction.reinit(this->dofs_per_cell, this->dofs_per_cell); @@ -1531,9 +1544,12 @@ std::pair, std::vector > FE_Q_Base::get_constant_modes () const { Table<2,bool> constant_modes(1, this->dofs_per_cell); - // FE_Q_DG0 should not use this function... - AssertDimension(this->dofs_per_cell, Utilities::fixed_power(this->degree+1)); - constant_modes.fill(true); + // We here just care for the constant mode due to the polynomial space + // without any enrichments + // As there may be more constant modes derived classes may to implement this + // themselves. An example for this is FE_Q_DG0. + for (unsigned int i=0; i(q_degree+1); ++i) + constant_modes(0, i) = true; return std::pair, std::vector > (constant_modes, std::vector(1, 0)); } diff --git a/source/fe/fe_q_base.inst.in b/source/fe/fe_q_base.inst.in index f0b13ce5c1..1e36e6bee3 100644 --- a/source/fe/fe_q_base.inst.in +++ b/source/fe/fe_q_base.inst.in @@ -20,6 +20,7 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS #if deal_II_dimension <= deal_II_space_dimension template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; + template class FE_Q_Base, deal_II_dimension, deal_II_space_dimension>; template class FE_Q_Base >, deal_II_dimension, deal_II_space_dimension>; #endif } diff --git a/source/fe/fe_q_bubbles.cc b/source/fe/fe_q_bubbles.cc new file mode 100644 index 0000000000..7616ac57a3 --- /dev/null +++ b/source/fe/fe_q_bubbles.cc @@ -0,0 +1,527 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 2012 - 2013 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + + +#include +#include + +DEAL_II_NAMESPACE_OPEN + + +namespace FE_Q_Bubbles_Helper +{ + namespace + { + template + inline + void + compute_embedding_matrices(const FE_Q_Bubbles &fe, + std::vector > > &matrices, + const bool isotropic_only) + { + const unsigned int dpc = fe.dofs_per_cell; + const unsigned int degree = fe.degree; + + // Initialize quadrature formula on fine cells + MappingQ1 mapping; + + Quadrature *q_fine; + Quadrature<1> q_dummy(std::vector >(1), std::vector (1,1.)); + switch (dim) + { + case 1: + if (spacedim==1) + q_fine = new QGauss (degree+1); + else if (spacedim==2) + q_fine = new QAnisotropic(QGauss<1>(degree+1), q_dummy); + else + q_fine = new QAnisotropic(QGauss<1>(degree+1), q_dummy, q_dummy); + break; + case 2: + if (spacedim==2) + q_fine = new QGauss (degree+1); + else + q_fine = new QAnisotropic(QGauss<1>(degree+1), QGauss<1>(degree+1), q_dummy); + break; + case 3: + q_fine = new QGauss (degree+1); + break; + default: + AssertThrow(false, ExcInternalError()); + } + + const unsigned int nq = q_fine->size(); + + // loop over all possible refinement cases + unsigned int ref_case = (isotropic_only) + ? RefinementCase::isotropic_refinement + : RefinementCase::cut_x; + for (; ref_case <= RefinementCase::isotropic_refinement; ++ref_case) + { + const unsigned int nc + = GeometryInfo::n_children(RefinementCase(ref_case)); + + for (unsigned int i=0; i tr; + GridGenerator::hyper_cube (tr, 0, 1); + tr.begin_active()->set_refine_flag(RefinementCase(ref_case)); + tr.execute_coarsening_and_refinement(); + + DoFHandler dh(tr); + dh.distribute_dofs(fe); + + FEValues fine (mapping, fe, *q_fine, + update_quadrature_points + | update_JxW_values | update_values); + + const unsigned int n_dofs = dh.n_dofs(); + + FullMatrix fine_mass(n_dofs); + FullMatrix coarse_rhs_matrix(n_dofs, dpc); + + std::vector > child_ldi + (nc, std::vector(fe.dofs_per_cell)); + + //now create the mass matrix and all the right_hand sides + unsigned int child_no = 0; + typename dealii::DoFHandler::active_cell_iterator cell + = dh.begin_active(); + for (; cell!=dh.end(); ++cell, ++child_no) + { + fine.reinit(cell); + cell->get_dof_indices(child_ldi[child_no]); + + for (unsigned int q=0; q quad_tmp; + for (unsigned int k=0; k solution (n_dofs, dpc); + fine_mass.gauss_jordan(); + fine_mass.mmult(solution, coarse_rhs_matrix); + + //and distribute to the fine cell matrices + for (unsigned int child_no=0; child_no 1.e-12) + matrices[ref_case-1][child_no](i,j)=solution(gdi, j); + } + } + } + } +} + + +template +FE_Q_Bubbles::FE_Q_Bubbles (const unsigned int q_degree) + : + FE_Q_Base, dim, spacedim> ( + TensorProductPolynomialsBubbles(Polynomials::LagrangeEquidistant::generate_complete_basis(q_degree)), + FiniteElementData(get_dpo_vector(q_degree), + 1, q_degree+1, + FiniteElementData::H1), + get_riaf_vector(q_degree)), + n_bubbles((q_degree<=1)?1:dim) +{ + Assert (q_degree > 0, + ExcMessage ("This element can only be used for polynomial degrees " + "greater than zero")); + + std::vector > support_points_1d(q_degree+1); + for (unsigned int i=0; i<=q_degree; ++i) + support_points_1d[i][0] = static_cast(i)/q_degree; + + this->initialize(support_points_1d); + + // adjust unit support point for discontinuous node + Point point; + for (unsigned int d=0; dunit_support_points.push_back(point); + AssertDimension(this->dofs_per_cell, this->unit_support_points.size()); + + this->reinit_restriction_and_prolongation_matrices(); + if (dim == spacedim) + { + FE_Q_Bubbles_Helper::compute_embedding_matrices + (*this, this->prolongation, false); + // Fill restriction matrices with L2-projection + FETools::compute_projection_matrices (*this, this->restriction); + } + +} + + + +template +FE_Q_Bubbles::FE_Q_Bubbles (const Quadrature<1> &points) + : + FE_Q_Base, dim, spacedim> ( + TensorProductPolynomialsBubbles(Polynomials::generate_complete_Lagrange_basis(points.get_points())), + FiniteElementData(get_dpo_vector(points.size()-1), + 1, points.size(), + FiniteElementData::H1), + get_riaf_vector(points.size()-1)), + n_bubbles((points.size()-1<=1)?1:dim) +{ + const unsigned int q_degree = points.size()-1; + (void) q_degree; + Assert (q_degree > 0, + ExcMessage ("This element can only be used for polynomial degrees " + "at least one")); + + this->initialize(points.get_points()); + + // adjust unit support point for discontinuous node + Point point; + for (unsigned int d=0; dunit_support_points.push_back(point); + AssertDimension(this->dofs_per_cell, this->unit_support_points.size()); + + this->reinit_restriction_and_prolongation_matrices(); + if (dim == spacedim) + { + FE_Q_Bubbles_Helper::compute_embedding_matrices + (*this, this->prolongation, false); + // Fill restriction matrices with L2-projection + FETools::compute_projection_matrices (*this, this->restriction); + } +} + + + +template +std::string +FE_Q_Bubbles::get_name () const +{ + // note that the FETools::get_fe_from_name function depends on the + // particular format of the string this function returns, so they have to be + // kept in synch + + std::ostringstream namebuf; + bool type = true; + const unsigned int n_points = this->degree; + std::vector points(n_points); + const unsigned int dofs_per_cell = this->dofs_per_cell; + const std::vector > &unit_support_points = this->unit_support_points; + unsigned int index = 0; + + // Decode the support points in one coordinate direction. + for (unsigned int j=0; j1) ? (unit_support_points[j](1)==0 && + ((dim>2) ? unit_support_points[j](2)==0: true)) : true) + { + if (index == 0) + points[index] = unit_support_points[j](0); + else if (index == 1) + points[n_points-1] = unit_support_points[j](0); + else + points[index-1] = unit_support_points[j](0); + + index++; + } + } + // Do not consider the discontinuous node for dimension 1 + Assert (index == n_points || (dim==1 && index == n_points+n_bubbles), + ExcMessage ("Could not decode support points in one coordinate direction.")); + + // Check whether the support points are equidistant. + for (unsigned int j=0; jdegree-1)) > 1e-15) + { + type = false; + break; + } + + if (type == true) + namebuf << "FE_Q_Bubbles<" << dim << ">(" << this->degree-1 << ")"; + else + { + + // Check whether the support points come from QGaussLobatto. + const QGaussLobatto<1> points_gl(n_points); + type = true; + for (unsigned int j=0; j(QGaussLobatto(" << this->degree << "))"; + else + namebuf << "FE_Q_Bubbles<" << dim << ">(QUnknownNodes(" << this->degree-1 << "))"; + } + return namebuf.str(); +} + + + +template +FiniteElement * +FE_Q_Bubbles::clone() const +{ + return new FE_Q_Bubbles(*this); +} + + + +template +void +FE_Q_Bubbles::interpolate(std::vector &local_dofs, + const std::vector &values) const +{ + Assert (values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (this->n_components() == 1, + ExcDimensionMismatch(this->n_components(), 1)); + + std::copy(values.begin(), values.end(), local_dofs.begin()); + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles::interpolate(std::vector &local_dofs, + const std::vector > &values, + unsigned int offset) const +{ + Assert (values.size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (values[0].size() >= offset+this->n_components(), + ExcDimensionMismatch(values[0].size(),offset+this->n_components())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[i](offset+index.first); + } + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles::interpolate( + std::vector &local_dofs, + const VectorSlice > > &values) const +{ + Assert (values[0].size() == this->unit_support_points.size(), + ExcDimensionMismatch(values.size(), + this->unit_support_points.size())); + Assert (local_dofs.size() == this->dofs_per_cell, + ExcDimensionMismatch(local_dofs.size(),this->dofs_per_cell)); + Assert (values.size() == this->n_components(), + ExcDimensionMismatch(values.size(), this->n_components())); + + for (unsigned int i=0; idofs_per_cell-1; ++i) + { + const std::pair index + = this->system_to_component_index(i); + local_dofs[i] = values[index.first][i]; + } + + // We don't use the bubble functions for local interpolation + for (unsigned int i = 0; i +void +FE_Q_Bubbles:: +get_interpolation_matrix (const FiniteElement &x_source_fe, + FullMatrix &interpolation_matrix) const +{ + // We don't know how to do this properly, yet. + // However, for SolutionTransfer to work we need to provide an implementation + // for the case that the x_source_fe is identical to this FE + typedef FE_Q_Bubbles FEQBUBBLES; + typedef FiniteElement FEL; + + AssertThrow ((x_source_fe.get_name().find ("FE_Q_Bubbles<") == 0) + || + (dynamic_cast(&x_source_fe) != 0) + , + typename FEL:: + ExcInterpolationNotImplemented()); + Assert (interpolation_matrix.m() == this->dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + this->dofs_per_cell)); + Assert (interpolation_matrix.n() == x_source_fe.dofs_per_cell, + ExcDimensionMismatch (interpolation_matrix.m(), + x_source_fe.dofs_per_cell)); + + //Provide a short cut in case we are just inquiring the identity + if (dynamic_cast(&x_source_fe)->degree == this->degree) + for (unsigned int i=0; i +std::vector +FE_Q_Bubbles::get_riaf_vector(const unsigned int q_deg) +{ + unsigned int n_cont_dofs = Utilities::fixed_power(q_deg+1); + const unsigned int n_bubbles = (q_deg<=1?1:dim); + return std::vector (n_cont_dofs+n_bubbles,true); +} + + + +template +std::vector +FE_Q_Bubbles::get_dpo_vector(const unsigned int q_deg) +{ + std::vector dpo(dim+1, 1U); + for (unsigned int i=1; i +bool +FE_Q_Bubbles::has_support_on_face +(const unsigned int shape_index, + const unsigned int face_index) const +{ + // discontinuous functions have no support on faces + if (shape_index >= this->n_dofs_per_cell()-n_bubbles) + return false; + else + return FE_Q_Base,dim,spacedim>::has_support_on_face(shape_index, face_index); +} + + + +template +const FullMatrix & +FE_Q_Bubbles::get_prolongation_matrix +(const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Prolongation matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + Assert (this->prolongation[refinement_case-1][child].n() != 0, + ExcMessage("This prolongation matrix has not been computed yet!")); + // finally return the matrix + return this->prolongation[refinement_case-1][child]; +} + + + +template +const FullMatrix & +FE_Q_Bubbles::get_restriction_matrix +(const unsigned int child, + const RefinementCase &refinement_case) const +{ + Assert (refinement_case::isotropic_refinement+1, + ExcIndexRange(refinement_case,0,RefinementCase::isotropic_refinement+1)); + Assert (refinement_case!=RefinementCase::no_refinement, + ExcMessage("Restriction matrices are only available for refined cells!")); + Assert (child::n_children(refinement_case), + ExcIndexRange(child,0,GeometryInfo::n_children(refinement_case))); + + Assert(this->restriction[refinement_case-1][child].n() != 0, + ExcMessage("This restriction matrix has not been computed yet!")); + + //finally return the matrix + return this->restriction[refinement_case-1][child]; +} + +// explicit instantiations +#include "fe_q_bubbles.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/fe/fe_q_bubbles.inst.in b/source/fe/fe_q_bubbles.inst.in new file mode 100644 index 0000000000..7cf26de8ed --- /dev/null +++ b/source/fe/fe_q_bubbles.inst.in @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// $Id$ +// +// Copyright (C) 1998 - 2013 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_Q_Bubbles; +#endif + } -- 2.39.5