From 90c438f3cc01e22f11c08b5cf0ad5bfdfdcdb2fa Mon Sep 17 00:00:00 2001 From: Guido Kanschat Date: Mon, 18 Apr 2005 18:24:00 +0000 Subject: [PATCH] use FETools::compute_embedding_matrices git-svn-id: https://svn.dealii.org/trunk@10510 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_dgp_monomial.cc | 51 +------ .../deal.II/source/fe/fe_raviart_thomas.cc | 135 +----------------- 2 files changed, 12 insertions(+), 174 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc index 72de096725..094e8871ea 100644 --- a/deal.II/deal.II/source/fe/fe_dgp_monomial.cc +++ b/deal.II/deal.II/source/fe/fe_dgp_monomial.cc @@ -12,6 +12,7 @@ //--------------------------------------------------------------------------- #include +#include #ifdef HAVE_STD_STRINGSTREAM # include @@ -126,7 +127,10 @@ FE_DGPMonomial::FE_DGPMonomial (const unsigned int degree) // initialize the interpolation // matrices - initialize_embedding (); + for (unsigned int i=0; i::children_per_cell; ++i) + this->prolongation[i].reinit (this->dofs_per_cell, + this->dofs_per_cell); + FETools::compute_embedding_matrices (*this, &this->prolongation[0]); // initialize_restriction (); // note, that these elements have @@ -223,51 +227,6 @@ get_interpolation_matrix (const FiniteElementBase &source_fe, -template -void -FE_DGPMonomial::initialize_embedding () -{ - std::vector > unit_points(this->dofs_per_cell); - generate_unit_points(this->degree, unit_points); - - FullMatrix cell_interpolation (this->dofs_per_cell, - this->dofs_per_cell); - FullMatrix subcell_interpolation (this->dofs_per_cell, - this->dofs_per_cell); - for (unsigned int child=0; child::children_per_cell; ++child) - this->prolongation[child].reinit (this->dofs_per_cell, - this->dofs_per_cell); - for (unsigned int child=0; child::children_per_cell; ++child) - { - for (unsigned int j=0; jdofs_per_cell; ++j) - { - const Point &p_subcell=unit_points[j]; - - const Point p_cell = - GeometryInfo::child_to_cell_coordinates (p_subcell, child); - - for (unsigned int i=0; idofs_per_cell; ++i) - { - cell_interpolation(j,i) = this->poly_space.compute_value (i, p_cell); - subcell_interpolation(j,i) = this->poly_space.compute_value (i, p_subcell); - } - } - - // then compute the embedding - // matrix for this child and - // this coordinate direction - subcell_interpolation.gauss_jordan (); - subcell_interpolation.mmult (this->prolongation[child], cell_interpolation); - - // cut off very small values - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - if (std::fabs(this->prolongation[child](i,j)) < 2e-14*this->degree*dim) - this->prolongation[child](i,j) = 0.; - } -} - - template void FE_DGPMonomial::initialize_restriction () diff --git a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc index a0c8fb9bc2..dda14f93f1 100644 --- a/deal.II/deal.II/source/fe/fe_raviart_thomas.cc +++ b/deal.II/deal.II/source/fe/fe_raviart_thomas.cc @@ -20,6 +20,7 @@ #include #include #include +#include #ifdef HAVE_STD_STRINGSTREAM # include @@ -155,7 +156,12 @@ FE_RaviartThomas::FE_RaviartThomas (const unsigned int degree) // initialize the various matrices initialize_constraints (); - initialize_embedding (); + + for (unsigned int i=0; i::children_per_cell; ++i) + this->prolongation[i].reinit (this->dofs_per_cell, + this->dofs_per_cell); + FETools::compute_embedding_matrices (*this, &this->prolongation[0]); + initialize_restriction (); // finally fill in support points @@ -605,133 +611,6 @@ FE_RaviartThomas<3>::initialize_constraints () #endif -#if deal_II_dimension == 1 - -template <> -void -FE_RaviartThomas<1>::initialize_embedding () -{ - Assert (false, ExcNotUsefulInThisDimension()); -} - -#endif - - -template -void -FE_RaviartThomas::initialize_embedding () -{ - // 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 weird - // renumbering stuff going on - // - // it is, however, slightly - // complicated by the fact that we - // have vector-valued elements - // here, so we do all the stuff for - // the degrees of freedom - // corresponding to each coordinate - // direction separately - const unsigned int dofs_per_coordinate = this->dofs_per_cell/dim; - Assert (dofs_per_coordinate*dim == this->dofs_per_cell, - ExcInternalError()); - for (unsigned int d=0; d cell_interpolation (dofs_per_coordinate, - dofs_per_coordinate); - FullMatrix subcell_interpolation (dofs_per_coordinate, - dofs_per_coordinate); - FullMatrix tmp (dofs_per_coordinate, - dofs_per_coordinate); - for (unsigned int child=0; child::children_per_cell; ++child) - this->prolongation[child].reinit (this->dofs_per_cell, - this->dofs_per_cell); - for (unsigned int child=0; child::children_per_cell; ++child) - for (unsigned int d=0; d p_subcell = generate_unit_point (j, dofs_per_coordinate, - d, int2type()); - const Point p_cell = - GeometryInfo::child_to_cell_coordinates (p_subcell, child); - - for (unsigned int i=0; idofs_per_cell; ++i) - if (renumber[i].first == d) - for (unsigned int j=0; jdofs_per_cell; ++j) - if (renumber[j].first == d) - if (std::fabs(tmp(renumber[i].second,renumber[j].second)) > 1e-15) - this->prolongation[child](i,j) = tmp(renumber[i].second, - renumber[j].second); - } - - // if this were a Lagrange - // interpolation element, we could - // make sure that the row sum of - // each of the matrices is 1 at - // this point. note that this won't - // work here, since we are working - // with hierarchical elements for - // which the shape functions don't - // sum up to 1 - // - // however, we can make sure that - // only components couple that have - // the same vector component - for (unsigned int child=0; child::children_per_cell; ++child) - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - Assert ((this->prolongation[child](i,j) == 0.) || - (renumber[i].first == renumber[j].first), - ExcInternalError()); - - - // there is one additional thing to - // be considered: since the shape - // functions on the real cell - // contain the Jacobian (actually, - // the determinant of the inverse), - // there is an additional factor of - // 2 when going from the big to the - // small cell: - for (unsigned int child=0; child::children_per_cell; ++child) - this->prolongation[child] *= 1./2; -} - - #if deal_II_dimension == 1 template <> -- 2.39.5