From: Wolfgang Bangerth Date: Wed, 7 Mar 2001 13:57:11 +0000 (+0000) Subject: Rather pervasive clean-ups. Move some functions to private sections, fix two bugs... X-Git-Tag: v8.0.0~19627 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=93c8e9057fffb9c7568c71dcd699b2bc32d2b777;p=dealii.git Rather pervasive clean-ups. Move some functions to private sections, fix two bugs where compute_2nd was called multiply. git-svn-id: https://svn.dealii.org/trunk@4144 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/fe/fe_q.h b/deal.II/deal.II/include/fe/fe_q.h index 4314c96e47..90fd7ca6ce 100644 --- a/deal.II/deal.II/include/fe/fe_q.h +++ b/deal.II/deal.II/include/fe/fe_q.h @@ -19,8 +19,11 @@ template class TensorProductPolynomials; template class MappingQ; + + /** * Tensor product elements based on equidistant support points. +//TODO: Document node numbering etc. copy from old documentation */ template class FE_Q : public FiniteElement @@ -37,18 +40,6 @@ class FE_Q : public FiniteElement */ ~FE_Q (); - /** - * Compute flags for initial - * update only. - */ - virtual UpdateFlags update_once (UpdateFlags flags) const; - - /** - * Compute flags for update on - * each cell. - */ - virtual UpdateFlags update_each (UpdateFlags flags) const; - /** * Return the support points of * the trial functions on the @@ -114,8 +105,7 @@ class FE_Q : public FiniteElement */ virtual unsigned int memory_consumption () const; - protected: - + protected: /** * @p{clone} function instead of * a copy constructor. @@ -178,6 +168,74 @@ class FE_Q : public FiniteElement FEValuesData& data) const ; private: + + /** + * Declare a nested class which + * will has static definitions of + * various matrices such as + * constraint and embedding + * matrices. The definition of + * the various static fields are + * in the files @p{mat_q.[123]} + * in the source directory. + */ + struct Matrices { + /** + * Embedding matrices. For + * each element type (the + * first index) there are as + * many embedding matrices as + * there are children per + * cell. The first index + * starts with linear + * elements and goes up in + * polynomial degree. The + * array may grow in the + * future with the number of + * elements for which these + * matrices have been + * computed. If for some + * element, the matrices have + * not been computed then you + * may use the element + * nevertheless but can not + * access the respective + * fields. + */ + static const double * const + embedding[][GeometryInfo::children_per_cell]; + + /** + * Number of elements (first + * index) the above field + * has. Equals the highest + * polynomial degree for + * which the embedding + * matrices have been + * computed. + */ + static const unsigned int n_embedding_matrices; + + /** + * As the + * @p{embedding_matrices} + * field, but for the + * interface constraints. One + * for each element for which + * it has been computed. + */ + static const double * const constraint_matrices[]; + + /** + * Like + * @p{n_embedding_matrices}, + * but for the number of + * interface constraint + * matrices. + */ + static const unsigned int n_constraint_matrices; + }; + /** * Only for internal use. Its * full name is @@ -221,8 +279,8 @@ class FE_Q : public FiniteElement * @p{FiniteElementData} argument. */ static void build_renumbering (const FiniteElementData &fe_data, - unsigned int degree, - std::vector& numbering); + const unsigned int degree, + std::vector &numbering); /** * Map tensor product data to @@ -235,18 +293,34 @@ class FE_Q : public FiniteElement * * It is used to compute the * order of face support points. + * + * In 1d, this function does + * nothing. */ static void build_face_renumbering (const FiniteElementData &fe_data, - unsigned int degree, - std::vector& numbering); + const unsigned int degree, + std::vector &numbering); /** * Compute support points. */ - static void compute_support_points (std::vector >& support_points, - unsigned int degree, - const std::vector& renumber); + static void compute_support_points (const unsigned int degree, + const std::vector &renumber, + std::vector > &support_points); + + /** + * Compute flags for initial + * update only. + */ + virtual UpdateFlags update_once (UpdateFlags flags) const; + + /** + * Compute flags for update on + * each cell. + */ + virtual UpdateFlags update_each (UpdateFlags flags) const; + /** * Degree of the polynomials. */ diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index a1f585675a..29cb80fd3a 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -24,7 +24,8 @@ #include -// Embedding matrices (produced by tests/fe/embedding with postprocessing +// Include embedding matrices (produced by tests/fe/embedding with +// postprocessing #if (deal_II_dimension == 1) #include "mat_q.1" @@ -39,8 +40,10 @@ #endif +//TODO: move build_renumbering to mapping class + template -FE_Q::FE_Q (unsigned int degree) +FE_Q::FE_Q (const unsigned int degree) : FiniteElement (FiniteElementData(get_dpo_vector(degree),1), std::vector (1,false)), @@ -49,539 +52,537 @@ FE_Q::FE_Q (unsigned int degree) face_renumber(dofs_per_face, 0), poly(0) { + // Q0 elements cannot be + // continuous, use FE_DGQ(0) + // instead + Assert (degree>0, ExcNotImplemented()); std::vector v; for (unsigned int i=0;i<=degree;++i) v.push_back(LagrangeEquidistant(degree,i)); poly = new TensorProductPolynomials (v); - build_renumbering (*this, degree, renumber); -#if (deal_II_dimension > 1) + // do some internal book-keeping + build_renumbering (*this, degree, renumber); +#if deal_II_dimension > 1 build_face_renumbering (FiniteElementData(FE_Q::get_dpo_vector(degree),1), degree, face_renumber); #endif - Assert (degree <= 10, ExcNotImplemented()); - - if ((dim>1) && (constraint_matrices[degree-1] != 0)) - interface_constraints.fill (constraint_matrices[degree-1]); + // copy constraint matrices if they + // are defined. otherwise set them + // to invalid size + if ((dim>1) && (degree1) - { - prolongation[2].fill (embedding2[degree-1]); - prolongation[3].fill (embedding3[degree-1]); - } - if (dim>2) - { - prolongation[4].fill (embedding4[degree-1]); - prolongation[5].fill (embedding5[degree-1]); - prolongation[6].fill (embedding6[degree-1]); - prolongation[7].fill (embedding7[degree-1]); - } - } + + // next copy over embedding + // matrices if they are defined + if ((degree < Matrices::n_embedding_matrices+1) && + (Matrices::embedding[degree-1][0] != 0)) + for (unsigned int c=0; c::children_per_cell; ++c) + prolongation[c].fill (Matrices::embedding[degree-1][c]); else for (unsigned int i=0; i::children_per_cell;++i) prolongation[i].reinit(0,0); - + + // then fill restriction + // matrices. they are hardcoded for + // the first few elements switch (dim) { - case 1: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](2,1) = 1; - restriction[1](1,1) = 1; - restriction[1](1,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[1](1,1) = 1; - restriction[1](3,2) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](2,3) = 1; - restriction[0](3,1) = 1; - restriction[1](1,1) = 1; - restriction[1](3,0) = 1; - restriction[1](4,3) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; - case 2: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](4,1) = 1; - restriction[0](7,3) = 1; - restriction[0](8,2) = 1; - restriction[1](1,1) = 1; - restriction[1](4,0) = 1; - restriction[1](5,2) = 1; - restriction[1](8,3) = 1; - restriction[2](2,2) = 1; - restriction[2](5,1) = 1; - restriction[2](6,3) = 1; - restriction[2](8,0) = 1; - restriction[3](3,3) = 1; - restriction[3](6,2) = 1; - restriction[3](7,0) = 1; - restriction[3](8,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](10,11) = 1; - restriction[0](12,15) = 1; - restriction[1](1,1) = 1; - restriction[1](5,4) = 1; - restriction[1](6,7) = 1; - restriction[1](13,14) = 1; - restriction[2](2,2) = 1; - restriction[2](7,6) = 1; - restriction[2](9,8) = 1; - restriction[2](15,12) = 1; - restriction[3](3,3) = 1; - restriction[3](8,9) = 1; - restriction[3](11,10) = 1; - restriction[3](14,13) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](4,5) = 1; - restriction[0](5,1) = 1; - restriction[0](13,14) = 1; - restriction[0](14,3) = 1; - restriction[0](16,20) = 1; - restriction[0](17,8) = 1; - restriction[0](19,11) = 1; - restriction[0](20,2) = 1; - restriction[1](1,1) = 1; - restriction[1](5,0) = 1; - restriction[1](6,5) = 1; - restriction[1](7,8) = 1; - restriction[1](8,2) = 1; - restriction[1](17,14) = 1; - restriction[1](18,20) = 1; - restriction[1](20,3) = 1; - restriction[1](21,11) = 1; - restriction[2](2,2) = 1; - restriction[2](8,1) = 1; - restriction[2](9,8) = 1; - restriction[2](11,3) = 1; - restriction[2](12,11) = 1; - restriction[2](20,0) = 1; - restriction[2](21,5) = 1; - restriction[2](23,14) = 1; - restriction[2](24,20) = 1; - restriction[3](3,3) = 1; - restriction[3](10,11) = 1; - restriction[3](11,2) = 1; - restriction[3](14,0) = 1; - restriction[3](15,14) = 1; - restriction[3](19,5) = 1; - restriction[3](20,1) = 1; - restriction[3](22,20) = 1; - restriction[3](23,8) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; - case 3: - switch (degree) - { - case 1: - restriction[0](0,0) = 1; - restriction[1](1,1) = 1; - restriction[2](2,2) = 1; - restriction[3](3,3) = 1; - restriction[4](4,4) = 1; - restriction[5](5,5) = 1; - restriction[6](6,6) = 1; - restriction[7](7,7) = 1; - break; - case 2: - restriction[0](0,0) = 1; - restriction[0](8,1) = 1; - restriction[0](11,3) = 1; - restriction[0](16,4) = 1; - restriction[0](20,2) = 1; - restriction[0](22,5) = 1; - restriction[0](25,7) = 1; - restriction[0](26,6) = 1; - restriction[1](1,1) = 1; - restriction[1](8,0) = 1; - restriction[1](9,2) = 1; - restriction[1](17,5) = 1; - restriction[1](20,3) = 1; - restriction[1](22,4) = 1; - restriction[1](23,6) = 1; - restriction[1](26,7) = 1; - restriction[2](2,2) = 1; - restriction[2](9,1) = 1; - restriction[2](10,3) = 1; - restriction[2](18,6) = 1; - restriction[2](20,0) = 1; - restriction[2](23,5) = 1; - restriction[2](24,7) = 1; - restriction[2](26,4) = 1; - restriction[3](3,3) = 1; - restriction[3](10,2) = 1; - restriction[3](11,0) = 1; - restriction[3](19,7) = 1; - restriction[3](20,1) = 1; - restriction[3](24,6) = 1; - restriction[3](25,4) = 1; - restriction[3](26,5) = 1; - restriction[4](4,4) = 1; - restriction[4](12,5) = 1; - restriction[4](15,7) = 1; - restriction[4](16,0) = 1; - restriction[4](21,6) = 1; - restriction[4](22,1) = 1; - restriction[4](25,3) = 1; - restriction[4](26,2) = 1; - restriction[5](5,5) = 1; - restriction[5](12,4) = 1; - restriction[5](13,6) = 1; - restriction[5](17,1) = 1; - restriction[5](21,7) = 1; - restriction[5](22,0) = 1; - restriction[5](23,2) = 1; - restriction[5](26,3) = 1; - restriction[6](6,6) = 1; - restriction[6](13,5) = 1; - restriction[6](14,7) = 1; - restriction[6](18,2) = 1; - restriction[6](21,4) = 1; - restriction[6](23,1) = 1; - restriction[6](24,3) = 1; - restriction[6](26,0) = 1; - restriction[7](7,7) = 1; - restriction[7](14,6) = 1; - restriction[7](15,4) = 1; - restriction[7](19,3) = 1; - restriction[7](21,5) = 1; - restriction[7](24,2) = 1; - restriction[7](25,0) = 1; - restriction[7](26,1) = 1; - break; - case 3: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](14,15) = 1; - restriction[0](24,25) = 1; - restriction[0](32,35) = 1; - restriction[0](40,43) = 1; - restriction[0](52,55) = 1; - restriction[0](56,63) = 1; - restriction[1](1,1) = 1; - restriction[1](9,8) = 1; - restriction[1](10,11) = 1; - restriction[1](26,27) = 1; - restriction[1](33,34) = 1; - restriction[1](41,42) = 1; - restriction[1](44,47) = 1; - restriction[1](57,62) = 1; - restriction[2](2,2) = 1; - restriction[2](11,10) = 1; - restriction[2](13,12) = 1; - restriction[2](28,29) = 1; - restriction[2](35,32) = 1; - restriction[2](46,45) = 1; - restriction[2](49,50) = 1; - restriction[2](61,58) = 1; - restriction[3](3,3) = 1; - restriction[3](12,13) = 1; - restriction[3](15,14) = 1; - restriction[3](30,31) = 1; - restriction[3](34,33) = 1; - restriction[3](48,51) = 1; - restriction[3](54,53) = 1; - restriction[3](60,59) = 1; - restriction[4](4,4) = 1; - restriction[4](16,17) = 1; - restriction[4](22,23) = 1; - restriction[4](25,24) = 1; - restriction[4](36,39) = 1; - restriction[4](42,41) = 1; - restriction[4](53,54) = 1; - restriction[4](58,61) = 1; - restriction[5](5,5) = 1; - restriction[5](17,16) = 1; - restriction[5](18,19) = 1; - restriction[5](27,26) = 1; - restriction[5](37,38) = 1; - restriction[5](43,40) = 1; - restriction[5](45,46) = 1; - restriction[5](59,60) = 1; - restriction[6](6,6) = 1; - restriction[6](19,18) = 1; - restriction[6](21,20) = 1; - restriction[6](29,28) = 1; - restriction[6](39,36) = 1; - restriction[6](47,44) = 1; - restriction[6](51,48) = 1; - restriction[6](63,56) = 1; - restriction[7](7,7) = 1; - restriction[7](20,21) = 1; - restriction[7](23,22) = 1; - restriction[7](31,30) = 1; - restriction[7](38,37) = 1; - restriction[7](50,49) = 1; - restriction[7](55,52) = 1; - restriction[7](62,57) = 1; - break; - case 4: - restriction[0](0,0) = 1; - restriction[0](8,9) = 1; - restriction[0](9,1) = 1; - restriction[0](17,18) = 1; - restriction[0](18,3) = 1; - restriction[0](32,33) = 1; - restriction[0](33,4) = 1; - restriction[0](44,48) = 1; - restriction[0](45,12) = 1; - restriction[0](47,15) = 1; - restriction[0](48,2) = 1; - restriction[0](62,66) = 1; - restriction[0](63,36) = 1; - restriction[0](65,21) = 1; - restriction[0](66,5) = 1; - restriction[0](89,93) = 1; - restriction[0](90,30) = 1; - restriction[0](92,42) = 1; - restriction[0](93,7) = 1; - restriction[0](98,111) = 1; - restriction[0](99,75) = 1; - restriction[0](101,57) = 1; - restriction[0](102,24) = 1; - restriction[0](107,84) = 1; - restriction[0](108,39) = 1; - restriction[0](110,27) = 1; - restriction[0](111,6) = 1; - restriction[1](1,1) = 1; - restriction[1](9,0) = 1; - restriction[1](10,9) = 1; - restriction[1](11,12) = 1; - restriction[1](12,2) = 1; - restriction[1](35,36) = 1; - restriction[1](36,5) = 1; - restriction[1](45,18) = 1; - restriction[1](46,48) = 1; - restriction[1](48,3) = 1; - restriction[1](49,15) = 1; - restriction[1](63,33) = 1; - restriction[1](64,66) = 1; - restriction[1](66,4) = 1; - restriction[1](67,21) = 1; - restriction[1](71,75) = 1; - restriction[1](72,24) = 1; - restriction[1](74,39) = 1; - restriction[1](75,6) = 1; - restriction[1](99,93) = 1; - restriction[1](100,111) = 1; - restriction[1](102,30) = 1; - restriction[1](103,57) = 1; - restriction[1](108,42) = 1; - restriction[1](109,84) = 1; - restriction[1](111,7) = 1; - restriction[1](112,27) = 1; - restriction[2](2,2) = 1; - restriction[2](12,1) = 1; - restriction[2](13,12) = 1; - restriction[2](15,3) = 1; - restriction[2](16,15) = 1; - restriction[2](38,39) = 1; - restriction[2](39,6) = 1; - restriction[2](48,0) = 1; - restriction[2](49,9) = 1; - restriction[2](51,18) = 1; - restriction[2](52,48) = 1; - restriction[2](74,36) = 1; - restriction[2](75,5) = 1; - restriction[2](77,75) = 1; - restriction[2](78,24) = 1; - restriction[2](81,42) = 1; - restriction[2](82,84) = 1; - restriction[2](84,7) = 1; - restriction[2](85,27) = 1; - restriction[2](108,33) = 1; - restriction[2](109,66) = 1; - restriction[2](111,4) = 1; - restriction[2](112,21) = 1; - restriction[2](117,93) = 1; - restriction[2](118,111) = 1; - restriction[2](120,30) = 1; - restriction[2](121,57) = 1; - restriction[3](3,3) = 1; - restriction[3](14,15) = 1; - restriction[3](15,2) = 1; - restriction[3](18,0) = 1; - restriction[3](19,18) = 1; - restriction[3](41,42) = 1; - restriction[3](42,7) = 1; - restriction[3](47,9) = 1; - restriction[3](48,1) = 1; - restriction[3](50,48) = 1; - restriction[3](51,12) = 1; - restriction[3](80,84) = 1; - restriction[3](81,39) = 1; - restriction[3](83,27) = 1; - restriction[3](84,6) = 1; - restriction[3](92,33) = 1; - restriction[3](93,4) = 1; - restriction[3](95,93) = 1; - restriction[3](96,30) = 1; - restriction[3](107,66) = 1; - restriction[3](108,36) = 1; - restriction[3](110,21) = 1; - restriction[3](111,5) = 1; - restriction[3](116,111) = 1; - restriction[3](117,75) = 1; - restriction[3](119,57) = 1; - restriction[3](120,24) = 1; - restriction[4](4,4) = 1; - restriction[4](20,21) = 1; - restriction[4](21,5) = 1; - restriction[4](29,30) = 1; - restriction[4](30,7) = 1; - restriction[4](33,0) = 1; - restriction[4](34,33) = 1; - restriction[4](53,57) = 1; - restriction[4](54,24) = 1; - restriction[4](56,27) = 1; - restriction[4](57,6) = 1; - restriction[4](65,9) = 1; - restriction[4](66,1) = 1; - restriction[4](68,66) = 1; - restriction[4](69,36) = 1; - restriction[4](90,18) = 1; - restriction[4](91,93) = 1; - restriction[4](93,3) = 1; - restriction[4](94,42) = 1; - restriction[4](101,48) = 1; - restriction[4](102,12) = 1; - restriction[4](104,111) = 1; - restriction[4](105,75) = 1; - restriction[4](110,15) = 1; - restriction[4](111,2) = 1; - restriction[4](113,84) = 1; - restriction[4](114,39) = 1; - restriction[5](5,5) = 1; - restriction[5](21,4) = 1; - restriction[5](22,21) = 1; - restriction[5](23,24) = 1; - restriction[5](24,6) = 1; - restriction[5](36,1) = 1; - restriction[5](37,36) = 1; - restriction[5](54,30) = 1; - restriction[5](55,57) = 1; - restriction[5](57,7) = 1; - restriction[5](58,27) = 1; - restriction[5](66,0) = 1; - restriction[5](67,9) = 1; - restriction[5](69,33) = 1; - restriction[5](70,66) = 1; - restriction[5](72,12) = 1; - restriction[5](73,75) = 1; - restriction[5](75,2) = 1; - restriction[5](76,39) = 1; - restriction[5](102,18) = 1; - restriction[5](103,48) = 1; - restriction[5](105,93) = 1; - restriction[5](106,111) = 1; - restriction[5](111,3) = 1; - restriction[5](112,15) = 1; - restriction[5](114,42) = 1; - restriction[5](115,84) = 1; - restriction[6](6,6) = 1; - restriction[6](24,5) = 1; - restriction[6](25,24) = 1; - restriction[6](27,7) = 1; - restriction[6](28,27) = 1; - restriction[6](39,2) = 1; - restriction[6](40,39) = 1; - restriction[6](57,4) = 1; - restriction[6](58,21) = 1; - restriction[6](60,30) = 1; - restriction[6](61,57) = 1; - restriction[6](75,1) = 1; - restriction[6](76,36) = 1; - restriction[6](78,12) = 1; - restriction[6](79,75) = 1; - restriction[6](84,3) = 1; - restriction[6](85,15) = 1; - restriction[6](87,42) = 1; - restriction[6](88,84) = 1; - restriction[6](111,0) = 1; - restriction[6](112,9) = 1; - restriction[6](114,33) = 1; - restriction[6](115,66) = 1; - restriction[6](120,18) = 1; - restriction[6](121,48) = 1; - restriction[6](123,93) = 1; - restriction[6](124,111) = 1; - restriction[7](7,7) = 1; - restriction[7](26,27) = 1; - restriction[7](27,6) = 1; - restriction[7](30,4) = 1; - restriction[7](31,30) = 1; - restriction[7](42,3) = 1; - restriction[7](43,42) = 1; - restriction[7](56,21) = 1; - restriction[7](57,5) = 1; - restriction[7](59,57) = 1; - restriction[7](60,24) = 1; - restriction[7](83,15) = 1; - restriction[7](84,2) = 1; - restriction[7](86,84) = 1; - restriction[7](87,39) = 1; - restriction[7](93,0) = 1; - restriction[7](94,33) = 1; - restriction[7](96,18) = 1; - restriction[7](97,93) = 1; - restriction[7](110,9) = 1; - restriction[7](111,1) = 1; - restriction[7](113,66) = 1; - restriction[7](114,36) = 1; - restriction[7](119,48) = 1; - restriction[7](120,12) = 1; - restriction[7](122,111) = 1; - restriction[7](123,75) = 1; - break; - default: - for (unsigned int i=0; i::children_per_cell;++i) - restriction[i].reinit(0,0); - } - break; - default: - Assert (false,ExcNotImplemented()); + case 1: + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](2,1) = 1; + restriction[1](1,1) = 1; + restriction[1](1,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](2,3) = 1; + restriction[1](1,1) = 1; + restriction[1](3,2) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](2,3) = 1; + restriction[0](3,1) = 1; + restriction[1](1,1) = 1; + restriction[1](3,0) = 1; + restriction[1](4,3) = 1; + break; + default: + for (unsigned int i=0; i::children_per_cell;++i) + restriction[i].reinit(0,0); + } + break; + case 2: + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + restriction[2](2,2) = 1; + restriction[3](3,3) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](4,1) = 1; + restriction[0](7,3) = 1; + restriction[0](8,2) = 1; + restriction[1](1,1) = 1; + restriction[1](4,0) = 1; + restriction[1](5,2) = 1; + restriction[1](8,3) = 1; + restriction[2](2,2) = 1; + restriction[2](5,1) = 1; + restriction[2](6,3) = 1; + restriction[2](8,0) = 1; + restriction[3](3,3) = 1; + restriction[3](6,2) = 1; + restriction[3](7,0) = 1; + restriction[3](8,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](4,5) = 1; + restriction[0](10,11) = 1; + restriction[0](12,15) = 1; + restriction[1](1,1) = 1; + restriction[1](5,4) = 1; + restriction[1](6,7) = 1; + restriction[1](13,14) = 1; + restriction[2](2,2) = 1; + restriction[2](7,6) = 1; + restriction[2](9,8) = 1; + restriction[2](15,12) = 1; + restriction[3](3,3) = 1; + restriction[3](8,9) = 1; + restriction[3](11,10) = 1; + restriction[3](14,13) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](4,5) = 1; + restriction[0](5,1) = 1; + restriction[0](13,14) = 1; + restriction[0](14,3) = 1; + restriction[0](16,20) = 1; + restriction[0](17,8) = 1; + restriction[0](19,11) = 1; + restriction[0](20,2) = 1; + restriction[1](1,1) = 1; + restriction[1](5,0) = 1; + restriction[1](6,5) = 1; + restriction[1](7,8) = 1; + restriction[1](8,2) = 1; + restriction[1](17,14) = 1; + restriction[1](18,20) = 1; + restriction[1](20,3) = 1; + restriction[1](21,11) = 1; + restriction[2](2,2) = 1; + restriction[2](8,1) = 1; + restriction[2](9,8) = 1; + restriction[2](11,3) = 1; + restriction[2](12,11) = 1; + restriction[2](20,0) = 1; + restriction[2](21,5) = 1; + restriction[2](23,14) = 1; + restriction[2](24,20) = 1; + restriction[3](3,3) = 1; + restriction[3](10,11) = 1; + restriction[3](11,2) = 1; + restriction[3](14,0) = 1; + restriction[3](15,14) = 1; + restriction[3](19,5) = 1; + restriction[3](20,1) = 1; + restriction[3](22,20) = 1; + restriction[3](23,8) = 1; + break; + default: + for (unsigned int i=0; i::children_per_cell;++i) + restriction[i].reinit(0,0); + } + break; + case 3: + switch (degree) + { + case 1: + restriction[0](0,0) = 1; + restriction[1](1,1) = 1; + restriction[2](2,2) = 1; + restriction[3](3,3) = 1; + restriction[4](4,4) = 1; + restriction[5](5,5) = 1; + restriction[6](6,6) = 1; + restriction[7](7,7) = 1; + break; + case 2: + restriction[0](0,0) = 1; + restriction[0](8,1) = 1; + restriction[0](11,3) = 1; + restriction[0](16,4) = 1; + restriction[0](20,2) = 1; + restriction[0](22,5) = 1; + restriction[0](25,7) = 1; + restriction[0](26,6) = 1; + restriction[1](1,1) = 1; + restriction[1](8,0) = 1; + restriction[1](9,2) = 1; + restriction[1](17,5) = 1; + restriction[1](20,3) = 1; + restriction[1](22,4) = 1; + restriction[1](23,6) = 1; + restriction[1](26,7) = 1; + restriction[2](2,2) = 1; + restriction[2](9,1) = 1; + restriction[2](10,3) = 1; + restriction[2](18,6) = 1; + restriction[2](20,0) = 1; + restriction[2](23,5) = 1; + restriction[2](24,7) = 1; + restriction[2](26,4) = 1; + restriction[3](3,3) = 1; + restriction[3](10,2) = 1; + restriction[3](11,0) = 1; + restriction[3](19,7) = 1; + restriction[3](20,1) = 1; + restriction[3](24,6) = 1; + restriction[3](25,4) = 1; + restriction[3](26,5) = 1; + restriction[4](4,4) = 1; + restriction[4](12,5) = 1; + restriction[4](15,7) = 1; + restriction[4](16,0) = 1; + restriction[4](21,6) = 1; + restriction[4](22,1) = 1; + restriction[4](25,3) = 1; + restriction[4](26,2) = 1; + restriction[5](5,5) = 1; + restriction[5](12,4) = 1; + restriction[5](13,6) = 1; + restriction[5](17,1) = 1; + restriction[5](21,7) = 1; + restriction[5](22,0) = 1; + restriction[5](23,2) = 1; + restriction[5](26,3) = 1; + restriction[6](6,6) = 1; + restriction[6](13,5) = 1; + restriction[6](14,7) = 1; + restriction[6](18,2) = 1; + restriction[6](21,4) = 1; + restriction[6](23,1) = 1; + restriction[6](24,3) = 1; + restriction[6](26,0) = 1; + restriction[7](7,7) = 1; + restriction[7](14,6) = 1; + restriction[7](15,4) = 1; + restriction[7](19,3) = 1; + restriction[7](21,5) = 1; + restriction[7](24,2) = 1; + restriction[7](25,0) = 1; + restriction[7](26,1) = 1; + break; + case 3: + restriction[0](0,0) = 1; + restriction[0](8,9) = 1; + restriction[0](14,15) = 1; + restriction[0](24,25) = 1; + restriction[0](32,35) = 1; + restriction[0](40,43) = 1; + restriction[0](52,55) = 1; + restriction[0](56,63) = 1; + restriction[1](1,1) = 1; + restriction[1](9,8) = 1; + restriction[1](10,11) = 1; + restriction[1](26,27) = 1; + restriction[1](33,34) = 1; + restriction[1](41,42) = 1; + restriction[1](44,47) = 1; + restriction[1](57,62) = 1; + restriction[2](2,2) = 1; + restriction[2](11,10) = 1; + restriction[2](13,12) = 1; + restriction[2](28,29) = 1; + restriction[2](35,32) = 1; + restriction[2](46,45) = 1; + restriction[2](49,50) = 1; + restriction[2](61,58) = 1; + restriction[3](3,3) = 1; + restriction[3](12,13) = 1; + restriction[3](15,14) = 1; + restriction[3](30,31) = 1; + restriction[3](34,33) = 1; + restriction[3](48,51) = 1; + restriction[3](54,53) = 1; + restriction[3](60,59) = 1; + restriction[4](4,4) = 1; + restriction[4](16,17) = 1; + restriction[4](22,23) = 1; + restriction[4](25,24) = 1; + restriction[4](36,39) = 1; + restriction[4](42,41) = 1; + restriction[4](53,54) = 1; + restriction[4](58,61) = 1; + restriction[5](5,5) = 1; + restriction[5](17,16) = 1; + restriction[5](18,19) = 1; + restriction[5](27,26) = 1; + restriction[5](37,38) = 1; + restriction[5](43,40) = 1; + restriction[5](45,46) = 1; + restriction[5](59,60) = 1; + restriction[6](6,6) = 1; + restriction[6](19,18) = 1; + restriction[6](21,20) = 1; + restriction[6](29,28) = 1; + restriction[6](39,36) = 1; + restriction[6](47,44) = 1; + restriction[6](51,48) = 1; + restriction[6](63,56) = 1; + restriction[7](7,7) = 1; + restriction[7](20,21) = 1; + restriction[7](23,22) = 1; + restriction[7](31,30) = 1; + restriction[7](38,37) = 1; + restriction[7](50,49) = 1; + restriction[7](55,52) = 1; + restriction[7](62,57) = 1; + break; + case 4: + restriction[0](0,0) = 1; + restriction[0](8,9) = 1; + restriction[0](9,1) = 1; + restriction[0](17,18) = 1; + restriction[0](18,3) = 1; + restriction[0](32,33) = 1; + restriction[0](33,4) = 1; + restriction[0](44,48) = 1; + restriction[0](45,12) = 1; + restriction[0](47,15) = 1; + restriction[0](48,2) = 1; + restriction[0](62,66) = 1; + restriction[0](63,36) = 1; + restriction[0](65,21) = 1; + restriction[0](66,5) = 1; + restriction[0](89,93) = 1; + restriction[0](90,30) = 1; + restriction[0](92,42) = 1; + restriction[0](93,7) = 1; + restriction[0](98,111) = 1; + restriction[0](99,75) = 1; + restriction[0](101,57) = 1; + restriction[0](102,24) = 1; + restriction[0](107,84) = 1; + restriction[0](108,39) = 1; + restriction[0](110,27) = 1; + restriction[0](111,6) = 1; + restriction[1](1,1) = 1; + restriction[1](9,0) = 1; + restriction[1](10,9) = 1; + restriction[1](11,12) = 1; + restriction[1](12,2) = 1; + restriction[1](35,36) = 1; + restriction[1](36,5) = 1; + restriction[1](45,18) = 1; + restriction[1](46,48) = 1; + restriction[1](48,3) = 1; + restriction[1](49,15) = 1; + restriction[1](63,33) = 1; + restriction[1](64,66) = 1; + restriction[1](66,4) = 1; + restriction[1](67,21) = 1; + restriction[1](71,75) = 1; + restriction[1](72,24) = 1; + restriction[1](74,39) = 1; + restriction[1](75,6) = 1; + restriction[1](99,93) = 1; + restriction[1](100,111) = 1; + restriction[1](102,30) = 1; + restriction[1](103,57) = 1; + restriction[1](108,42) = 1; + restriction[1](109,84) = 1; + restriction[1](111,7) = 1; + restriction[1](112,27) = 1; + restriction[2](2,2) = 1; + restriction[2](12,1) = 1; + restriction[2](13,12) = 1; + restriction[2](15,3) = 1; + restriction[2](16,15) = 1; + restriction[2](38,39) = 1; + restriction[2](39,6) = 1; + restriction[2](48,0) = 1; + restriction[2](49,9) = 1; + restriction[2](51,18) = 1; + restriction[2](52,48) = 1; + restriction[2](74,36) = 1; + restriction[2](75,5) = 1; + restriction[2](77,75) = 1; + restriction[2](78,24) = 1; + restriction[2](81,42) = 1; + restriction[2](82,84) = 1; + restriction[2](84,7) = 1; + restriction[2](85,27) = 1; + restriction[2](108,33) = 1; + restriction[2](109,66) = 1; + restriction[2](111,4) = 1; + restriction[2](112,21) = 1; + restriction[2](117,93) = 1; + restriction[2](118,111) = 1; + restriction[2](120,30) = 1; + restriction[2](121,57) = 1; + restriction[3](3,3) = 1; + restriction[3](14,15) = 1; + restriction[3](15,2) = 1; + restriction[3](18,0) = 1; + restriction[3](19,18) = 1; + restriction[3](41,42) = 1; + restriction[3](42,7) = 1; + restriction[3](47,9) = 1; + restriction[3](48,1) = 1; + restriction[3](50,48) = 1; + restriction[3](51,12) = 1; + restriction[3](80,84) = 1; + restriction[3](81,39) = 1; + restriction[3](83,27) = 1; + restriction[3](84,6) = 1; + restriction[3](92,33) = 1; + restriction[3](93,4) = 1; + restriction[3](95,93) = 1; + restriction[3](96,30) = 1; + restriction[3](107,66) = 1; + restriction[3](108,36) = 1; + restriction[3](110,21) = 1; + restriction[3](111,5) = 1; + restriction[3](116,111) = 1; + restriction[3](117,75) = 1; + restriction[3](119,57) = 1; + restriction[3](120,24) = 1; + restriction[4](4,4) = 1; + restriction[4](20,21) = 1; + restriction[4](21,5) = 1; + restriction[4](29,30) = 1; + restriction[4](30,7) = 1; + restriction[4](33,0) = 1; + restriction[4](34,33) = 1; + restriction[4](53,57) = 1; + restriction[4](54,24) = 1; + restriction[4](56,27) = 1; + restriction[4](57,6) = 1; + restriction[4](65,9) = 1; + restriction[4](66,1) = 1; + restriction[4](68,66) = 1; + restriction[4](69,36) = 1; + restriction[4](90,18) = 1; + restriction[4](91,93) = 1; + restriction[4](93,3) = 1; + restriction[4](94,42) = 1; + restriction[4](101,48) = 1; + restriction[4](102,12) = 1; + restriction[4](104,111) = 1; + restriction[4](105,75) = 1; + restriction[4](110,15) = 1; + restriction[4](111,2) = 1; + restriction[4](113,84) = 1; + restriction[4](114,39) = 1; + restriction[5](5,5) = 1; + restriction[5](21,4) = 1; + restriction[5](22,21) = 1; + restriction[5](23,24) = 1; + restriction[5](24,6) = 1; + restriction[5](36,1) = 1; + restriction[5](37,36) = 1; + restriction[5](54,30) = 1; + restriction[5](55,57) = 1; + restriction[5](57,7) = 1; + restriction[5](58,27) = 1; + restriction[5](66,0) = 1; + restriction[5](67,9) = 1; + restriction[5](69,33) = 1; + restriction[5](70,66) = 1; + restriction[5](72,12) = 1; + restriction[5](73,75) = 1; + restriction[5](75,2) = 1; + restriction[5](76,39) = 1; + restriction[5](102,18) = 1; + restriction[5](103,48) = 1; + restriction[5](105,93) = 1; + restriction[5](106,111) = 1; + restriction[5](111,3) = 1; + restriction[5](112,15) = 1; + restriction[5](114,42) = 1; + restriction[5](115,84) = 1; + restriction[6](6,6) = 1; + restriction[6](24,5) = 1; + restriction[6](25,24) = 1; + restriction[6](27,7) = 1; + restriction[6](28,27) = 1; + restriction[6](39,2) = 1; + restriction[6](40,39) = 1; + restriction[6](57,4) = 1; + restriction[6](58,21) = 1; + restriction[6](60,30) = 1; + restriction[6](61,57) = 1; + restriction[6](75,1) = 1; + restriction[6](76,36) = 1; + restriction[6](78,12) = 1; + restriction[6](79,75) = 1; + restriction[6](84,3) = 1; + restriction[6](85,15) = 1; + restriction[6](87,42) = 1; + restriction[6](88,84) = 1; + restriction[6](111,0) = 1; + restriction[6](112,9) = 1; + restriction[6](114,33) = 1; + restriction[6](115,66) = 1; + restriction[6](120,18) = 1; + restriction[6](121,48) = 1; + restriction[6](123,93) = 1; + restriction[6](124,111) = 1; + restriction[7](7,7) = 1; + restriction[7](26,27) = 1; + restriction[7](27,6) = 1; + restriction[7](30,4) = 1; + restriction[7](31,30) = 1; + restriction[7](42,3) = 1; + restriction[7](43,42) = 1; + restriction[7](56,21) = 1; + restriction[7](57,5) = 1; + restriction[7](59,57) = 1; + restriction[7](60,24) = 1; + restriction[7](83,15) = 1; + restriction[7](84,2) = 1; + restriction[7](86,84) = 1; + restriction[7](87,39) = 1; + restriction[7](93,0) = 1; + restriction[7](94,33) = 1; + restriction[7](96,18) = 1; + restriction[7](97,93) = 1; + restriction[7](110,9) = 1; + restriction[7](111,1) = 1; + restriction[7](113,66) = 1; + restriction[7](114,36) = 1; + restriction[7](119,48) = 1; + restriction[7](120,12) = 1; + restriction[7](122,111) = 1; + restriction[7](123,75) = 1; + break; + default: + for (unsigned int i=0; i::children_per_cell;++i) + restriction[i].reinit(0,0); + } + break; + default: + Assert (false,ExcNotImplemented()); } } @@ -596,6 +597,7 @@ FE_Q::~FE_Q () } + template FiniteElement * FE_Q::clone() const @@ -604,19 +606,21 @@ FE_Q::clone() const } + template void FE_Q::get_unit_support_points (std::vector > &points) const { - compute_support_points (points, degree, renumber); + compute_support_points (degree, renumber, points); } + template void FE_Q::get_unit_face_support_points (std::vector > &points) const { - FE_Q::compute_support_points (points, degree, face_renumber); + FE_Q::compute_support_points (degree, face_renumber, points); } @@ -630,31 +634,34 @@ FE_Q<1>::get_unit_face_support_points (std::vector > &points) const } #endif + + + //---------------------------------------------------------------------- -// Auxilliary functions +// Auxiliary functions //---------------------------------------------------------------------- template void -FE_Q::compute_support_points (std::vector >& support_points, - unsigned int degree, - const std::vector& renumber) +FE_Q::compute_support_points (const unsigned int degree, + const std::vector &renumber, + std::vector > &support_points) { // number of points: (degree+1)^dim - unsigned int n= degree+1; - for (unsigned int i=1;i p; unsigned int k=0; - for (unsigned int iz=0;iz <= ((dim>2) ? degree : 0) ; ++iz) - for (unsigned int iy=0;iy <= ((dim>1) ? degree : 0) ; ++iy) - for (unsigned int ix=0;ix<=degree;++ix) + for (unsigned int iz=0; iz <= ((dim>2) ? degree : 0) ; ++iz) + for (unsigned int iy=0; iy <= ((dim>1) ? degree : 0) ; ++iy) + for (unsigned int ix=0; ix<=degree; ++ix) { p(0) = ix * step; if (dim>1) @@ -669,7 +676,7 @@ FE_Q::compute_support_points (std::vector >& support_points, template std::vector -FE_Q::get_dpo_vector(unsigned int deg) +FE_Q::get_dpo_vector(const unsigned int deg) { std::vector dpo(dim+1, 1); for (unsigned int i=1; i::get_dpo_vector(unsigned int deg) template void FE_Q::build_renumbering (const FiniteElementData &fe_data, - unsigned int degree, - std::vector &renumber) + const unsigned int degree, + std::vector &renumber) { const unsigned int n = degree+1; @@ -702,58 +709,59 @@ FE_Q::build_renumbering (const FiniteElementData &fe_data, // sophisticated. switch (dim) { - case 1: - if (i==1) - index += degree; - break; - case 2: - switch (i) - { - case 1: - index += degree; - break; - case 3: - index += n*degree; - break; - case 2: - index += n*degree+degree; - break; - } - break; - case 3: - switch (i) - { - case 1: - index += degree; - break; - case 4: - index += n*degree; - break; - case 5: - index += n*degree+degree; - break; - case 3: - index += n*n*degree; - break; - case 2: - index += n*n*degree + degree; - break; - case 7: - index += n*n*degree + n*degree; - break; - case 6: - index += n*n*degree + n*degree+degree; - break; - } - break; - - default: - Assert(false, ExcNotImplemented()); + case 1: + if (i==1) + index += degree; + break; + case 2: + switch (i) + { + case 1: + index += degree; + break; + case 3: + index += n*degree; + break; + case 2: + index += n*degree+degree; + break; + } + break; + case 3: + switch (i) + { + case 1: + index += degree; + break; + case 4: + index += n*degree; + break; + case 5: + index += n*degree+degree; + break; + case 3: + index += n*n*degree; + break; + case 2: + index += n*n*degree + degree; + break; + case 7: + index += n*n*degree + n*degree; + break; + case 6: + index += n*n*degree + n*degree+degree; + break; + } + break; + + default: + Assert(false, ExcNotImplemented()); } renumber[index] = i; } else + // degree == 0 renumber[0] = 0; // Lines and higher @@ -769,62 +777,62 @@ FE_Q::build_renumbering (const FiniteElementData &fe_data, // numbering scheme. switch (i+100*dim) { - // lines in x-direction - case 100: - case 200: case 202: - case 300: case 302: case 304: case 306: - incr = 1; - break; - // lines in y-direction - case 201: case 203: - case 308: case 309: case 310: case 311: - incr = n; - break; - // lines in z-direction - case 301: case 303: case 305: case 307: - incr = n*n; - break; - default: - Assert(false, ExcNotImplemented()); + // lines in x-direction + case 100: + case 200: case 202: + case 300: case 302: case 304: case 306: + incr = 1; + break; + // lines in y-direction + case 201: case 203: + case 308: case 309: case 310: case 311: + incr = n; + break; + // lines in z-direction + case 301: case 303: case 305: case 307: + incr = n*n; + break; + default: + Assert(false, ExcNotImplemented()); } switch (i+100*dim) { - // x=y=z=0 - case 100: - case 200: case 203: - case 300: case 303: case 308: - tensorstart = 0; - break; - // x=1 y=z=0 - case 201: - case 301: case 309: - tensorstart = degree; - break; - // y=1 x=z=0 - case 202: - case 304: case 307: - tensorstart = n*degree; - break; - // x=z=1 y=0 - case 310: - tensorstart = n*n*degree+degree; - break; - // z=1 x=y=0 - case 302: case 311: - tensorstart = n*n*degree; - break; - // x=y=1 z=0 - case 305: - tensorstart = n*degree+degree; - break; - // y=z=1 x=0 - case 306: - tensorstart = n*n*n-n; - break; - default: - Assert(false, ExcNotImplemented()); + // x=y=z=0 + case 100: + case 200: case 203: + case 300: case 303: case 308: + tensorstart = 0; + break; + // x=1 y=z=0 + case 201: + case 301: case 309: + tensorstart = degree; + break; + // y=1 x=z=0 + case 202: + case 304: case 307: + tensorstart = n*degree; + break; + // x=z=1 y=0 + case 310: + tensorstart = n*n*degree+degree; + break; + // z=1 x=y=0 + case 302: case 311: + tensorstart = n*n*degree; + break; + // x=y=1 z=0 + case 305: + tensorstart = n*degree+degree; + break; + // y=z=1 x=0 + case 306: + tensorstart = n*n*n-n; + break; + default: + Assert(false, ExcNotImplemented()); } - + for (unsigned int jx = 1; jx::build_renumbering (const FiniteElementData &fe_data, switch (i) { case 0: - tensorstart = 0; incx = 1; - if (dim==2) - incy = n; - else - incy = n*n; - break; + tensorstart = 0; incx = 1; + if (dim==2) + incy = n; + else + incy = n*n; + break; case 1: - tensorstart = n*degree; incx = 1; incy = n*n; - break; + tensorstart = n*degree; incx = 1; incy = n*n; + break; case 2: - tensorstart = 0; incx = 1; incy = n; - break; + tensorstart = 0; incx = 1; incy = n; + break; case 3: - tensorstart = degree; incx = n; incy = n*n; - break; + tensorstart = degree; incx = n; incy = n*n; + break; case 4: - tensorstart = n*n*degree; incx = 1; incy = n; - break; + tensorstart = n*n*degree; incx = 1; incy = n; + break; case 5: - tensorstart = 0; incx = n; incy = n*n; - break; - default: - Assert(false, ExcNotImplemented()); + tensorstart = 0; incx = n; incy = n*n; + break; + default: + Assert(false, ExcNotImplemented()); } for (unsigned int jy = 1; jy::build_renumbering (const FiniteElementData &fe_data, for (unsigned int jy = 1; jy::build_renumbering (const FiniteElementData &fe_data, } + template void FE_Q::build_face_renumbering (const FiniteElementData &fe_data, - unsigned int degree, - std::vector& numbering) + const unsigned int degree, + std::vector &numbering) { FE_Q::build_renumbering (fe_data, degree, numbering); } + #if (deal_II_dimension == 1) template <> void FE_Q<1>::build_face_renumbering (const FiniteElementData<0> &, - unsigned int, + const unsigned int, std::vector&) {} + #endif template UpdateFlags -FE_Q::update_once (UpdateFlags flags) const +FE_Q::update_once (const UpdateFlags flags) const { UpdateFlags out = update_default; @@ -926,6 +936,7 @@ FE_Q::update_once (UpdateFlags flags) const } + template UpdateFlags FE_Q::update_each (UpdateFlags flags) const @@ -941,13 +952,15 @@ FE_Q::update_each (UpdateFlags flags) const } + //---------------------------------------------------------------------- // Data field initialization //---------------------------------------------------------------------- + template Mapping::InternalDataBase* -FE_Q::get_data (const UpdateFlags update_flags, - const Mapping& mapping, +FE_Q::get_data (const UpdateFlags update_flags, + const Mapping &mapping, const Quadrature &quadrature) const { InternalData* data = new InternalData; @@ -960,29 +973,30 @@ FE_Q::get_data (const UpdateFlags update_flags, data->update_flags = data->update_once | data->update_each; const UpdateFlags flags(data->update_flags); + const unsigned int n_q_points = quadrature.n_quadrature_points; if (flags & update_values) { values.resize (dofs_per_cell); data->shape_values.resize(dofs_per_cell, - std::vector(quadrature.n_quadrature_points)); + std::vector(n_q_points)); } if (flags & update_gradients) { grads.resize (dofs_per_cell); data->shape_gradients.resize(dofs_per_cell, - std::vector >(quadrature.n_quadrature_points)); + std::vector >(n_q_points)); } if (flags & update_second_derivatives) data->initialize (this, mapping, quadrature); if (flags & (update_values | update_gradients)) - for (unsigned int i=0;icompute(quadrature.point(i), values, grads, grad_grads); - for (unsigned int k=0;kshape_values[renumber[k]][i] = values[k]; @@ -995,36 +1009,40 @@ FE_Q::get_data (const UpdateFlags update_flags, + //---------------------------------------------------------------------- // Fill data of FEValues //---------------------------------------------------------------------- + template void FE_Q::fill_fe_values (const Mapping &mapping, const DoFHandler::cell_iterator &cell, const Quadrature &quadrature, - Mapping::InternalDataBase &mapping_data, - Mapping::InternalDataBase &fedata, + Mapping::InternalDataBase &mapping_data, + Mapping::InternalDataBase &fedata, FEValuesData &data) const { + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible InternalData &fe_data = dynamic_cast (fedata); const UpdateFlags flags(fe_data.current_update_flags()); - for (unsigned int k=0;k::fill_fe_values (const Mapping &mapping, template void -FE_Q::fill_fe_face_values (const Mapping &mapping, +FE_Q::fill_fe_face_values (const Mapping &mapping, const DoFHandler::cell_iterator &cell, - const unsigned int face, - const Quadrature& quadrature, - Mapping::InternalDataBase& mapping_data, - Mapping::InternalDataBase& fedata, - FEValuesData& data) const + const unsigned int face, + const Quadrature &quadrature, + Mapping::InternalDataBase &mapping_data, + Mapping::InternalDataBase &fedata, + FEValuesData &data) const { + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible InternalData &fe_data = dynamic_cast (fedata); - unsigned int offset = face * quadrature.n_quadrature_points; + + // offset determines which data set + // to take (all data sets for all + // faces are stored contiguously) + const unsigned int offset = face * quadrature.n_quadrature_points; const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0;k::fill_fe_face_values (const Mapping &mapping, template void -FE_Q::fill_fe_subface_values (const Mapping &mapping, +FE_Q::fill_fe_subface_values (const Mapping &mapping, const DoFHandler::cell_iterator &cell, - const unsigned int face, - const unsigned int subface, - const Quadrature& quadrature, - Mapping::InternalDataBase& mapping_data, - Mapping::InternalDataBase& fedata, - FEValuesData& data) const + const unsigned int face, + const unsigned int subface, + const Quadrature &quadrature, + Mapping::InternalDataBase &mapping_data, + Mapping::InternalDataBase &fedata, + FEValuesData &data) const { + // convert data object to internal + // data for this class. fails with + // an exception if that is not + // possible InternalData &fe_data = dynamic_cast (fedata); - unsigned int offset = (face * GeometryInfo::subfaces_per_face - + subface) * quadrature.n_quadrature_points; + + // offset determines which data set + // to take (all data sets for all + // sub-faces are stored contiguously) + const unsigned int offset = (face * GeometryInfo::subfaces_per_face + subface) + * quadrature.n_quadrature_points; const UpdateFlags flags(fe_data.update_once | fe_data.update_each); - for (unsigned int k=0;k::Matrices::n_embedding_matrices = 4; + +const double * const +FE_Q<1>::Matrices::embedding[][GeometryInfo<1>::children_per_cell] = { - q1_into_q1_refined_1, - q2_into_q2_refined_1, - q3_into_q3_refined_1, - q4_into_q4_refined_1, - 0, - 0, - 0, - 0, - 0, - 0 + {q1_into_q1_refined_0, q1_into_q1_refined_1}, + {q2_into_q2_refined_0, q2_into_q2_refined_1}, + {q3_into_q3_refined_0, q3_into_q3_refined_1}, + {q4_into_q4_refined_0, q4_into_q4_refined_1}, }; -static double** embedding2; -static double** embedding3; -static double** embedding4; -static double** embedding5; -static double** embedding6; -static double** embedding7; -// No constrants in 1d +// No constraints in 1d +const unsigned int +FE_Q<1>::Matrices::n_constraint_matrices = 0; -static double ** constraint_matrices = 0; +const double * const +FE_Q<1>::Matrices::constraint_matrices[] = { 0 }; diff --git a/deal.II/deal.II/source/fe/mat_q.2 b/deal.II/deal.II/source/fe/mat_q.2 index 089a416c63..c52d9cc7bd 100644 --- a/deal.II/deal.II/source/fe/mat_q.2 +++ b/deal.II/deal.II/source/fe/mat_q.2 @@ -13,7 +13,7 @@ // Transfer matrices for finite elements -static double q1_into_q1_refined_0[] = +static const double q1_into_q1_refined_0[] = { 1., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., @@ -21,7 +21,7 @@ static double q1_into_q1_refined_0[] = 13.5/27., 0., 0., 13.5/27., }; -static double q1_into_q1_refined_1[] = +static const double q1_into_q1_refined_1[] = { 13.5/27., 13.5/27., 0., 0., 0., 1., 0., 0., @@ -29,7 +29,7 @@ static double q1_into_q1_refined_1[] = 6.75/27., 6.75/27., 6.75/27., 6.75/27., }; -static double q1_into_q1_refined_2[] = +static const double q1_into_q1_refined_2[] = { 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 13.5/27., 13.5/27., 0., @@ -37,7 +37,7 @@ static double q1_into_q1_refined_2[] = 0., 0., 13.5/27., 13.5/27., }; -static double q1_into_q1_refined_3[] = +static const double q1_into_q1_refined_3[] = { 13.5/27., 0., 0., 13.5/27., 6.75/27., 6.75/27., 6.75/27., 6.75/27., @@ -45,7 +45,7 @@ static double q1_into_q1_refined_3[] = 0., 0., 0., 1., }; -static double q2_into_q2_refined_0[] = +static const double q2_into_q2_refined_0[] = { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., @@ -58,7 +58,7 @@ static double q2_into_q2_refined_0[] = 3.796875/27., -1.265625/27., 0.421875/27., -1.265625/27., 7.59375/27., -2.53125/27., -2.53125/27., 7.59375/27., 15.1875/27., }; -static double q2_into_q2_refined_1[] = +static const double q2_into_q2_refined_1[] = { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., @@ -71,7 +71,7 @@ static double q2_into_q2_refined_1[] = -1.265625/27., 3.796875/27., -1.265625/27., 0.421875/27., 7.59375/27., 7.59375/27., -2.53125/27., -2.53125/27., 15.1875/27., }; -static double q2_into_q2_refined_2[] = +static const double q2_into_q2_refined_2[] = { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 1., 0., 0., 0., @@ -84,7 +84,7 @@ static double q2_into_q2_refined_2[] = 0.421875/27., -1.265625/27., 3.796875/27., -1.265625/27., -2.53125/27., 7.59375/27., 7.59375/27., -2.53125/27., 15.1875/27., }; -static double q2_into_q2_refined_3[] = +static const double q2_into_q2_refined_3[] = { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., @@ -97,7 +97,7 @@ static double q2_into_q2_refined_3[] = -1.265625/27., 0.421875/27., -1.265625/27., 3.796875/27., -2.53125/27., -2.53125/27., 7.59375/27., 7.59375/27., 15.1875/27., }; -static double q3_into_q3_refined_0[] = +static const double q3_into_q3_refined_0[] = { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -117,7 +117,7 @@ static double q3_into_q3_refined_0[] = 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., }; -static double q3_into_q3_refined_1[] = +static const double q3_into_q3_refined_1[] = { -1.6875/27., -1.6875/27., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -137,7 +137,7 @@ static double q3_into_q3_refined_1[] = 0., 0., 0., 0., 0., 0., 8.4375/27., 0., 0., 0., 1.6875/27., 0., -8.4375/27., 25.3125/27., 0., 0., }; -static double q3_into_q3_refined_2[] = +static const double q3_into_q3_refined_2[] = { 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 0., -1.6875/27., -1.6875/27., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0., 0., 0., 0., @@ -157,7 +157,7 @@ static double q3_into_q3_refined_2[] = 0.10546875/27., 0.52734375/27., 2.63671875/27., 0.52734375/27., -0.52734375/27., 1.58203125/27., -2.63671875/27., 7.91015625/27., -2.63671875/27., 7.91015625/27., -0.52734375/27., 1.58203125/27., 2.63671875/27., -7.91015625/27., -7.91015625/27., 23.73046875/27., }; -static double q3_into_q3_refined_3[] = +static const double q3_into_q3_refined_3[] = { -1.6875/27., 0., 0., -1.6875/27., 0., 0., 0., 0., 0., 0., 15.1875/27., 15.1875/27., 0., 0., 0., 0., 0.10546875/27., 0.10546875/27., 0.10546875/27., 0.10546875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., 8.54296875/27., @@ -178,82 +178,34 @@ static double q3_into_q3_refined_3[] = }; -static double * embedding0[] = -{ - q1_into_q1_refined_0, - q2_into_q2_refined_0, - q3_into_q3_refined_0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; +// embedding matrices -static double * embedding1[] = -{ - q1_into_q1_refined_1, - q2_into_q2_refined_1, - q3_into_q3_refined_1, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; +const unsigned int FE_Q<2>::Matrices::n_embedding_matrices = 3; -static double * embedding2[] = +const double * const +FE_Q<2>::Matrices::embedding[][GeometryInfo<2>::children_per_cell] = { - q1_into_q1_refined_2, - q2_into_q2_refined_2, - q3_into_q3_refined_2, - 0, - 0, - 0, - 0, - 0, - 0, - 0 + { q1_into_q1_refined_0, q1_into_q1_refined_1, q1_into_q1_refined_2, q1_into_q1_refined_3 }, + { q2_into_q2_refined_0, q2_into_q2_refined_1, q2_into_q2_refined_2, q2_into_q2_refined_3 }, + { q3_into_q3_refined_0, q3_into_q3_refined_1, q3_into_q3_refined_2, q3_into_q3_refined_3 } }; -static double * embedding3[] = -{ - q1_into_q1_refined_3, - q2_into_q2_refined_3, - q3_into_q3_refined_3, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; -static double** embedding4; -static double** embedding5; -static double** embedding6; -static double** embedding7; // Constraint matrices taken from Wolfgangs old version - -static double constraint_q1[] = +static const double constraint_q1[] = { .5, .5 }; -static double constraint_q2[] = +static const double constraint_q2[] = { 0., 0., 1., .375, -.125, .75, -.125, .375, .75 }; -static double constraint_q3[] = +static const double constraint_q3[] = { -.0625, -.0625, .5625, .5625, .3125, .0625, .9375, -.3125, @@ -262,7 +214,7 @@ static double constraint_q3[] = .0625, .3125, -.3125, 0.9375 }; -static double constraint_q4[] = +static const double constraint_q4[] = { 0., 0., 0., 1., 0., 0.2734375, -0.0390625, 1.09375, -0.546875, 0.21875, @@ -273,20 +225,15 @@ static double constraint_q4[] = -0.0390625, 0.2734375, 0.21875, -0.546875, 1.09375 }; -static double* constraint_matrices[] = +const double * const +FE_Q<2>::Matrices::constraint_matrices[] = { constraint_q1, constraint_q2, constraint_q3, constraint_q4, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 }; +const unsigned int +FE_Q<2>::Matrices::n_constraint_matrices = 4; + diff --git a/deal.II/deal.II/source/fe/mat_q.3 b/deal.II/deal.II/source/fe/mat_q.3 index 38791acf43..6159077b43 100644 --- a/deal.II/deal.II/source/fe/mat_q.3 +++ b/deal.II/deal.II/source/fe/mat_q.3 @@ -14,7 +14,7 @@ // Transfer matrices for finite elements -static double q1_into_q1_refined_0[] = +static const double q1_into_q1_refined_0[] = { 1., 0., 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., @@ -26,7 +26,7 @@ static double q1_into_q1_refined_0[] = 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., }; -static double q1_into_q1_refined_1[] = +static const double q1_into_q1_refined_1[] = { 13.5/27., 13.5/27., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., @@ -38,7 +38,7 @@ static double q1_into_q1_refined_1[] = 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., }; -static double q1_into_q1_refined_2[] = +static const double q1_into_q1_refined_2[] = { 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., 0., 0., 0., 0., 0., @@ -50,7 +50,7 @@ static double q1_into_q1_refined_2[] = 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., }; -static double q1_into_q1_refined_3[] = +static const double q1_into_q1_refined_3[] = { 13.5/27., 0., 0., 13.5/27., 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., 0., 0., 0., 0., @@ -62,7 +62,7 @@ static double q1_into_q1_refined_3[] = 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., }; -static double q1_into_q1_refined_4[] = +static const double q1_into_q1_refined_4[] = { 13.5/27., 0., 0., 0., 13.5/27., 0., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., @@ -74,7 +74,7 @@ static double q1_into_q1_refined_4[] = 0., 0., 0., 0., 13.5/27., 0., 0., 13.5/27., }; -static double q1_into_q1_refined_5[] = +static const double q1_into_q1_refined_5[] = { 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 0., 13.5/27., 0., 0., 0., 13.5/27., 0., 0., @@ -86,7 +86,7 @@ static double q1_into_q1_refined_5[] = 0., 0., 0., 0., 6.75/27., 6.75/27., 6.75/27., 6.75/27., }; -static double q1_into_q1_refined_6[] = +static const double q1_into_q1_refined_6[] = { 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., @@ -98,7 +98,7 @@ static double q1_into_q1_refined_6[] = 0., 0., 0., 0., 0., 0., 13.5/27., 13.5/27., }; -static double q1_into_q1_refined_7[] = +static const double q1_into_q1_refined_7[] = { 6.75/27., 0., 0., 6.75/27., 6.75/27., 0., 0., 6.75/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., 3.375/27., @@ -110,7 +110,7 @@ static double q1_into_q1_refined_7[] = 0., 0., 0., 0., 0., 0., 0., 1., }; -static double q2_into_q2_refined_0[] = +static const double q2_into_q2_refined_0[] = { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -141,7 +141,7 @@ static double q2_into_q2_refined_0[] = 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27., }; -static double q2_into_q2_refined_1[] = +static const double q2_into_q2_refined_1[] = { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -172,7 +172,7 @@ static double q2_into_q2_refined_1[] = -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27., }; -static double q2_into_q2_refined_2[] = +static const double q2_into_q2_refined_2[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -203,7 +203,7 @@ static double q2_into_q2_refined_2[] = 0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27., }; -static double q2_into_q2_refined_3[] = +static const double q2_into_q2_refined_3[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., @@ -234,7 +234,7 @@ static double q2_into_q2_refined_3[] = -0.474609375/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 11.390625/27., }; -static double q2_into_q2_refined_4[] = +static const double q2_into_q2_refined_4[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., @@ -265,7 +265,7 @@ static double q2_into_q2_refined_4[] = -0.474609375/27., 0.158203125/27., -0.052734375/27., 0.158203125/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.474609375/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 5.6953125/27., 11.390625/27., }; -static double q2_into_q2_refined_5[] = +static const double q2_into_q2_refined_5[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., @@ -296,7 +296,7 @@ static double q2_into_q2_refined_5[] = 0.158203125/27., -0.474609375/27., 0.158203125/27., -0.052734375/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.158203125/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., 0.31640625/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., -0.94921875/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., -1.8984375/27., 11.390625/27., }; -static double q2_into_q2_refined_6[] = +static const double q2_into_q2_refined_6[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., @@ -327,7 +327,7 @@ static double q2_into_q2_refined_6[] = -0.052734375/27., 0.158203125/27., -0.474609375/27., 0.158203125/27., 0.158203125/27., -0.474609375/27., 1.423828125/27., -0.474609375/27., 0.31640625/27., -0.94921875/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., 2.84765625/27., -0.94921875/27., 0.31640625/27., -0.94921875/27., 2.84765625/27., -0.94921875/27., -1.8984375/27., 5.6953125/27., -1.8984375/27., 5.6953125/27., 5.6953125/27., -1.8984375/27., 11.390625/27., }; -static double q2_into_q2_refined_7[] = +static const double q2_into_q2_refined_7[] = { 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., @@ -360,146 +360,27 @@ static double q2_into_q2_refined_7[] = +// embedding matrices -static double * embedding0[] = -{ - q1_into_q1_refined_0, - q2_into_q2_refined_0, -// q3_into_q3_refined_0, -// q4_into_q4_refined_0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; - -static double * embedding1[] = -{ - q1_into_q1_refined_1, - q2_into_q2_refined_1, -// q3_into_q3_refined_1, -// q4_into_q4_refined_1, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; +const unsigned int +FE_Q<3>::Matrices::n_embedding_matrices = 2; -static double * embedding2[] = +const double * const +FE_Q<3>::Matrices::embedding[][GeometryInfo<3>::children_per_cell] = { - q1_into_q1_refined_2, - q2_into_q2_refined_2, -// q3_into_q3_refined_2, -// q4_into_q4_refined_2, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 + { q1_into_q1_refined_0, q1_into_q1_refined_1, q1_into_q1_refined_2, q1_into_q1_refined_3, + q1_into_q1_refined_4, q1_into_q1_refined_5, q1_into_q1_refined_6, q1_into_q1_refined_7 }, + { q2_into_q2_refined_0, q2_into_q2_refined_1, q2_into_q2_refined_2, q2_into_q2_refined_3, + q2_into_q2_refined_4, q2_into_q2_refined_5, q2_into_q2_refined_6, q2_into_q2_refined_7 } }; -static double * embedding3[] = -{ - q1_into_q1_refined_3, - q2_into_q2_refined_3, -// q3_into_q3_refined_3, -// q4_into_q4_refined_3, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; - -static double * embedding4[] = -{ - q1_into_q1_refined_4, - q2_into_q2_refined_4, -// q3_into_q3_refined_4, -// q4_into_q4_refined_4, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; - -static double * embedding5[] = -{ - q1_into_q1_refined_5, - q2_into_q2_refined_5, -// q3_into_q3_refined_5, -// q4_into_q4_refined_5, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; - -static double * embedding6[] = -{ - q1_into_q1_refined_6, - q2_into_q2_refined_6, -// q3_into_q3_refined_6, -// q4_into_q4_refined_6, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; - -static double * embedding7[] = -{ - q1_into_q1_refined_7, - q2_into_q2_refined_7, -// q3_into_q3_refined_7, -// q4_into_q4_refined_7, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -}; +// Constraint matrices taken from Wolfgangs old version +const unsigned int +FE_Q<3>::Matrices::n_constraint_matrices = 2; -static double constraint_q1[] = +static const double constraint_q1[] = { .25,.25,.25,.25, .5,.5,0.,0., @@ -508,7 +389,7 @@ static double constraint_q1[] = .5,0.,0.,.5 }; -static double constraint_q2[] = +static const double constraint_q2[] = { 0,0,0,0,0,0,0,0,1, 0,0,0,0,1,0,0,0,0, @@ -533,20 +414,10 @@ static double constraint_q2[] = -.046875,.015625,-.046875,.140625,-.09375,-.09375,.28125,.28125,.5625 }; -static double* constraint_matrices[] = +const double * const +FE_Q<3>::Matrices::constraint_matrices[] = { constraint_q1, - constraint_q2, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0, - 0 -// constraint_q3, -// constraint_q4 + constraint_q2 };