From: kronbichler Date: Sun, 21 Apr 2013 20:04:34 +0000 (+0000) Subject: Avoid using FETools::get_fe_from_name and evaluate on 1D line instead. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7d5b72ae69d19f97f29a7d7afecb6f7869671c46;p=dealii-svn.git Avoid using FETools::get_fe_from_name and evaluate on 1D line instead. git-svn-id: https://svn.dealii.org/trunk@29356 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/matrix_free/shape_info.h b/deal.II/include/deal.II/matrix_free/shape_info.h index 438b418254..654db99e48 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.h @@ -46,30 +46,17 @@ namespace internal ShapeInfo (); /** - * Initializes the data fields. Takes a - * one-dimensional quadrature formula and a - * finite element as arguments and evaluates - * the shape functions, gradients and Hessians - * on the one-dimensional unit cell. This - * function assumes that the finite element is - * derived from a one-dimensional element by a - * tensor product. It uses FETools::get_name() - * and FETools::get_fe_from_name() to find the - * one-dimensional element corresponding to - * the input element in @p dim dimensions. + * Initializes the data fields. Takes a one-dimensional quadrature + * formula and a finite element as arguments and evaluates the shape + * functions, gradients and Hessians on the one-dimensional unit + * cell. This function assumes that the finite element is derived from a + * one-dimensional element by a tensor product and that the zeroth shape + * function in zero evaluates to one. */ template void reinit (const Quadrature<1> &quad, const FiniteElement &fe_dim); - /** - * Internal helper function for initialization - * that does the main work. - */ - void do_initialize (const Quadrature<1> &quad, - const FiniteElement<1> &fe, - const unsigned int dim); - /** * Returns the memory consumption of this * class in bytes. diff --git a/deal.II/include/deal.II/matrix_free/shape_info.templates.h b/deal.II/include/deal.II/matrix_free/shape_info.templates.h index 491b1851b2..c1f251e56d 100644 --- a/deal.II/include/deal.II/matrix_free/shape_info.templates.h +++ b/deal.II/include/deal.II/matrix_free/shape_info.templates.h @@ -44,63 +44,37 @@ namespace internal template void ShapeInfo::reinit (const Quadrature<1> &quad, - const FiniteElement &fe_dim) + const FiniteElement &fe) { - Assert (fe_dim.n_components() == 1, + Assert (fe.n_components() == 1, ExcMessage("FEEvaluation only works for scalar finite elements.")); - // take the name of the finite element - // and generate a 1d element. read the - // name, change the template argument - // to one and construct an element - std::string fe_name = fe_dim.get_name(); - const std::size_t template_starts = fe_name.find_first_of('<'); - Assert (fe_name[template_starts+1] == (dim==1?'1':(dim==2?'2':'3')), - ExcInternalError()); - fe_name[template_starts+1] = '1'; - std_cxx1x::shared_ptr > fe_1d - (FETools::get_fe_from_name<1>(fe_name)); - const FiniteElement<1> &fe = *fe_1d; - do_initialize (quad, fe, dim); - } - - - template - void - ShapeInfo::do_initialize (const Quadrature<1> &quad, - const FiniteElement<1> &fe, - const unsigned int dim) - { - const unsigned int n_dofs_1d = fe.dofs_per_cell, + const unsigned int n_dofs_1d = fe.degree+1, n_q_points_1d = quad.size(); - std::vector lexicographic (n_dofs_1d); + AssertDimension(fe.dofs_per_cell, Utilities::fixed_power(n_dofs_1d)); + std::vector lexicographic (fe.dofs_per_cell); - // renumber (this is necessary for FE_Q, for - // example, since there the vertex DoFs come - // first, which is incompatible with the - // lexicographic ordering necessary to apply - // tensor products efficiently) + // renumber (this is necessary for FE_Q, for example, since there the + // vertex DoFs come first, which is incompatible with the lexicographic + // ordering necessary to apply tensor products efficiently) { - const FE_Poly,1,1> *fe_poly = - dynamic_cast,1,1>*>(&fe); + const FE_Poly,dim,dim> *fe_poly = + dynamic_cast,dim,dim>*>(&fe); Assert (fe_poly != 0, ExcNotImplemented()); - lexicographic = fe_poly->get_poly_space_numbering(); + lexicographic = fe_poly->get_poly_space_numbering_inverse(); + + // to evaluate 1D polynomials, evaluate along the line where y=z=0, + // assuming that shape_value(0,Point()) == 1. otherwise, need + // other entry point (e.g. generating a 1D element by reading the + // name, as done before r29356) + Assert(std::fabs(fe.shape_value(lexicographic[0], Point())-1) < 1e-13, + ExcInternalError()); } - n_q_points = 1; - dofs_per_cell = 1; - n_q_points_face = 1; - dofs_per_face = 1; - for (unsigned int d=0; d(dim)-1; ++d) - { - n_q_points_face *= n_q_points_1d; - dofs_per_face *= n_dofs_1d; - } + n_q_points = Utilities::fixed_power(n_q_points_1d); + dofs_per_cell = Utilities::fixed_power(n_dofs_1d); + n_q_points_face = dim>1?Utilities::fixed_power(n_q_points_1d):1; + dofs_per_face = dim>1?Utilities::fixed_power(n_dofs_1d):1; const unsigned int array_size = n_dofs_1d*n_q_points_1d; this->shape_gradients.resize_fast (array_size); @@ -116,8 +90,8 @@ namespace internal for (unsigned int i=0; i::n_array_elements // copies for the shape information and // non-vectorized fields - const Point<1> q_point = quad.get_points()[q]; - shape_values_number[my_i*n_q_points_1d+q] = fe.shape_value(i,q_point); - shape_gradient_number[my_i*n_q_points_1d+q] = fe.shape_grad (i,q_point)[0]; - shape_values [my_i*n_q_points_1d+q] = - shape_values_number [my_i*n_q_points_1d+q]; - shape_gradients[my_i*n_q_points_1d+q] = - shape_gradient_number[my_i*n_q_points_1d+q]; - shape_hessians[my_i*n_q_points_1d+q] = - fe.shape_grad_grad(i,q_point)[0][0]; - face_value[0][my_i*n_q_points_1d+q] = fe.shape_value(i,q_point*0.5); - face_value[1][my_i*n_q_points_1d+q] = fe.shape_value(i,Point<1>(0.5)+q_point*0.5); + Point q_point; + q_point[0] = quad.get_points()[q][0]; + shape_values_number[i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); + shape_gradient_number[i*n_q_points_1d+q] = fe.shape_grad (my_i,q_point)[0]; + shape_values [i*n_q_points_1d+q] = + shape_values_number [i*n_q_points_1d+q]; + shape_gradients[i*n_q_points_1d+q] = + shape_gradient_number[i*n_q_points_1d+q]; + shape_hessians[i*n_q_points_1d+q] = + fe.shape_grad_grad(my_i,q_point)[0][0]; + q_point[0] *= 0.5; + face_value[0][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); + q_point[0] += 0.5; + face_value[1][i*n_q_points_1d+q] = fe.shape_value(my_i,q_point); } - this->face_gradient[0][my_i] = fe.shape_grad(i,Point<1>(0.))[0]; - this->face_gradient[1][my_i] = fe.shape_grad(i,Point<1>(1.))[0]; + Point q_point; + this->face_gradient[0][i] = fe.shape_grad(my_i,q_point)[0]; + q_point[0] = 1; + this->face_gradient[1][i] = fe.shape_grad(my_i,q_point)[0]; } // face information - unsigned int n_faces = 1; - for (unsigned int d=0; d::faces_per_cell; this->face_indices.reinit(n_faces, this->dofs_per_face); switch (dim) {