From 2654e248ad8af21182cbcd03366c2902195fa2b0 Mon Sep 17 00:00:00 2001 From: hartmann Date: Tue, 23 Dec 2008 14:09:12 +0000 Subject: [PATCH] Bring patch deal.II_DLR@5588 to deal.II (leicht): Supply anisotropic restriction and embedding matrices for FE_Q. git-svn-id: https://svn.dealii.org/trunk@18004 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/source/fe/fe_q.cc | 329 +++++++++++++++--------------- 1 file changed, 162 insertions(+), 167 deletions(-) diff --git a/deal.II/deal.II/source/fe/fe_q.cc b/deal.II/deal.II/source/fe/fe_q.cc index 6c1790dd48..e3fa057a36 100644 --- a/deal.II/deal.II/source/fe/fe_q.cc +++ b/deal.II/deal.II/source/fe/fe_q.cc @@ -413,7 +413,9 @@ struct FE_Q::Implementation for (unsigned int i=0, iy=1; iy<=n; ++iy) for (unsigned int ix=1; ix<=n; ++ix) inner_points[i++] = Point (ix*step, iy*step); - + + // at the moment do this for + // isotropic face refinement only for (unsigned int child=0; child::max_children_per_cell; ++child) for (unsigned int i=0; i::FE_Q (const unsigned int degree) // compute constraint, embedding // and restriction matrices initialize_constraints (); + this->reinit_restriction_and_prolongation_matrices(); initialize_embedding (); initialize_restriction (); @@ -1514,8 +1517,6 @@ template void FE_Q::initialize_embedding () { - unsigned int iso=RefinementCase::isotropic_refinement-1; - // compute the interpolation // matrices in much the same way as // we do for the constraints. it's @@ -1528,126 +1529,124 @@ FE_Q::initialize_embedding () this->dofs_per_cell); const std::vector &index_map= this->poly_space.get_numbering(); - - for (unsigned int child=0; child::max_children_per_cell; ++child) - this->prolongation[iso][child].reinit (this->dofs_per_cell, - this->dofs_per_cell); - for (unsigned int child=0; child::max_children_per_cell; ++child) - { - for (unsigned int j=0; jdofs_per_cell; ++j) - { - // generate a point on - // the child cell and - // evaluate the shape - // functions there - const Point p_subcell - = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, - dealii::internal::int2type()); - const Point p_cell = - GeometryInfo::child_to_cell_coordinates (p_subcell, child); - - for (unsigned int i=0; idofs_per_cell; ++i) - { - const double - cell_value = this->poly_space.compute_value (i, p_cell), - subcell_value = this->poly_space.compute_value (i, p_subcell); - - // cut off values that - // are too small. note - // that we have here - // Lagrange - // interpolation - // functions, so they - // should be zero at - // almost all points, - // and one at the - // others, at least on - // the subcells. so set - // them to their exact - // values - // - // the actual cut-off - // value is somewhat - // fuzzy, but it works - // for - // 1e-14*degree*dim, - // which is kind of - // reasonable given - // that we compute the - // values of the - // polynomials via an - // degree-step - // recursion and then - // multiply the - // 1d-values. this - // gives us a linear - // growth in - // degree*dim, times a - // small constant. - if (std::fabs(cell_value) < 2e-13*this->degree*this->degree*dim) - cell_interpolation(j, i) = 0.; - else - cell_interpolation(j, i) = cell_value; - - if (std::fabs(subcell_value) < 2e-13*this->degree*this->degree*dim) - subcell_interpolation(j, i) = 0.; - else - if (std::fabs(subcell_value-1) < 2e-13*this->degree*this->degree*dim) - subcell_interpolation(j, i) = 1.; - else - // we have put our - // evaluation - // points onto the - // interpolation - // points, so we - // should either - // get zeros or - // ones! - Assert (false, ExcInternalError()); - } - } - // then compute the embedding - // matrix for this child and - // this coordinate - // direction. by the same trick - // as with the constraint - // matrices, don't compute the - // inverse of - // subcell_interpolation, but - // use the fact that we have - // put our interpolation points - // onto the interpolation - // points of the Lagrange - // polynomials used here. then, - // the subcell_interpolation - // matrix is just a permutation - // of the identity matrix and - // its inverse is also its - // transpose - subcell_interpolation.Tmmult (this->prolongation[iso][child], - cell_interpolation); + for (unsigned int ref=0; ref::isotropic_refinement; ++ref) + for (unsigned int child=0; child::n_children(RefinementCase(ref+1)); ++child) + { + for (unsigned int j=0; jdofs_per_cell; ++j) + { + // generate a point on + // the child cell and + // evaluate the shape + // functions there + const Point p_subcell + = FE_Q_Helper::generate_unit_point (index_map[j], this->dofs_per_cell, + dealii::internal::int2type()); + const Point p_cell = + GeometryInfo::child_to_cell_coordinates (p_subcell, child, RefinementCase(ref+1)); + + for (unsigned int i=0; idofs_per_cell; ++i) + { + const double + cell_value = this->poly_space.compute_value (i, p_cell), + subcell_value = this->poly_space.compute_value (i, p_subcell); + + // cut off values that + // are too small. note + // that we have here + // Lagrange + // interpolation + // functions, so they + // should be zero at + // almost all points, + // and one at the + // others, at least on + // the subcells. so set + // them to their exact + // values + // + // the actual cut-off + // value is somewhat + // fuzzy, but it works + // for + // 1e-14*degree*dim, + // which is kind of + // reasonable given + // that we compute the + // values of the + // polynomials via an + // degree-step + // recursion and then + // multiply the + // 1d-values. this + // gives us a linear + // growth in + // degree*dim, times a + // small constant. + if (std::fabs(cell_value) < 2e-13*this->degree*this->degree*dim) + cell_interpolation(j, i) = 0.; + else + cell_interpolation(j, i) = cell_value; + + if (std::fabs(subcell_value) < 2e-13*this->degree*this->degree*dim) + subcell_interpolation(j, i) = 0.; + else + if (std::fabs(subcell_value-1) < 2e-13*this->degree*this->degree*dim) + subcell_interpolation(j, i) = 1.; + else + // we have put our + // evaluation + // points onto the + // interpolation + // points, so we + // should either + // get zeros or + // ones! + Assert (false, ExcInternalError()); + } + } + + // then compute the embedding + // matrix for this child and + // this coordinate + // direction. by the same trick + // as with the constraint + // matrices, don't compute the + // inverse of + // subcell_interpolation, but + // use the fact that we have + // put our interpolation points + // onto the interpolation + // points of the Lagrange + // polynomials used here. then, + // the subcell_interpolation + // matrix is just a permutation + // of the identity matrix and + // its inverse is also its + // transpose + subcell_interpolation.Tmmult (this->prolongation[ref][child], + cell_interpolation); // cut off very small values // here - for (unsigned int i=0; idofs_per_cell; ++i) - for (unsigned int j=0; jdofs_per_cell; ++j) - if (std::fabs(this->prolongation[iso][child](i,j)) < 2e-13*this->degree*dim) - this->prolongation[iso][child](i,j) = 0.; - - // and make sure that the row - // sum is 1. this must be so - // since for this element, the - // shape functions add up to on - for (unsigned int row=0; rowdofs_per_cell; ++row) - { - double sum = 0; - for (unsigned int col=0; coldofs_per_cell; ++col) - sum += this->prolongation[iso][child](row,col); - Assert (std::fabs(sum-1.) < 2e-13*this->degree*this->degree*dim, - ExcInternalError()); - } - } + for (unsigned int i=0; idofs_per_cell; ++i) + for (unsigned int j=0; jdofs_per_cell; ++j) + if (std::fabs(this->prolongation[ref][child](i,j)) < 2e-13*this->degree*dim) + this->prolongation[ref][child](i,j) = 0.; + + // and make sure that the row + // sum is 1. this must be so + // since for this element, the + // shape functions add up to on + for (unsigned int row=0; rowdofs_per_cell; ++row) + { + double sum = 0; + for (unsigned int col=0; coldofs_per_cell; ++col) + sum += this->prolongation[ref][child](row,col); + Assert (std::fabs(sum-1.) < 2e-13*this->degree*this->degree*dim, + ExcInternalError()); + } + } } @@ -1656,8 +1655,6 @@ template void FE_Q::initialize_restriction () { - unsigned int iso=RefinementCase::isotropic_refinement-1; - // for these Lagrange interpolation // polynomials, construction of the // restriction matrices is @@ -1703,8 +1700,6 @@ FE_Q::initialize_restriction () // one child) by the same value // (compute on a later child), so // we don't have to care about this - for (unsigned int c=0; c::max_children_per_cell; ++c) - this->restriction[iso][c].reinit (this->dofs_per_cell, this->dofs_per_cell); for (unsigned int i=0; idofs_per_cell; ++i) { const Point p_cell @@ -1736,50 +1731,50 @@ FE_Q::initialize_restriction () // then find the children on // which the interpolation // point is located - for (unsigned int child=0; child::max_children_per_cell; - ++child) - { - // first initialize this - // column of the matrix - for (unsigned int j=0; jdofs_per_cell; ++j) - this->restriction[iso][child](mother_dof, j) = 0.; - - // then check whether this - // interpolation point is - // inside this child cell - const Point p_subcell - = GeometryInfo::cell_to_child_coordinates (p_cell, child); - if (GeometryInfo::is_inside_unit_cell (p_subcell)) - { - // find the one child - // shape function - // corresponding to - // this point. do it in - // the same way as - // above - unsigned int child_dof = 0; - for (; child_dofdofs_per_cell; ++child_dof) - { - const double val - = this->poly_space.compute_value(child_dof, p_subcell); - if (std::fabs (val-1.) < 2e-13*this->degree*this->degree*dim) - break; - else - Assert (std::fabs(val) < 2e-13*this->degree*this->degree*dim, - ExcInternalError()); - } - for (unsigned int j=child_dof+1; jdofs_per_cell; ++j) - Assert (std::fabs (this->poly_space.compute_value(j, p_subcell)) - < 2e-13*this->degree*this->degree*dim, - ExcInternalError()); - - // so now that we have - // it, set the - // corresponding value - // in the matrix - this->restriction[iso][child](mother_dof, child_dof) = 1.; - } - } + for (unsigned int ref=RefinementCase::cut_x; ref<=RefinementCase::isotropic_refinement; ++ref) + for (unsigned int child=0; child::n_children(RefinementCase(ref)); ++child) + { + // first initialize this + // column of the matrix + for (unsigned int j=0; jdofs_per_cell; ++j) + this->restriction[ref-1][child](mother_dof, j) = 0.; + + // then check whether this + // interpolation point is + // inside this child cell + const Point p_subcell + = GeometryInfo::cell_to_child_coordinates (p_cell, child, RefinementCase(ref)); + if (GeometryInfo::is_inside_unit_cell (p_subcell)) + { + // find the one child + // shape function + // corresponding to + // this point. do it in + // the same way as + // above + unsigned int child_dof = 0; + for (; child_dofdofs_per_cell; ++child_dof) + { + const double val + = this->poly_space.compute_value(child_dof, p_subcell); + if (std::fabs (val-1.) < 2e-13*this->degree*this->degree*dim) + break; + else + Assert (std::fabs(val) < 2e-13*this->degree*this->degree*dim, + ExcInternalError()); + } + for (unsigned int j=child_dof+1; jdofs_per_cell; ++j) + Assert (std::fabs (this->poly_space.compute_value(j, p_subcell)) + < 2e-13*this->degree*this->degree*dim, + ExcInternalError()); + + // so now that we have + // it, set the + // corresponding value + // in the matrix + this->restriction[ref-1][child](mother_dof, child_dof) = 1.; + } + } } } -- 2.39.5