From: heister Date: Wed, 19 Dec 2012 17:52:41 +0000 (+0000) Subject: merge from upstream X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d251ee32a0bb0ddd40ca0e513bda8d1142c340f1;p=dealii-svn.git merge from upstream git-svn-id: https://svn.dealii.org/branches/branch_bigger_global_dof_indices_3@27830 0785d39b-7218-0410-832d-ea1e28bc413d --- d251ee32a0bb0ddd40ca0e513bda8d1142c340f1 diff --cc deal.II/include/deal.II/lac/constraint_matrix.h index e665929d2b,103433b068..03143da3cc --- a/deal.II/include/deal.II/lac/constraint_matrix.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.h @@@ -1903,9 -1903,18 +1903,18 @@@ ConstraintMatrix::set_inhomogeneity (co + inline + unsigned int + ConstraintMatrix::n_constraints () const + { + return lines.size(); + } + + + inline bool -ConstraintMatrix::is_constrained (const unsigned int index) const +ConstraintMatrix::is_constrained (const types::global_dof_index index) const { const unsigned int line_index = calculate_line_index(index); return ((line_index < lines_cache.size()) diff --cc deal.II/include/deal.II/lac/sparse_matrix.h index 2d77baf4b0,9092530630..add3c01108 --- a/deal.II/include/deal.II/lac/sparse_matrix.h +++ b/deal.II/include/deal.II/lac/sparse_matrix.h @@@ -610,24 -484,17 +484,17 @@@ public explicit SparseMatrix (const SparsityPattern &sparsity); /** - * Copy constructor: initialize - * the matrix with the identity - * matrix. This constructor will - * throw an exception if the - * sizes of the sparsity pattern - * and the identity matrix do not - * coincide, or if the sparsity - * pattern does not provide for - * nonzero entries on the entire - * diagonal. + * Copy constructor: initialize the matrix with the identity matrix. This + * constructor will throw an exception if the sizes of the sparsity pattern + * and the identity matrix do not coincide, or if the sparsity pattern does + * not provide for nonzero entries on the entire diagonal. */ SparseMatrix (const SparsityPattern &sparsity, - const IdentityMatrix &id); + const IdentityMatrix &id); /** - * Destructor. Free all memory, but do not - * release the memory of the sparsity - * structure. + * Destructor. Free all memory, but do not release the memory of the + * sparsity structure. */ virtual ~SparseMatrix (); @@@ -653,34 -510,21 +510,21 @@@ SparseMatrix &operator = (const SparseMatrix &); /** - * Copy operator: initialize - * the matrix with the identity - * matrix. This operator will - * throw an exception if the - * sizes of the sparsity pattern - * and the identity matrix do not - * coincide, or if the sparsity - * pattern does not provide for - * nonzero entries on the entire - * diagonal. + * Copy operator: initialize the matrix with the identity matrix. This + * operator will throw an exception if the sizes of the sparsity pattern and + * the identity matrix do not coincide, or if the sparsity pattern does not + * provide for nonzero entries on the entire diagonal. */ SparseMatrix & - operator= (const IdentityMatrix &id); + operator= (const IdentityMatrix &id); /** - * This operator assigns a scalar to - * a matrix. Since this does usually - * not make much sense (should we set - * all matrix entries to this value? - * Only the nonzero entries of the - * sparsity pattern?), this operation - * is only allowed if the actual - * value to be assigned is zero. This - * operator only exists to allow for - * the obvious notation - * matrix=0, which sets all - * elements of the matrix to zero, - * but keep the sparsity pattern + * This operator assigns a scalar to a matrix. Since this does usually not + * make much sense (should we set all matrix entries to this value? Only + * the nonzero entries of the sparsity pattern?), this operation is only + * allowed if the actual value to be assigned is zero. This operator only + * exists to allow for the obvious notation matrix=0, which sets + * all elements of the matrix to zero, but keep the sparsity pattern * previously used. */ SparseMatrix &operator = (const double d); diff --cc deal.II/include/deal.II/lac/sparsity_pattern.h index 0eb8828807,fe87da3f19..5c13c278e0 --- a/deal.II/include/deal.II/lac/sparsity_pattern.h +++ b/deal.II/include/deal.II/lac/sparsity_pattern.h @@@ -496,46 -385,26 +385,26 @@@ public /** * Make a copy with extra off-diagonals. * - * This constructs objects intended for - * the application of the ILU(n)-method - * or other incomplete decompositions. - * Therefore, additional to the original - * entry structure, space for - * extra_off_diagonals - * side-diagonals is provided on both - * sides of the main diagonal. - * - * max_per_row is the - * maximum number of nonzero - * elements per row which this - * structure is to hold. It is - * assumed that this number is - * sufficiently large to - * accommodate both the elements - * in original as well - * as the new off-diagonal - * elements created by this - * constructor. You will usually - * want to give the same number - * as you gave for - * original plus the - * number of side diagonals times - * two. You may however give a - * larger value if you wish to - * add further nonzero entries - * for the decomposition based on - * other criteria than their - * being on side-diagonals. - * - * This function requires that - * original refers to a - * quadratic matrix structure. - * It must be compressed. The - * matrix structure is not - * compressed after this function - * finishes. + * This constructs objects intended for the application of the ILU(n)-method + * or other incomplete decompositions. Therefore, additional to the + * original entry structure, space for extra_off_diagonals + * side-diagonals is provided on both sides of the main diagonal. + * + * max_per_row is the maximum number of nonzero elements per row + * which this structure is to hold. It is assumed that this number is + * sufficiently large to accommodate both the elements in original + * as well as the new off-diagonal elements created by this constructor. You + * will usually want to give the same number as you gave for + * original plus the number of side diagonals times two. You may + * however give a larger value if you wish to add further nonzero entries + * for the decomposition based on other criteria than their being on + * side-diagonals. + * + * This function requires that original refers to a quadratic + * matrix structure. It must be compressed. The matrix structure is not + * compressed after this function finishes. */ - SparsityPattern (const SparsityPattern &original, + SparsityPattern (const SparsityPattern &original, const unsigned int max_per_row, const unsigned int extra_off_diagonals); diff --cc deal.II/source/dofs/dof_tools.cc index 30697f7ee0,619d478814..fcfea2260e --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@@ -382,19 -367,14 +367,14 @@@ namespace DoFTool }; #endif - std::vector dofs_on_this_face; + std::vector dofs_on_this_face; dofs_on_this_face.reserve (max_dofs_per_face(dof)); - // loop over all faces to check - // whether they are at a - // boundary. note that we need not - // take special care of single - // lines (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every boundary line is - // also part of a boundary face. + // loop over all faces to check whether they are at a boundary. note + // that we need not take special care of single lines (using + // @p{cell->has_boundary_lines}), since we do not support boundaries of + // dimension dim-2, and so every boundary line is also part of a + // boundary face. typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); for (; cell!=endc; ++cell) @@@ -445,10 -422,9 +422,9 @@@ cell = cell->child(direction); const unsigned int dofs_per_vertex = cell->get_fe().dofs_per_vertex; - std::vector boundary_dof_boundary_indices (dofs_per_vertex); + std::vector boundary_dof_boundary_indices (dofs_per_vertex); - // next get boundary mapped dof - // indices of boundary dofs + // next get boundary mapped dof indices of boundary dofs for (unsigned int i=0; ivertex_dof_index(direction,i)]; @@@ -1217,10 -1172,9 +1172,9 @@@ { namespace { - inline - bool + inline bool check_master_dof_list (const FullMatrix &face_interpolation_matrix, - const std::vector &master_dof_list) + const std::vector &master_dof_list) { const unsigned int N = master_dof_list.size(); @@@ -1387,43 -1302,28 +1302,28 @@@ (fe2.dofs_per_quad <= fe1.dofs_per_quad), ExcInternalError()); - // the idea here is to designate as - // many DoFs in fe1 per object - // (vertex, line, quad) as master as - // there are such dofs in fe2 - // (indices are int, because we want - // to avoid the 'unsigned int < 0 is - // always false warning for the cases - // at the bottom in 1d and 2d) + // the idea here is to designate as many DoFs in fe1 per object + // (vertex, line, quad) as master as there are such dofs in fe2 + // (indices are int, because we want to avoid the 'unsigned int < 0 + // is always false warning for the cases at the bottom in 1d and + // 2d) // - // as mentioned in the paper, it is - // not always easy to find a set of - // master dofs that produces an - // invertible matrix. to this end, we - // check in each step whether the - // matrix is still invertible and - // simply discard this dof if the - // matrix is not invertible anymore. + // as mentioned in the paper, it is not always easy to find a set + // of master dofs that produces an invertible matrix. to this end, + // we check in each step whether the matrix is still invertible and + // simply discard this dof if the matrix is not invertible anymore. // - // the cases where we did have - // trouble in the past were with - // adding more quad dofs when Q3 and - // Q4 elements meet at a refined face - // in 3d (see the hp/crash_12 test - // that tests that we can do exactly - // this, and failed before we had - // code to compensate for this - // case). the other case are system - // elements: if we have say a Q1Q2 vs - // a Q2Q3 element, then we can't just - // take all master dofs on a line - // from a single base element, since - // the shape functions of that base - // element are independent of that of - // the other one. this latter case - // shows up when running - // hp/hp_constraints_q_system_06 + // the cases where we did have trouble in the past were with adding + // more quad dofs when Q3 and Q4 elements meet at a refined face in + // 3d (see the hp/crash_12 test that tests that we can do exactly + // this, and failed before we had code to compensate for this + // case). the other case are system elements: if we have say a Q1Q2 + // vs a Q2Q3 element, then we can't just take all master dofs on a + // line from a single base element, since the shape functions of + // that base element are independent of that of the other one. this + // latter case shows up when running hp/hp_constraints_q_system_06 + - std::vector master_dof_list; + std::vector master_dof_list; unsigned int index = 0; for (int v=0; v(GeometryInfo::vertices_per_face); @@@ -1523,10 -1411,9 +1411,9 @@@ AssertDimension (index, fe1.dofs_per_face); AssertDimension (master_dof_list.size(), fe2.dofs_per_face); - // finally copy the list into the - // mask + // finally copy the list into the mask std::fill (master_dof_mask.begin(), master_dof_mask.end(), false); - for (std::vector::const_iterator i=master_dof_list.begin(); + for (std::vector::const_iterator i=master_dof_list.begin(); i!=master_dof_list.end(); ++i) master_dof_mask[*i] = true; } @@@ -1802,24 -1661,18 +1661,18 @@@ /** - * Copy constraints into a constraint - * matrix object. + * Copy constraints into a constraint matrix object. * - * This function removes zero - * constraints and those, which - * constrain a DoF which was - * already eliminated in one of - * the previous steps of the hp - * hanging node procedure. + * This function removes zero constraints and those, which constrain + * a DoF which was already eliminated in one of the previous steps of + * the hp hanging node procedure. * - * It also suppresses very small - * entries in the constraint matrix to - * avoid making the sparsity pattern - * fuller than necessary. + * It also suppresses very small entries in the constraint matrix to + * avoid making the sparsity pattern fuller than necessary. */ void - filter_constraints (const std::vector &master_dofs, - const std::vector &slave_dofs, + filter_constraints (const std::vector &master_dofs, + const std::vector &slave_dofs, const FullMatrix &face_constraints, ConstraintMatrix &constraints) { @@@ -1947,28 -1781,26 +1781,26 @@@ void - make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/, - ConstraintMatrix & /*constraints*/) + make_hp_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, + ConstraintMatrix &/*constraints*/) { - // we may have to compute - // constraints for - // vertices. gotta think about + // we may have to compute constraints for vertices. gotta think about // that a bit more - //TODO[WB]: think about what to do here... + + //TODO[WB]: think about what to do here... } void - make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> & /*dof_handler*/, - ConstraintMatrix & /*constraints*/, + make_oldstyle_hanging_node_constraints (const dealii::hp::DoFHandler<1> &/*dof_handler*/, + ConstraintMatrix &/*constraints*/, dealii::internal::int2type<1>) { - // we may have to compute - // constraints for - // vertices. gotta think about + // we may have to compute constraints for vertices. gotta think about // that a bit more - //TODO[WB]: think about what to do here... + + //TODO[WB]: think about what to do here... } @@@ -2057,25 -1886,18 +1886,18 @@@ const unsigned int spacedim = DH::space_dimension; - std::vector dofs_on_mother; - std::vector dofs_on_children; + std::vector dofs_on_mother; + std::vector dofs_on_children; - // loop over all lines; only on - // lines there can be constraints. - // We do so by looping over all - // active cells and checking - // whether any of the faces are - // refined which can only be from - // the neighboring cell because - // this one is active. In that - // case, the face is subject to + // loop over all lines; only on lines there can be constraints. We do + // so by looping over all active cells and checking whether any of + // the faces are refined which can only be from the neighboring cell + // because this one is active. In that case, the face is subject to // constraints // - // note that even though we may - // visit a face twice if the - // neighboring cells are equally - // refined, we can only visit each - // face with hanging nodes once + // note that even though we may visit a face twice if the neighboring + // cells are equally refined, we can only visit each face with + // hanging nodes once typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@@ -2211,25 -2010,18 +2010,18 @@@ { const unsigned int dim = 3; - std::vector dofs_on_mother; - std::vector dofs_on_children; + std::vector dofs_on_mother; + std::vector dofs_on_children; - // loop over all quads; only on - // quads there can be constraints. - // We do so by looping over all - // active cells and checking - // whether any of the faces are - // refined which can only be from - // the neighboring cell because - // this one is active. In that - // case, the face is subject to + // loop over all quads; only on quads there can be constraints. We do + // so by looping over all active cells and checking whether any of + // the faces are refined which can only be from the neighboring cell + // because this one is active. In that case, the face is subject to // constraints // - // note that even though we may - // visit a face twice if the - // neighboring cells are equally - // refined, we can only visit each - // face with hanging nodes once + // note that even though we may visit a face twice if the neighboring + // cells are equally refined, we can only visit each face with + // hanging nodes once typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@@ -2459,30 -2221,21 +2221,21 @@@ const unsigned int spacedim = DH::space_dimension; - // a matrix to be used for - // constraints below. declared - // here and simply resized down - // below to avoid permanent - // re-allocation of memory + // a matrix to be used for constraints below. declared here and + // simply resized down below to avoid permanent re-allocation of + // memory FullMatrix constraint_matrix; - // similarly have arrays that - // will hold master and slave - // dof numbers, as well as a - // scratch array needed for the - // complicated case below + // similarly have arrays that will hold master and slave dof numbers, + // as well as a scratch array needed for the complicated case below - std::vector master_dofs; - std::vector slave_dofs; - std::vector scratch_dofs; + std::vector master_dofs; + std::vector slave_dofs; + std::vector scratch_dofs; - // caches for the face and - // subface interpolation - // matrices between different - // (or the same) finite - // elements. we compute them - // only once, namely the first - // time they are needed, and - // then just reuse them + // caches for the face and subface interpolation matrices between + // different (or the same) finite elements. we compute them only + // once, namely the first time they are needed, and then just reuse + // them Table<2,std_cxx1x::shared_ptr > > face_interpolation_matrices (n_finite_elements (dof_handler), n_finite_elements (dof_handler)); @@@ -3836,9 -3318,8 +3318,8 @@@ ExcInternalError()); } - // then loop over all cells and do - // the work + // then loop over all cells and do the work - std::vector indices; + std::vector indices; for (typename DH::active_cell_iterator c=dof.begin_active(); c!=dof.end(); ++ c) if (c->is_locally_owned()) @@@ -4112,10 -3573,9 +3573,9 @@@ for (unsigned int i=0; i indices(fe.dofs_per_cell); - typename MGDoFHandler::cell_iterator c; + std::vector indices(fe.dofs_per_cell); + typename DH::cell_iterator c; for (c = dof.begin(level) ; c != dof.end(level) ; ++ c) { c->get_mg_dof_indices(indices); @@@ -4199,26 -3656,20 +3656,20 @@@ (component_mask.n_selected_components(n_components(dof_handler)) != n_components(dof_handler))); - std::vector face_dof_indices; + std::vector face_dof_indices; face_dof_indices.reserve (max_dofs_per_face(dof_handler)); - // now loop over all cells and - // check whether their faces are at - // the boundary. note that we need - // not take special care of single - // lines being at the boundary - // (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every isolated boundary - // line is also part of a boundary - // face which we will be visiting - // sooner or later + // now loop over all cells and check whether their faces are at the + // boundary. note that we need not take special care of single lines + // being at the boundary (using @p{cell->has_boundary_lines}), since we + // do not support boundaries of dimension dim-2, and so every isolated + // boundary line is also part of a boundary face which we will be + // visiting sooner or later for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) - // only work on cells that are either locally owned or at - // least ghost cells + + // only work on cells that are either locally owned or at least ghost + // cells if (cell->is_artificial() == false) for (unsigned int face=0; face::faces_per_cell; ++face) @@@ -4318,26 -3755,18 +3755,18 @@@ const bool check_vector_component = (component_mask.represents_the_all_selected_mask() == false); - // clear and reset array by default - // values + // clear and reset array by default values selected_dofs.clear (); selected_dofs.resize (dof_handler.n_dofs(), false); - std::vector cell_dof_indices; + std::vector cell_dof_indices; cell_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - // now loop over all cells and - // check whether their faces are at - // the boundary. note that we need - // not take special care of single - // lines being at the boundary - // (using - // @p{cell->has_boundary_lines}), - // since we do not support - // boundaries of dimension dim-2, - // and so every isolated boundary - // line is also part of a boundary - // face which we will be visiting - // sooner or later + // now loop over all cells and check whether their faces are at the + // boundary. note that we need not take special care of single lines + // being at the boundary (using @p{cell->has_boundary_lines}), since we + // do not support boundaries of dimension dim-2, and so every isolated + // boundary line is also part of a boundary face which we will be + // visiting sooner or later for (typename DH::active_cell_iterator cell=dof_handler.begin_active(); cell!=dof_handler.end(); ++cell) for (unsigned int face=0; @@@ -4537,12 -3955,11 +3955,11 @@@ // preset all values by false std::fill_n (selected_dofs.begin(), dof_handler.n_dofs(), false); - std::vector local_dof_indices; + std::vector local_dof_indices; local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - // this function is similar to the - // make_sparsity_pattern function, - // see there for more information + // this function is similar to the make_sparsity_pattern function, see + // there for more information typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@@ -4579,14 -3996,11 +3996,11 @@@ // collect all the locally owned dofs dof_set = dof_handler.locally_owned_dofs(); - // add the DoF on the adjacent ghost cells - // to the IndexSet, cache them in a - // set. need to check each dof manually - // because we can't be sure that the dof - // range of locally_owned_dofs is really - // contiguous. + // add the DoF on the adjacent ghost cells to the IndexSet, cache them + // in a set. need to check each dof manually because we can't be sure + // that the dof range of locally_owned_dofs is really contiguous. - std::vector dof_indices; - std::set global_dof_indices; + std::vector dof_indices; + std::set global_dof_indices; typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@@ -4618,14 -4032,11 +4032,11 @@@ // collect all the locally owned dofs dof_set = dof_handler.locally_owned_dofs(); - // add the DoF on the adjacent ghost cells - // to the IndexSet, cache them in a - // set. need to check each dof manually - // because we can't be sure that the dof - // range of locally_owned_dofs is really - // contiguous. + // add the DoF on the adjacent ghost cells to the IndexSet, cache them + // in a set. need to check each dof manually because we can't be sure + // that the dof range of locally_owned_dofs is really contiguous. - std::vector dof_indices; - std::set global_dof_indices; + std::vector dof_indices; + std::set global_dof_indices; typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@@ -4733,19 -4142,15 +4142,15 @@@ std::fill_n (subdomain_association.begin(), dof_handler.n_dofs(), types::invalid_subdomain_id); - std::vector local_dof_indices; + std::vector local_dof_indices; local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - // pseudo-randomly assign variables - // which lie on the interface - // between subdomains to each of - // the two or more + // pseudo-randomly assign variables which lie on the interface between + // subdomains to each of the two or more bool coin_flip = true; - // loop over all cells and record - // which subdomain a DoF belongs - // to. toss a coin in case it is on - // an interface + // loop over all cells and record which subdomain a DoF belongs to. + // toss a coin in case it is on an interface typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@@ -4824,18 -4223,15 +4223,15 @@@ IndexSet index_set (dof_handler.n_dofs()); - std::vector local_dof_indices; + std::vector local_dof_indices; local_dof_indices.reserve (max_dofs_per_cell(dof_handler)); - // first generate an unsorted list of all - // indices which we fill from the back. could - // also insert them directly into the - // IndexSet, but that inserts indices in the - // middle, which is an O(n^2) algorithm and - // hence too expensive. Could also use - // std::set, but that is in general more - // expensive than a vector + // first generate an unsorted list of all indices which we fill from + // the back. could also insert them directly into the IndexSet, but + // that inserts indices in the middle, which is an O(n^2) algorithm and + // hence too expensive. Could also use std::set, but that is in general + // more expensive than a vector - std::vector subdomain_indices; + std::vector subdomain_indices; typename DH::active_cell_iterator cell = dof_handler.begin_active(), @@@ -4965,14 -4361,14 +4361,14 @@@ void resolve_components (const hp::FECollection &fe_collection, const std::vector &dofs_by_component, - const std::vector &target_component, + const std::vector &target_component, const bool only_once, - std::vector &dofs_per_component, + std::vector &dofs_per_component, unsigned int &component) { - // assert that all elements in the collection have the same - // structure (base elements and multiplicity, components per base - // element) and then simply call the function above + // assert that all elements in the collection have the same structure + // (base elements and multiplicity, components per base element) and + // then simply call the function above for (unsigned int fe=1; fe dof_is_interesting (fine_grid.n_dofs(), false); - std::vector local_dof_indices (fine_fe.dofs_per_cell); + std::vector local_dof_indices (fine_fe.dofs_per_cell); for (typename dealii::DoFHandler::active_cell_iterator cell=fine_grid.begin_active(); diff --cc deal.II/source/fe/fe_values.cc index a1f22a30b2,5a63aa75f2..028fccb660 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@@ -2036,85 -2371,34 +2371,34 @@@ template void FEValuesBase::get_function_values ( const InputVector &fe_function, - const VectorSlice > &indices, + const VectorSlice > &indices, std::vector > &values) const { - // One value per quadrature point - Assert (n_quadrature_points == values.size(), - ExcDimensionMismatch(values.size(), n_quadrature_points)); - - const unsigned int n_components = fe->n_components(); - - // Size of indices must be a - // multiple of dofs_per_cell such - // that an integer number of - // function values is generated in - // each point. + // Size of indices must be a multiple of dofs_per_cell such that an integer + // number of function values is generated in each point. Assert (indices.size() % dofs_per_cell == 0, ExcNotMultiple(indices.size(), dofs_per_cell)); - - // The number of components of the - // result may be a multiple of the - // number of components of the - // finite element - const unsigned int result_components = indices.size() * n_components / dofs_per_cell; - - for (unsigned int i=0; iupdate_flags & update_values, ExcAccessToUninitializedField()); - // initialize with zero - for (unsigned int i=0; iis_primitive(shape_func)) - { - const unsigned int comp = fe->system_to_component_index(shape_func).first - + mc * n_components; - const unsigned int - row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp]; - - const double *shape_value_ptr = &this->shape_values(row, 0); - for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) - continue; - - const unsigned int - row = this->shape_function_to_row_table[shape_func * fe->n_components() + c]; - - const double *shape_value_ptr = &this->shape_values(row, 0); - const unsigned int comp = c + mc * n_components; - - for (unsigned int point=0; point > > val(values); + if (indices.size() <= 100) + { + double dof_values[100]; + for (unsigned int i=0; ishape_values, *fe, + this->shape_function_to_row_table, val, + false, indices.size()/dofs_per_cell); + } + else + { + Vector dof_values(100); + for (unsigned int i=0; ishape_values, *fe, + this->shape_function_to_row_table, val, + false, indices.size()/dofs_per_cell); + } } @@@ -2587,41 -2581,28 +2581,28 @@@ template void FEValuesBase::get_function_hessians ( const InputVector &fe_function, - const VectorSlice > &indices, + const VectorSlice > &indices, std::vector > &hessians) const { - Assert (this->update_flags & update_second_derivatives, ExcAccessToUninitializedField()); - // This function fills a single - // component only - Assert (fe->n_components() == 1, - ExcDimensionMismatch(fe->n_components(), 1)); - // One index for each dof - Assert (indices.size() == dofs_per_cell, - ExcDimensionMismatch(indices.size(), dofs_per_cell)); - // This vector has one entry for - // each quadrature point - Assert (hessians.size() == n_quadrature_points, - ExcDimensionMismatch(hessians.size(), n_quadrature_points)); - - // initialize with zero - std::fill_n (hessians.begin(), n_quadrature_points, Tensor<2,spacedim>()); - - // add up contributions of trial - // functions. note that here we - // deal with scalar finite - // elements, so no need to check - // for non-primitivity of shape - // functions - for (unsigned int shape_func=0; shape_funcupdate_flags & update_second_derivatives, + ExcAccessToUninitializedField()); + AssertDimension (fe_function.size(), present_cell->n_dofs_for_dof_handler()); + AssertDimension (indices.size(), dofs_per_cell); + if (dofs_per_cell <= 100) + { + double dof_values[100]; + for (unsigned int i=0; ishape_hessians, + hessians); + } + else { - const double value = get_vector_element (fe_function, indices[shape_func]); - if (value == 0.) - continue; - - const Tensor<2,spacedim> *shape_hessians_ptr - = &this->shape_hessians[shape_func][0]; - for (unsigned int point=0; point dof_values(dofs_per_cell); + for (unsigned int i=0; ishape_hessians, + hessians); } } @@@ -3000,87 -2745,34 +2745,34 @@@ template void FEValuesBase::get_function_laplacians ( const InputVector &fe_function, - const VectorSlice > &indices, + const VectorSlice > &indices, std::vector > &laplacians) const { - // One value per quadrature point - Assert (n_quadrature_points == laplacians.size(), - ExcDimensionMismatch(laplacians.size(), n_quadrature_points)); - - const unsigned int n_components = fe->n_components(); - - // Size of indices must be a - // multiple of dofs_per_cell such - // that an integer number of - // function values is generated in - // each point. + // Size of indices must be a multiple of dofs_per_cell such that an integer + // number of function values is generated in each point. Assert (indices.size() % dofs_per_cell == 0, ExcNotMultiple(indices.size(), dofs_per_cell)); - - // The number of components of the - // result may be a multiple of the - // number of components of the - // finite element - const unsigned int result_components = indices.size() * n_components / dofs_per_cell; - - for (unsigned int i=0; iupdate_flags & update_hessians, ExcAccessToUninitializedField()); - - // initialize with zero - for (unsigned int i=0; iis_primitive(shape_func)) - { - const unsigned int comp = fe->system_to_component_index(shape_func).first - + mc * n_components; - const unsigned int - row = this->shape_function_to_row_table[shape_func * fe->n_components() + comp]; - - const Tensor<2,spacedim> *shape_hessian_ptr - = &this->shape_hessians[row][0]; - for (unsigned int point=0; pointget_nonzero_components(shape_func)[c] == false) - continue; - - const unsigned int - row = this->shape_function_to_row_table[shape_func * fe->n_components() + c]; - - const Tensor<2,spacedim> *shape_hessian_ptr - = &this->shape_hessians[row][0]; - const unsigned int comp = c + mc * n_components; - - for (unsigned int point=0; pointshape_hessians, + *fe, this->shape_function_to_row_table, + laplacians, false, + indices.size()/dofs_per_cell); + } + else + { + Vector dof_values(indices.size()); + for (unsigned int i=0; ishape_hessians, + *fe, this->shape_function_to_row_table, + laplacians, false, + indices.size()/dofs_per_cell); + } }