NGSolve
4.9
|
00001 #ifndef FILE_L2HOFE 00002 #define FILE_L2HOFE 00003 00004 /*********************************************************************/ 00005 /* File: l2hofe.hpp */ 00006 /* Author: Start */ 00007 /* Date: 6. Feb. 2003 */ 00008 /*********************************************************************/ 00009 00010 00011 #include "tscalarfe.hpp" 00012 #include "precomp.hpp" 00013 00014 00015 namespace ngfem 00016 { 00017 00022 template<int D> 00023 class L2HighOrderFiniteElement : virtual public ScalarFiniteElement<D> 00024 { 00025 protected: 00026 00027 enum { DIM = D }; 00028 int vnums[8]; 00029 INT<DIM> order_inner; 00030 00031 using ScalarFiniteElement<DIM>::ndof; 00032 using ScalarFiniteElement<DIM>::order; 00033 using ScalarFiniteElement<DIM>::eltype; 00034 00035 00036 public: 00038 template <typename TA> 00039 void SetVertexNumbers (const TA & avnums) 00040 { for (int i = 0; i < avnums.Size(); i++) vnums[i] = avnums[i]; } 00041 00042 void SetVertexNumber (int nr, int vnum) { vnums[nr] = vnum; } 00043 00045 // void SetOrder (int p) { order_inner = p; } 00046 00048 void SetOrder (INT<DIM> p) { order_inner = p; } 00049 00051 virtual void ComputeNDof () = 0; 00052 00053 virtual void GetInternalDofs (Array<int> & idofs) const; 00054 00055 virtual void PrecomputeTrace () = 0; 00056 00057 void CalcTraceMatrix (int facet, FlatMatrix<> & trace) const; 00058 00059 virtual void GetTrace (int facet, FlatVector<> coefs, FlatVector<> fcoefs) const 00060 { 00061 Matrix<> trace(fcoefs.Size(), coefs.Size()); 00062 CalcTraceMatrix(facet, trace); 00063 fcoefs = trace * coefs; 00064 } 00065 00066 virtual void GetTraceTrans (int facet, FlatVector<> fcoefs, FlatVector<> coefs) const 00067 { 00068 Matrix<> trace(fcoefs.Size(), coefs.Size()); 00069 CalcTraceMatrix(facet, trace); 00070 coefs = Trans (trace) * fcoefs; 00071 } 00072 }; 00073 00074 00075 00080 // template <ELEMENT_TYPE ET> class L2HighOrderFE; 00081 00082 00086 template <ELEMENT_TYPE ET> 00087 class T_L2HighOrderFiniteElement : 00088 public L2HighOrderFiniteElement<ET_trait<ET>::DIM>, 00089 public ET_trait<ET> 00090 { 00091 protected: 00092 enum { DIM = ET_trait<ET>::DIM }; 00093 00094 using ScalarFiniteElement<DIM>::ndof; 00095 using ScalarFiniteElement<DIM>::order; 00096 using ScalarFiniteElement<DIM>::eltype; 00097 00098 using L2HighOrderFiniteElement<DIM>::vnums; 00099 using L2HighOrderFiniteElement<DIM>::order_inner; 00100 00101 using ET_trait<ET>::N_VERTEX; 00102 using ET_trait<ET>::N_EDGE; 00103 using ET_trait<ET>::N_FACE; 00104 using ET_trait<ET>::FaceType; 00105 using ET_trait<ET>::GetEdgeSort; 00106 using ET_trait<ET>::GetFaceSort; 00107 using ET_trait<ET>::PolDimension; 00108 00109 00110 public: 00111 00112 T_L2HighOrderFiniteElement () 00113 { eltype = ET; } 00114 00115 T_L2HighOrderFiniteElement (int aorder) 00116 { 00117 for (int i = 0; i < ET_trait<ET>::N_VERTEX; i++) vnums[i] = i; 00118 eltype = ET; 00119 00120 order = aorder; 00121 order_inner = aorder; 00122 ndof = PolDimension (aorder); 00123 } 00124 00125 virtual void ComputeNDof(); 00126 }; 00127 00128 00129 00130 00131 00132 00133 template <ELEMENT_TYPE ET> class L2HighOrderFE_Shape; 00134 00135 template <int DIM> 00136 class PrecomputedScalShapes 00137 { 00138 public: 00139 Matrix<> shapes; 00140 Matrix<> dshapes; 00141 00142 PrecomputedScalShapes (int nip, int ndof) 00143 : shapes(nip, ndof), dshapes (DIM*nip, ndof) 00144 { ; } 00145 00146 }; 00147 00148 00149 template <ELEMENT_TYPE ET, template <ELEMENT_TYPE ET> class SHAPES = L2HighOrderFE_Shape> 00150 class NGS_DLL_HEADER L2HighOrderFE : public T_L2HighOrderFiniteElement<ET>, 00151 public T_ScalarFiniteElement2<SHAPES<ET>, ET > 00152 { 00153 protected: 00154 using T_L2HighOrderFiniteElement<ET>::ndof; 00155 using T_L2HighOrderFiniteElement<ET>::order; 00156 using T_L2HighOrderFiniteElement<ET>::vnums; 00157 00158 enum { DIM = ET_trait<ET>::DIM }; 00159 typedef PrecomputedShapesContainer<PrecomputedScalShapes<DIM> > TPRECOMP; 00160 static TPRECOMP precomp; 00161 00162 typedef HashTable<INT<2>, Matrix<>*> TPRECOMP_TRACE; 00163 static TPRECOMP_TRACE precomp_trace; 00164 public: 00165 L2HighOrderFE () { ; } 00166 00167 L2HighOrderFE (int aorder) 00168 : T_L2HighOrderFiniteElement<ET> (aorder) { ; } 00169 00170 00171 virtual void PrecomputeTrace (); 00172 virtual void PrecomputeShapes (const IntegrationRule & ir); 00173 00174 virtual void Evaluate (const IntegrationRule & ir, FlatVector<double> coefs, FlatVector<double> vals) const; 00175 00176 virtual void EvaluateGrad (const IntegrationRule & ir, FlatVector<> coefs, FlatMatrixFixWidth<DIM> values) const; 00177 00178 virtual void EvaluateGradTrans (const IntegrationRule & ir, FlatMatrixFixWidth<DIM> values, FlatVector<> coefs) const; 00179 00180 virtual void GetTrace (int facet, FlatVector<> coefs, FlatVector<> fcoefs) const; 00181 00182 virtual void GetTraceTrans (int facet, FlatVector<> fcoefs, FlatVector<> coefs) const; 00183 }; 00184 00185 00186 00187 00188 00189 00190 00191 00192 00193 00194 00198 template <> 00199 class L2HighOrderFE_Shape<ET_SEGM> : public L2HighOrderFE<ET_SEGM> 00200 { 00201 public: 00202 template<typename Tx, typename TFA> 00203 void T_CalcShape (Tx x[], TFA & shape) const 00204 { 00205 Tx lam[2] = { x[0], 1-x[0] }; 00206 INT<2> e = GetEdgeSort (0, vnums); 00207 LegendrePolynomial (order, lam[e[1]]-lam[e[0]], shape); 00208 } 00209 }; 00210 00211 00215 template <> 00216 class L2HighOrderFE_Shape<ET_TRIG> : public L2HighOrderFE<ET_TRIG> 00217 { 00218 public: 00219 template<typename Tx, typename TFA> 00220 void T_CalcShape (Tx x[], TFA & shape) const 00221 { 00222 Tx lam[3] = { x[0], x[1], 1-x[0]-x[1] }; 00223 INT<4> f = GetFaceSort (0, vnums); 00224 int p = order_inner[0]; 00225 DubinerBasis::Eval (p, lam[f[0]], lam[f[1]], shape); 00226 } 00227 }; 00228 00229 00233 template <> 00234 class L2HighOrderFE_Shape<ET_QUAD> : public L2HighOrderFE<ET_QUAD> 00235 { 00236 public: 00237 template<typename Tx, typename TFA> 00238 void T_CalcShape (Tx hx[], TFA & shape) const 00239 { 00240 Tx x = hx[0], y = hx[1]; 00241 00242 Tx sigma[4] = {(1-x)+(1-y),x+(1-y),x+y,(1-x)+y}; 00243 00244 INT<4> f = GetFaceSort (0, vnums); 00245 00246 Tx xi = sigma[f[0]]-sigma[f[1]]; 00247 Tx eta = sigma[f[0]]-sigma[f[3]]; 00248 00249 int n = max(order_inner[0],order_inner[1]); 00250 ArrayMem<Tx, 20> polx(n+1), poly(n+1); 00251 00252 LegendrePolynomial (n, xi, polx); 00253 LegendrePolynomial (n, eta, poly); 00254 00255 for (int i = 0, ii = 0; i <= order_inner[0]; i++) 00256 for (int j = 0; j <= order_inner[1]; j++) 00257 shape[ii++] = polx[i] * poly[j]; 00258 } 00259 }; 00260 00261 00265 template <> 00266 class L2HighOrderFE_Shape<ET_TET> : public L2HighOrderFE<ET_TET> 00267 { 00268 public: 00269 template<typename Tx, typename TFA> 00270 void T_CalcShape (Tx x[], TFA & shape) const 00271 { 00272 Tx lami[4] = { x[0], x[1], x[2], 1-x[0]-x[1]-x[2] }; 00273 00274 int sort[4]; 00275 for (int i = 0; i < 4; i++) sort[i] = i; 00276 00277 if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]); 00278 if (vnums[sort[2]] > vnums[sort[3]]) Swap (sort[2], sort[3]); 00279 if (vnums[sort[0]] > vnums[sort[2]]) Swap (sort[0], sort[2]); 00280 if (vnums[sort[1]] > vnums[sort[3]]) Swap (sort[1], sort[3]); 00281 if (vnums[sort[1]] > vnums[sort[2]]) Swap (sort[1], sort[2]); 00282 00283 Tx lamis[4]; 00284 for (int i = 0; i < 4; i++) 00285 lamis[i] = lami[sort[i]]; 00286 00287 ArrayMem<Tx, 20> memx(sqr(order+1)); 00288 ArrayMem<Tx, 20> memy(sqr(order+1)); 00289 00290 FlatMatrix<Tx> polsx(order+1, &memx[0]); 00291 FlatMatrix<Tx> polsy(order+1, &memy[0]); 00292 VectorMem<10, Tx> polsz(order+1); 00293 00294 for (int i = 0; i <= order; i++) 00295 JacobiPolynomial (order, 2*lamis[0]-1, 2*i+2, 0, polsx.Row(i)); 00296 for (int i = 0; i <= order; i++) 00297 ScaledJacobiPolynomial (order, lamis[1]-lamis[2]-lamis[3], 1-lamis[0], 2*i+1, 0, polsy.Row(i)); 00298 00299 ScaledLegendrePolynomial (order, lamis[2]-lamis[3], lamis[2]+lamis[3], polsz); 00300 00301 for (int i = 0, ii = 0; i <= order; i++) 00302 for (int j = 0; j <= order-i; j++) 00303 for (int k = 0; k <= order-i-j; k++, ii++) 00304 shape[ii] = polsz(k) * polsy(k, j) * polsx(j+k, i); 00305 } 00306 00307 }; 00308 00309 00313 template <> 00314 class L2HighOrderFE_Shape<ET_PRISM> : public L2HighOrderFE<ET_PRISM> 00315 { 00316 public: 00317 template<typename Tx, typename TFA> 00318 void T_CalcShape (Tx hx[], TFA & shape) const 00319 { 00320 Tx lami[3] = { hx[0], hx[1], 1-hx[0]-hx[1] }; 00321 00322 int sort[3]; 00323 for (int i = 0; i < 3; i++) sort[i] = i; 00324 00325 if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]); 00326 if (vnums[sort[1]] > vnums[sort[2]]) Swap (sort[1], sort[2]); 00327 if (vnums[sort[0]] > vnums[sort[1]]) Swap (sort[0], sort[1]); 00328 00329 Tx lamis[3]; 00330 for (int i = 0; i < 3; i++) 00331 lamis[i] = lami[sort[i]]; 00332 00333 Tx x = lamis[0]; 00334 // Tx y = lamis[1]; 00335 Tx z = hx[2]; 00336 00337 int p=order_inner[0]; 00338 int q=order_inner[1]; 00339 00340 ArrayMem<Tx, 20> memx(sqr(p+1)); 00341 FlatMatrix<Tx> polsx(p+1, &memx[0]); 00342 00343 VectorMem<10, Tx> polsy(p+1); 00344 VectorMem<10, Tx> polsz(q+1); 00345 00346 for (int i = 0; i <= p; i++) 00347 JacobiPolynomial (p, 2*x-1, 2*i+1, 0, polsx.Row(i)); 00348 00349 ScaledLegendrePolynomial (order, lamis[1]-lamis[2], lamis[1]+lamis[2], polsy); 00350 LegendrePolynomial (order, 2*z-1, polsz); 00351 00352 int ii = 0; 00353 for (int k = 0; k <= q; k++) 00354 for (int i = 0; i <= p; i++) 00355 for (int j = 0; j <= p-i; j++) 00356 shape[ii++] = polsx(j,i) * polsy(j) * polsz(k); 00357 } 00358 00359 }; 00360 00361 00365 template <> 00366 class L2HighOrderFE_Shape<ET_PYRAMID> : public L2HighOrderFE<ET_PYRAMID> 00367 { 00368 public: 00369 template<typename Tx, typename TFA> 00370 void T_CalcShape (Tx x[], TFA & shape) const; 00371 }; 00372 00373 00377 template <> 00378 class L2HighOrderFE_Shape<ET_HEX> : public L2HighOrderFE<ET_HEX> 00379 { 00380 public: 00381 template<typename Tx, typename TFA> 00382 void T_CalcShape (Tx x[], TFA & shape) const; 00383 }; 00384 } 00385 00386 #endif