From: maier Date: Mon, 3 Dec 2012 22:58:22 +0000 (+0000) Subject: Deuglify X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=0aed7c27ece33c382f8bf72cdb9c18585a512ee5;p=dealii-svn.git Deuglify Also, try to find a fixpoint (*) for initializer lists under astyle (*) as usual, with the help of a fixpoint iteration... git-svn-id: https://svn.dealii.org/trunk@27742 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/grid/tria_objects.h b/deal.II/include/deal.II/grid/tria_objects.h index a0cf6f1a85..c90e306a62 100644 --- a/deal.II/include/deal.II/grid/tria_objects.h +++ b/deal.II/include/deal.II/grid/tria_objects.h @@ -763,7 +763,7 @@ namespace internal template inline - void*& + void *& TriaObjects::user_pointer (const unsigned int i) { #ifdef DEBUG diff --git a/deal.II/source/dofs/dof_tools.cc b/deal.II/source/dofs/dof_tools.cc index aecaa3dd48..02d1a637de 100644 --- a/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/source/dofs/dof_tools.cc @@ -69,10 +69,9 @@ namespace DoFTools Assert (sparsity.n_cols() == n_dofs, ExcDimensionMismatch (sparsity.n_cols(), n_dofs)); - // If we have a distributed::Triangulation only - // allow locally_owned subdomain. Not setting a - // subdomain is also okay, because we skip ghost - // cells in the loop below. + // If we have a distributed::Triangulation only allow locally_owned + // subdomain. Not setting a subdomain is also okay, because we skip + // ghost cells in the loop below. Assert ( (dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id) || @@ -88,11 +87,9 @@ namespace DoFTools typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); - // In case we work with a distributed - // sparsity pattern of Trilinos type, we - // only have to do the work if the - // current cell is owned by the calling - // processor. Otherwise, just continue. + // In case we work with a distributed sparsity pattern of Trilinos + // type, we only have to do the work if the current cell is owned by + // the calling processor. Otherwise, just continue. for (; cell!=endc; ++cell) if (((subdomain_id == types::invalid_subdomain_id) || @@ -104,10 +101,9 @@ namespace DoFTools dofs_on_this_cell.resize (dofs_per_cell); cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed + // make sparsity pattern for this cell. if no constraints pattern + // was given, then the following call acts as if simply no + // constraints existed constraints.add_entries_local_to_global (dofs_on_this_cell, sparsity, keep_constrained_dofs); @@ -136,10 +132,9 @@ namespace DoFTools Assert (couplings.n_cols() == dof.get_fe().n_components(), ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); - // If we have a distributed::Triangulation only - // allow locally_owned subdomain. Not setting a - // subdomain is also okay, because we skip ghost - // cells in the loop below. + // If we have a distributed::Triangulation only allow locally_owned + // subdomain. Not setting a subdomain is also okay, because we skip + // ghost cells in the loop below. Assert ( (dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id) || @@ -152,20 +147,15 @@ namespace DoFTools const hp::FECollection fe_collection (dof.get_fe()); - // first, for each finite element, build a - // mask for each dof, not like the one - // given which represents components. make - // sure we do the right thing also with - // respect to non-primitive shape - // functions, which takes some additional - // thought + // first, for each finite element, build a mask for each dof, not like + // the one given which represents components. make sure we do the right + // thing also with respect to non-primitive shape functions, which + // takes some additional thought std::vector > dof_mask(fe_collection.size()); - // check whether the table of couplings - // contains only true arguments, i.e., we - // do not exclude any index. that is the - // easy case, since we don't have to set - // up the tables + // check whether the table of couplings contains only true arguments, + // i.e., we do not exclude any index. that is the easy case, since we + // don't have to set up the tables bool need_dof_mask = false; for (unsigned int i=0; iget_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed + // make sparsity pattern for this cell. if no constraints pattern + // was given, then the following call acts as if simply no + // constraints existed constraints.add_entries_local_to_global (dofs_on_this_cell, sparsity, keep_constrained_dofs, @@ -350,10 +337,8 @@ namespace DoFTools { if (DH::dimension == 1) { - // there are only 2 boundary - // indicators in 1d, so it is no - // performance problem to call the - // other function + // there are only 2 boundary indicators in 1d, so it is no + // performance problem to call the other function typename DH::FunctionMap boundary_indicators; boundary_indicators[0] = 0; boundary_indicators[1] = 0; @@ -385,16 +370,11 @@ namespace DoFTools 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) @@ -425,19 +405,16 @@ namespace DoFTools { if (DH::dimension == 1) { - // first check left, then right - // boundary point + // first check left, then right boundary point for (unsigned int direction=0; direction<2; ++direction) { - // if this boundary is not - // requested, then go on with next one + // if this boundary is not requested, then go on with next one if (boundary_indicators.find(direction) == boundary_indicators.end()) continue; - // find active cell at that - // boundary: first go to - // left/right, then to children + // find active cell at that boundary: first go to left/right, + // then to children typename DH::cell_iterator cell = dof.begin(0); while (!cell->at_boundary(direction)) cell = cell->neighbor(direction); @@ -447,8 +424,7 @@ namespace DoFTools const unsigned int dofs_per_vertex = cell->get_fe().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)]; @@ -509,21 +485,23 @@ namespace DoFTools template void - make_flux_sparsity_pattern (const DH &dof, - SparsityPattern &sparsity, - const ConstraintMatrix &constraints, - const bool keep_constrained_dofs, + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity, + const ConstraintMatrix &constraints, + const bool keep_constrained_dofs, const types::subdomain_id subdomain_id) + + // TODO: QA: reduce the indentation level of this method..., Maier 2012 + { const unsigned int n_dofs = dof.n_dofs(); AssertDimension (sparsity.n_rows(), n_dofs); AssertDimension (sparsity.n_cols(), n_dofs); - // If we have a distributed::Triangulation only - // allow locally_owned subdomain. Not setting a - // subdomain is also okay, because we skip ghost - // cells in the loop below. + // If we have a distributed::Triangulation only allow locally_owned + // subdomain. Not setting a subdomain is also okay, because we skip + // ghost cells in the loop below. Assert ( (dof.get_tria().locally_owned_subdomain() == types::invalid_subdomain_id) || @@ -541,18 +519,14 @@ namespace DoFTools typename DH::active_cell_iterator cell = dof.begin_active(), endc = dof.end(); - // TODO: in an old implementation, we used - // user flags before to tag faces that were - // already touched. this way, we could reduce - // the work a little bit. now, we instead add - // only data from one side. this should be OK, - // but we need to actually verify it. - - // In case we work with a distributed - // sparsity pattern of Trilinos type, we - // only have to do the work if the - // current cell is owned by the calling - // processor. Otherwise, just continue. + // TODO: in an old implementation, we used user flags before to tag + // faces that were already touched. this way, we could reduce the work + // a little bit. now, we instead add only data from one side. this + // should be OK, but we need to actually verify it. + + // In case we work with a distributed sparsity pattern of Trilinos + // type, we only have to do the work if the current cell is owned by + // the calling processor. Otherwise, just continue. for (; cell!=endc; ++cell) if (((subdomain_id == types::invalid_subdomain_id) || @@ -564,10 +538,9 @@ namespace DoFTools dofs_on_this_cell.resize (n_dofs_on_this_cell); cell->get_dof_indices (dofs_on_this_cell); - // make sparsity pattern for this - // cell. if no constraints pattern was - // given, then the following call acts - // as if simply no constraints existed + // make sparsity pattern for this cell. if no constraints pattern + // was given, then the following call acts as if simply no + // constraints existed constraints.add_entries_local_to_global (dofs_on_this_cell, sparsity, keep_constrained_dofs); @@ -614,13 +587,12 @@ namespace DoFTools } else { - // Refinement edges are taken care of - // by coarser cells + // Refinement edges are taken care of by coarser + // cells - // TODO: in the distributed case, we miss out - // the constraints when the neighbor cell is - // coarser, but only the current cell is owned - // locally! + // TODO: in the distributed case, we miss out the + // constraints when the neighbor cell is coarser, but + // only the current cell is owned locally! if (cell->neighbor_is_coarser(face)) continue; @@ -634,10 +606,10 @@ namespace DoFTools (dofs_on_this_cell, dofs_on_other_cell, sparsity, keep_constrained_dofs); - // only need to add these in case the neighbor - // cell is not locally owned - otherwise, we - // touch each face twice and hence put the - // indices the other way around + // only need to add these in case the neighbor cell + // is not locally owned - otherwise, we touch each + // face twice and hence put the indices the other way + // around if (cell->neighbor(face)->subdomain_id() != cell->subdomain_id()) constraints.add_entries_local_to_global @@ -650,11 +622,11 @@ namespace DoFTools } + template void - make_flux_sparsity_pattern ( - const DH &dof, - SparsityPattern &sparsity) + make_flux_sparsity_pattern (const DH &dof, + SparsityPattern &sparsity) { ConstraintMatrix constraints; make_flux_sparsity_pattern (dof, sparsity, constraints); @@ -664,9 +636,8 @@ namespace DoFTools template Table<2,Coupling> - dof_couplings_from_component_couplings - (const FiniteElement &fe, - const Table<2,Coupling> &component_couplings) + dof_couplings_from_component_couplings (const FiniteElement &fe, + const Table<2,Coupling> &component_couplings) { Assert(component_couplings.n_rows() == fe.n_components(), ExcDimensionMismatch(component_couplings.n_rows(), @@ -706,6 +677,7 @@ namespace DoFTools } + template std::vector > dof_couplings_from_component_couplings @@ -722,15 +694,13 @@ namespace DoFTools - namespace internal { namespace { - // implementation of the same function in - // namespace DoFTools for non-hp - // DoFHandlers + // implementation of the same function in namespace DoFTools for + // non-hp DoFHandlers template void make_flux_sparsity_pattern (const DH &dof, @@ -801,8 +771,8 @@ namespace DoFTools { typename DH::cell_iterator neighbor = cell->neighbor(face); - // Refinement edges are taken care of - // by coarser cells + // Refinement edges are taken care of by coarser + // cells if (cell->neighbor_is_coarser(face)) continue; @@ -962,9 +932,8 @@ namespace DoFTools } - // implementation of the same function in - // namespace DoFTools for non-hp - // DoFHandlers + // implementation of the same function in namespace DoFTools for + // non-hp DoFHandlers template void make_flux_sparsity_pattern (const dealii::hp::DoFHandler &dof, @@ -972,17 +941,12 @@ namespace DoFTools const Table<2,Coupling> &int_mask, const Table<2,Coupling> &flux_mask) { - // while the implementation above is - // quite optimized and caches a lot of - // data (see e.g. the int/flux_dof_mask - // tables), this is no longer practical - // for the hp version since we would - // have to have it for all combinations - // of elements in the - // hp::FECollection. consequently, the - // implementation here is simpler and - // probably less efficient but at least - // readable... + // while the implementation above is quite optimized and caches a + // lot of data (see e.g. the int/flux_dof_mask tables), this is no + // longer practical for the hp version since we would have to have + // it for all combinations of elements in the hp::FECollection. + // consequently, the implementation here is simpler and probably + // less efficient but at least readable... const dealii::hp::FECollection &fe = dof.get_fe(); @@ -1000,6 +964,7 @@ namespace DoFTools { dofs_on_this_cell.resize (cell->get_fe().dofs_per_cell); cell->get_dof_indices (dofs_on_this_cell); + // make sparsity pattern for this cell for (unsigned int i=0; iget_fe().dofs_per_cell; ++i) for (unsigned int j=0; jget_fe().dofs_per_cell; ++j) @@ -1035,8 +1000,8 @@ namespace DoFTools { typename dealii::hp::DoFHandler::cell_iterator neighbor = cell->neighbor(face); - // Refinement edges are taken care of - // by coarser cells + + // Refinement edges are taken care of by coarser cells if (cell->neighbor_is_coarser(face)) continue; @@ -1167,10 +1132,8 @@ namespace DoFTools const Table<2,Coupling> &int_mask, const Table<2,Coupling> &flux_mask) { - // do the error checking and frame code - // here, and then pass on to more - // specialized functions in the internal - // namespace + // do the error checking and frame code here, and then pass on to more + // specialized functions in the internal namespace const unsigned int n_dofs = dof.n_dofs(); const unsigned int n_comp = dof.get_fe().n_components(); @@ -1187,15 +1150,10 @@ namespace DoFTools Assert (flux_mask.n_cols() == n_comp, ExcDimensionMismatch (flux_mask.n_cols(), n_comp)); - - // Clear user flags because we will - // need them. But first we save - // them and make sure that we - // restore them later such that at - // the end of this function the - // Triangulation will be in the - // same state as it was at the - // beginning of this function. + // Clear user flags because we will need them. But first we save them + // and make sure that we restore them later such that at the end of + // this function the Triangulation will be in the same state as it was + // at the beginning of this function. std::vector user_flags; dof.get_tria().save_user_flags(user_flags); const_cast &>(dof.get_tria()).clear_user_flags (); @@ -1210,15 +1168,11 @@ namespace DoFTools - - - namespace internal { namespace { - inline - bool + inline bool check_master_dof_list (const FullMatrix &face_interpolation_matrix, const std::vector &master_dof_list) { @@ -1229,36 +1183,23 @@ namespace DoFTools for (unsigned int j=0; j p(N); for (unsigned int i=0; i void @@ -1387,42 +1302,27 @@ namespace DoFTools (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; unsigned int index = 0; for (int v=0; @@ -1433,38 +1333,26 @@ namespace DoFTools unsigned int i = 0; while (dofs_added < fe2.dofs_per_vertex) { - // make sure that we - // were able to find - // a set of master - // dofs and that the - // code down below - // didn't just reject - // all our efforts + // make sure that we were able to find a set of master dofs + // and that the code down below didn't just reject all our + // efforts Assert (i < fe1.dofs_per_vertex, ExcInternalError()); - // tentatively push - // this vertex dof + + // tentatively push this vertex dof master_dof_list.push_back (index+i); - // then see what - // happens. if it - // succeeds, fine + // then see what happens. if it succeeds, fine if (check_master_dof_list (face_interpolation_matrix, master_dof_list) == true) ++dofs_added; else - // well, it - // didn't. simply - // pop that dof - // from the list - // again and try - // with the next - // dof + // well, it didn't. simply pop that dof from the list + // again and try with the next dof master_dof_list.pop_back (); - // forward counter by - // one + // forward counter by one ++i; } index += fe1.dofs_per_vertex; @@ -1523,8 +1411,7 @@ namespace DoFTools 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(); i!=master_dof_list.end(); ++i) @@ -1534,10 +1421,8 @@ namespace DoFTools /** - * Make sure that the mask exists that - * determines which dofs will be the - * masters on refined faces where an - * fe1 and a fe2 meet. + * Make sure that the mask exists that determines which dofs will be + * the masters on refined faces where an fe1 and a fe2 meet. */ template void @@ -1558,14 +1443,12 @@ namespace DoFTools } + /** - * Make sure that the given @p - * face_interpolation_matrix - * pointer points to a valid - * matrix. If the pointer is zero - * beforehand, create an entry - * with the correct data. If it - * is nonzero, don't touch it. + * Make sure that the given @p face_interpolation_matrix pointer + * points to a valid matrix. If the pointer is zero beforehand, + * create an entry with the correct data. If it is nonzero, don't + * touch it. */ template void @@ -1586,8 +1469,7 @@ namespace DoFTools /** - * Same, but for subface - * interpolation matrices. + * Same, but for subface interpolation matrices. */ template void @@ -1608,13 +1490,11 @@ namespace DoFTools } + /** - * Given the face interpolation - * matrix between two elements, - * split it into its master and - * slave parts and invert the - * master part as explained in - * the @ref hp_paper "hp paper". + * Given the face interpolation matrix between two elements, split it + * into its master and slave parts and invert the master part as + * explained in the @ref hp_paper "hp paper". */ void ensure_existence_of_split_face_matrix (const FullMatrix &face_interpolation_matrix, @@ -1666,17 +1546,14 @@ namespace DoFTools AssertDimension (nth_master_dof, n_master_dofs); AssertDimension (nth_slave_dof, n_dofs-n_master_dofs); -//TODO[WB]: We should make sure very small entries are removed after inversion + //TODO[WB]: We should make sure very small entries are removed after inversion split_matrix->first.gauss_jordan (); } } - // a template that can - // determine statically whether - // a given DoFHandler class - // supports different finite - // element elements + // a template that can determine statically whether a given + // DoFHandler class supports different finite element elements template struct DoFHandlerSupportsDifferentFEs { @@ -1699,13 +1576,9 @@ namespace DoFTools /** - * A function that returns how - * many different finite - * elements a dof handler - * uses. This is one for non-hp - * DoFHandlers and - * dof_handler.get_fe().size() - * for the hp-versions. + * A function that returns how many different finite elements a dof + * handler uses. This is one for non-hp DoFHandlers and + * dof_handler.get_fe().size() for the hp-versions. */ template unsigned int @@ -1724,13 +1597,9 @@ namespace DoFTools /** - * For a given face belonging - * to an active cell that - * borders to a more refined - * cell, return the fe_index of - * the most dominating finite - * element used on any of the - * face's subfaces. + * For a given face belonging to an active cell that borders to a + * more refined cell, return the fe_index of the most dominating + * finite element used on any of the face's subfaces. */ template unsigned int @@ -1745,12 +1614,9 @@ namespace DoFTools for (; dominating_subface_non_children(); ++dominating_subface_no) { - // each of the subfaces - // can have only a single - // fe_index associated - // with them, since there - // is no cell on the - // other side + // each of the subfaces can have only a single fe_index + // associated with them, since there is no cell on the other + // side Assert (face->child(dominating_subface_no) ->n_active_fe_indices() == 1, @@ -1775,47 +1641,34 @@ namespace DoFTools this_subface_fe.compare_for_face_domination(that_subface_fe); } - // see if the element - // on this subface is - // able to dominate - // the ones on all - // other subfaces, - // and if so take it + // see if the element on this subface is able to dominate the + // ones on all other subfaces, and if so take it if ((domination == FiniteElementDomination::this_element_dominates) || (domination == FiniteElementDomination::either_element_can_dominate)) break; } - // check that we have - // found one such subface + // check that we have found one such subface Assert (dominating_subface_no < face->n_children(), ExcNotImplemented()); - // return the finite element - // index used on it. note - // that only a single fe can - // be active on such subfaces + // return the finite element index used on it. note that only a + // single fe can be active on such subfaces return face->child (dominating_subface_no)->nth_active_fe_index(0); } /** - * 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, @@ -1833,10 +1686,8 @@ namespace DoFTools const unsigned int n_master_dofs = master_dofs.size (); const unsigned int n_slave_dofs = slave_dofs.size (); - // check for a couple - // conditions that happened - // in parallel distributed - // mode + // check for a couple conditions that happened in parallel + // distributed mode for (unsigned int row=0; row!=n_slave_dofs; ++row) Assert (slave_dofs[row] != numbers::invalid_unsigned_int, ExcInternalError()); @@ -1850,10 +1701,8 @@ namespace DoFTools { bool constraint_already_satisfied = false; - // Check if we have an identity - // constraint, which is already - // satisfied by unification of - // the corresponding global dof + // Check if we have an identity constraint, which is already + // satisfied by unification of the corresponding global dof // indices for (unsigned int i=0; i &, ConstraintMatrix &) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1921,8 +1758,7 @@ namespace DoFTools ConstraintMatrix &, dealii::internal::int2type<1>) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1930,8 +1766,7 @@ namespace DoFTools make_hp_hanging_node_constraints (const dealii::MGDoFHandler<1> &, ConstraintMatrix &) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1941,8 +1776,7 @@ namespace DoFTools ConstraintMatrix &, dealii::internal::int2type<1>) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1950,11 +1784,10 @@ namespace DoFTools 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... } @@ -1964,11 +1797,10 @@ namespace DoFTools 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... } @@ -1976,8 +1808,7 @@ namespace DoFTools make_hp_hanging_node_constraints (const dealii::DoFHandler<1,2> &, ConstraintMatrix &) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1987,8 +1818,7 @@ namespace DoFTools ConstraintMatrix &, dealii::internal::int2type<1>) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -1996,8 +1826,7 @@ namespace DoFTools make_hp_hanging_node_constraints (const dealii::DoFHandler<1,3> &, ConstraintMatrix &) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } void @@ -2005,8 +1834,7 @@ namespace DoFTools ConstraintMatrix &, dealii::internal::int2type<1>) { - // nothing to do for regular - // dof handlers in 1d + // nothing to do for regular dof handlers in 1d } @@ -2041,6 +1869,7 @@ namespace DoFTools // // constraints for // // vertices. gotta think about // // that a bit more +// // //TODO[WB]: think about what to do here... // } //#endif @@ -2060,45 +1889,28 @@ namespace DoFTools 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) - // artificial cells can at - // best neighbor ghost cells, - // but we're not interested - // in these interfaces + // artificial cells can at best neighbor ghost cells, but we're not + // interested in these interfaces if (!cell->is_artificial ()) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->has_children()) { - // in any case, faces - // can have at most two - // active fe indices, - // but here the face - // can have only one - // (namely the same as - // that from the cell - // we're sitting on), - // and each of the - // children can have - // only one as - // well. check this + // in any case, faces can have at most two active fe + // indices, but here the face can have only one (namely the + // same as that from the cell we're sitting on), and each + // of the children can have only one as well. check this Assert (cell->face(face)->n_active_fe_indices() == 1, ExcInternalError()); Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) @@ -2109,11 +1921,8 @@ namespace DoFTools Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, ExcInternalError()); - // right now, all that - // is implemented is - // the case that both - // sides use the same - // fe + // right now, all that is implemented is the case that both + // sides use the same fe for (unsigned int c=0; cface(face)->n_children(); ++c) if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) Assert (cell->face(face)->child(c) @@ -2140,8 +1949,7 @@ namespace DoFTools const typename DH::line_iterator this_face = cell->face(face); - // fill the dofs indices. Use same - // enumeration scheme as in + // fill the dofs indices. Use same enumeration scheme as in // @p{FiniteElement::constraints()} unsigned int next_index = 0; for (unsigned int vertex=0; vertex<2; ++vertex) @@ -2162,8 +1970,7 @@ namespace DoFTools = this_face->child(child)->dof_index(dof, fe_index); AssertDimension (next_index, dofs_on_children.size()); - // for each row in the constraint - // matrix for this line: + // for each row in the constraint matrix for this line: for (unsigned int row=0; row!=dofs_on_children.size(); ++row) { constraints.add_line (dofs_on_children[row]); @@ -2177,18 +1984,10 @@ namespace DoFTools } else { - // this face has no - // children, but it - // could still be - // that it is shared - // by two cells that - // use a different fe - // index. check a - // couple of things, - // but ignore the - // case that the - // neighbor is an - // artificial cell + // this face has no children, but it could still be that it + // is shared by two cells that use a different fe index. + // check a couple of things, but ignore the case that the + // neighbor is an artificial cell if (!cell->at_boundary(face) && !cell->neighbor(face)->is_artificial()) { @@ -2214,37 +2013,26 @@ namespace DoFTools 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) - // artificial cells can at - // best neighbor ghost cells, - // but we're not interested - // in these interfaces + // artificial cells can at best neighbor ghost cells, but we're not + // interested in these interfaces if (!cell->is_artificial ()) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->has_children()) { - // first of all, make sure that - // we treat a case which is - // possible, i.e. either no dofs - // on the face at all or no + // first of all, make sure that we treat a case which is + // possible, i.e. either no dofs on the face at all or no // anisotropic refinement if (cell->get_fe().dofs_per_face == 0) continue; @@ -2252,18 +2040,10 @@ namespace DoFTools Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, ExcNotImplemented()); - // in any case, faces - // can have at most two - // active fe indices, - // but here the face - // can have only one - // (namely the same as - // that from the cell - // we're sitting on), - // and each of the - // children can have - // only one as - // well. check this + // in any case, faces can have at most two active fe + // indices, but here the face can have only one (namely the + // same as that from the cell we're sitting on), and each + // of the children can have only one as well. check this AssertDimension (cell->face(face)->n_active_fe_indices(), 1); Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) == true, @@ -2271,15 +2051,9 @@ namespace DoFTools for (unsigned int c=0; cface(face)->n_children(); ++c) AssertDimension (cell->face(face)->child(c)->n_active_fe_indices(), 1); - // right now, all that - // is implemented is - // the case that both - // sides use the same - // fe, and not only - // that but also that - // all lines bounding - // this face and the - // children have the + // right now, all that is implemented is the case that both + // sides use the same fe, and not only that but also that + // all lines bounding this face and the children have the // same fe for (unsigned int c=0; cface(face)->n_children(); ++c) if (!cell->neighbor_child_on_subface(face,c)->is_artificial()) @@ -2315,7 +2089,8 @@ namespace DoFTools const unsigned int n_dofs_on_children = (5*fe.dofs_per_vertex+ 12*fe.dofs_per_line+ 4*fe.dofs_per_quad); -//TODO[TL]: think about this and the following in case of anisotropic refinement + + //TODO[TL]: think about this and the following in case of anisotropic refinement dofs_on_mother.resize (n_dofs_on_mother); dofs_on_children.resize (n_dofs_on_children); @@ -2329,8 +2104,7 @@ namespace DoFTools const typename DH::face_iterator this_face = cell->face(face); - // fill the dofs indices. Use same - // enumeration scheme as in + // fill the dofs indices. Use same enumeration scheme as in // @p{FiniteElement::constraints()} unsigned int next_index = 0; for (unsigned int vertex=0; vertex<4; ++vertex) @@ -2347,9 +2121,11 @@ namespace DoFTools next_index = 0; - // assert some consistency - // assumptions - //TODO[TL]: think about this in case of anisotropic refinement + // assert some consistency assumptions + + //TODO[TL]: think about this in case of anisotropic + //refinement + Assert (dof_handler.get_tria().get_anisotropic_refinement_flag() || ((this_face->child(0)->vertex_index(3) == this_face->child(1)->vertex_index(2)) && @@ -2362,17 +2138,16 @@ namespace DoFTools dofs_on_children[next_index++] = this_face->child(0)->vertex_dof_index(3,dof); - // dof numbers on the centers of - // the lines bounding this face + // dof numbers on the centers of the lines bounding this + // face for (unsigned int line=0; line<4; ++line) for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) dofs_on_children[next_index++] = this_face->line(line)->child(0)->vertex_dof_index(1,dof, fe_index); - // next the dofs on the lines interior - // to the face; the order of these - // lines is laid down in the - // FiniteElement class documentation + // next the dofs on the lines interior to the face; the + // order of these lines is laid down in the FiniteElement + // class documentation for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof, fe_index); @@ -2393,16 +2168,14 @@ namespace DoFTools dofs_on_children[next_index++] = this_face->line(line)->child(child)->dof_index(dof, fe_index); - // finally, for the dofs interior - // to the four child faces + // finally, for the dofs interior to the four child faces for (unsigned int child=0; child<4; ++child) for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) dofs_on_children[next_index++] = this_face->child(child)->dof_index(dof, fe_index); AssertDimension (next_index, dofs_on_children.size()); - // for each row in the constraint - // matrix for this line: + // for each row in the constraint matrix for this line: for (unsigned int row=0; row!=dofs_on_children.size(); ++row) { constraints.add_line (dofs_on_children[row]); @@ -2416,18 +2189,10 @@ namespace DoFTools } else { - // this face has no - // children, but it - // could still be - // that it is shared - // by two cells that - // use a different fe - // index. check a - // couple of things, - // but ignore the - // case that the - // neighbor is an - // artificial cell + // this face has no children, but it could still be that it + // is shared by two cells that use a different fe index. + // check a couple of things, but ignore the case that the + // neighbor is an artificial cell if (!cell->at_boundary(face) && !cell->neighbor(face)->is_artificial()) { @@ -2446,43 +2211,31 @@ namespace DoFTools make_hp_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints) { - // note: this function is going - // to be hard to understand if - // you haven't read the hp - // paper. however, we try to - // follow the notation laid out - // there, so go read the paper - // before you try to understand + // note: this function is going to be hard to understand if you + // haven't read the hp paper. however, we try to follow the notation + // laid out there, so go read the paper before you try to understand // what is going on here + const unsigned int dim = DH::dimension; 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; - // 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)); @@ -2491,50 +2244,37 @@ namespace DoFTools n_finite_elements (dof_handler), GeometryInfo::max_children_per_face); - // similarly have a cache for - // the matrices that are split - // into their master and slave - // parts, and for which the - // master part is - // inverted. these two matrices - // are derived from the face - // interpolation matrix as - // described in the @ref hp_paper "hp paper" + // similarly have a cache for the matrices that are split into their + // master and slave parts, and for which the master part is inverted. + // these two matrices are derived from the face interpolation matrix + // as described in the @ref hp_paper "hp paper" Table<2,std_cxx1x::shared_ptr,FullMatrix > > > split_face_interpolation_matrices (n_finite_elements (dof_handler), n_finite_elements (dof_handler)); - // finally, for each pair of finite - // elements, have a mask that states - // which of the degrees of freedom on - // the coarse side of a refined face - // will act as master dofs. + // finally, for each pair of finite elements, have a mask that states + // which of the degrees of freedom on the coarse side of a refined + // face will act as master dofs. Table<2,std_cxx1x::shared_ptr > > master_dof_masks (n_finite_elements (dof_handler), n_finite_elements (dof_handler)); // loop over all faces // - // 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) - // artificial cells can at - // best neighbor ghost cells, - // but we're not interested - // in these interfaces + // artificial cells can at best neighbor ghost cells, but we're not + // interested in these interfaces if (!cell->is_artificial ()) for (unsigned int face=0; face::faces_per_cell; ++face) if (cell->face(face)->has_children()) { - // first of all, make sure that - // we treat a case which is - // possible, i.e. either no dofs - // on the face at all or no + // first of all, make sure that we treat a case which is + // possible, i.e. either no dofs on the face at all or no // anisotropic refinement if (cell->get_fe().dofs_per_face == 0) continue; @@ -2542,25 +2282,14 @@ namespace DoFTools Assert(cell->face(face)->refinement_case()==RefinementCase::isotropic_refinement, ExcNotImplemented()); - // so now we've found a - // face of an active - // cell that has - // children. that means - // that there are - // hanging nodes here. - - // in any case, faces - // can have at most two - // sets of active fe - // indices, but here - // the face can have - // only one (namely the - // same as that from - // the cell we're - // sitting on), and - // each of the children - // can have only one as - // well. check this + // so now we've found a face of an active cell that has + // children. that means that there are hanging nodes here. + + // in any case, faces can have at most two sets of active + // fe indices, but here the face can have only one (namely + // the same as that from the cell we're sitting on), and + // each of the children can have only one as well. check + // this Assert (cell->face(face)->n_active_fe_indices() == 1, ExcInternalError()); Assert (cell->face(face)->fe_index_is_active(cell->active_fe_index()) @@ -2570,23 +2299,13 @@ namespace DoFTools Assert (cell->face(face)->child(c)->n_active_fe_indices() == 1, ExcInternalError()); - // first find out - // whether we can - // constrain each of - // the subfaces to the - // mother face. in the - // lingo of the hp - // paper, this would be - // the simple - // case. note that we - // can short-circuit - // this decision if the - // dof_handler doesn't + // first find out whether we can constrain each of the + // subfaces to the mother face. in the lingo of the hp + // paper, this would be the simple case. note that we can + // short-circuit this decision if the dof_handler doesn't // support hp at all // - // ignore all - // interfaces with - // artificial cells + // ignore all interfaces with artificial cells FiniteElementDomination::Domination mother_face_dominates = FiniteElementDomination::either_element_can_dominate; @@ -2602,71 +2321,31 @@ namespace DoFTools case FiniteElementDomination::this_element_dominates: case FiniteElementDomination::either_element_can_dominate: { - // Case 1 (the - // simple case - // and the only - // case that can - // happen for - // non-hp - // DoFHandlers): - // The coarse - // element - // dominates the - // elements on - // the subfaces - // (or they are + // Case 1 (the simple case and the only case that can + // happen for non-hp DoFHandlers): The coarse element + // dominates the elements on the subfaces (or they are // all the same) // - // so we are - // going to - // constrain - // the DoFs on - // the face - // children - // against the - // DoFs on the - // face itself + // so we are going to constrain the DoFs on the face + // children against the DoFs on the face itself master_dofs.resize (cell->get_fe().dofs_per_face); cell->face(face)->get_dof_indices (master_dofs, cell->active_fe_index ()); - // Now create - // constraint matrix - // for the subfaces and - // assemble it. ignore - // all interfaces with - // artificial cells - // because we can only - // get to such - // interfaces if the - // current cell is a - // ghost cell. also - // ignore the interface - // if the neighboring - // cell is a ghost cell - // in 2d: what we would - // compute here are the - // constraints on the - // ghost cell's DoFs, - // but we are not - // interested in those: - // we only want - // constraints on - // *locally active* - // DoFs, not on - // *locally relevant* - // DoFs. However, in 3d - // we must still - // compute those - // constraints because - // it might happen that - // a constraint is - // related to an edge - // where the hanging - // node is only - // detected if we also - // look between ghosts + // Now create constraint matrix for the subfaces and + // assemble it. ignore all interfaces with artificial + // cells because we can only get to such interfaces if + // the current cell is a ghost cell. also ignore the + // interface if the neighboring cell is a ghost cell in + // 2d: what we would compute here are the constraints + // on the ghost cell's DoFs, but we are not interested + // in those: we only want constraints on *locally + // active* DoFs, not on *locally relevant* DoFs. + // However, in 3d we must still compute those + // constraints because it might happen that a + // constraint is related to an edge where the hanging + // node is only detected if we also look between ghosts for (unsigned int c=0; cface(face)->n_children(); ++c) { if (cell->neighbor_child_on_subface (face, c)->is_artificial() @@ -2683,43 +2362,23 @@ namespace DoFTools const unsigned int subface_fe_index = subface->nth_active_fe_index(0); - // we sometime run - // into the - // situation where - // for example on - // one big cell we - // have a FE_Q(1) - // and on the - // subfaces we have - // a mixture of - // FE_Q(1) and - // FE_Nothing. In - // that case, the - // face domination - // is - // either_element_can_dominate - // for the whole - // collection of - // subfaces, but on - // the particular - // subface between - // FE_Q(1) and - // FE_Nothing, - // there are no - // constraints that - // we need to take - // care of. in that - // case, just - // continue + // we sometime run into the situation where for + // example on one big cell we have a FE_Q(1) and on + // the subfaces we have a mixture of FE_Q(1) and + // FE_Nothing. In that case, the face domination is + // either_element_can_dominate for the whole + // collection of subfaces, but on the particular + // subface between FE_Q(1) and FE_Nothing, there + // are no constraints that we need to take care of. + // in that case, just continue if (cell->get_fe().compare_for_face_domination (subface->get_fe(subface_fe_index)) == FiniteElementDomination::no_requirements) continue; - // Same procedure as for the - // mother cell. Extract the face - // DoFs from the cell DoFs. + // Same procedure as for the mother cell. Extract + // the face DoFs from the cell DoFs. slave_dofs.resize (subface->get_fe(subface_fe_index) .dofs_per_face); subface->get_dof_indices (slave_dofs, subface_fe_index); @@ -2728,60 +2387,28 @@ namespace DoFTools Assert (slave_dofs[i] != numbers::invalid_unsigned_int, ExcInternalError()); - // Now create the - // element constraint - // for this subface. + // Now create the element constraint for this + // subface. // - // As a side remark, - // one may wonder the - // following: - // neighbor_child is - // clearly computed - // correctly, - // i.e. taking into - // account - // face_orientation - // (just look at the - // implementation of - // that - // function). however, - // we don't care about - // this here, when we - // ask for - // subface_interpolation - // on subface c. the - // question rather is: - // do we have to - // translate 'c' here - // as well? + // As a side remark, one may wonder the following: + // neighbor_child is clearly computed correctly, + // i.e. taking into account face_orientation (just + // look at the implementation of that function). + // however, we don't care about this here, when we + // ask for subface_interpolation on subface c. the + // question rather is: do we have to translate 'c' + // here as well? // - // the answer is in - // fact 'no'. if one - // does that, results - // are wrong: - // constraints are - // added twice for the - // same pair of nodes - // but with differing - // weights. in - // addition, one can - // look at the - // deal.II/project_*_03 - // tests that look at - // exactly this case: - // there, we have a - // mesh with at least - // one - // face_orientation==false - // and hanging nodes, - // and the results of - // those tests show - // that the result of - // projection verifies - // the approximation - // properties of a - // finite element onto - // that mesh + // the answer is in fact 'no'. if one does that, + // results are wrong: constraints are added twice + // for the same pair of nodes but with differing + // weights. in addition, one can look at the + // deal.II/project_*_03 tests that look at exactly + // this case: there, we have a mesh with at least + // one face_orientation==false and hanging nodes, + // and the results of those tests show that the + // result of projection verifies the approximation + // properties of a finite element onto that mesh ensure_existence_of_subface_matrix (cell->get_fe(), subface->get_fe(subface_fe_index), @@ -2789,8 +2416,7 @@ namespace DoFTools subface_interpolation_matrices [cell->active_fe_index()][subface_fe_index][c]); - // Add constraints to global constraint - // matrix. + // Add constraints to global constraint matrix. filter_constraints (master_dofs, slave_dofs, *(subface_interpolation_matrices @@ -2804,54 +2430,28 @@ namespace DoFTools case FiniteElementDomination::other_element_dominates: case FiniteElementDomination::neither_element_dominates: { - // Case 2 (the "complex" - // case): at least one - // (the neither_... case) - // of the finer elements - // or all of them (the - // other_... case) is - // dominating. See the hp - // paper for a way how to - // deal with this - // situation + // Case 2 (the "complex" case): at least one (the + // neither_... case) of the finer elements or all of + // them (the other_... case) is dominating. See the hp + // paper for a way how to deal with this situation // - // since this is - // something that - // can only - // happen for hp - // dof handlers, - // add a check - // here... + // since this is something that can only happen for hp + // dof handlers, add a check here... Assert (DoFHandlerSupportsDifferentFEs::value == true, ExcInternalError()); - // we first have - // to find the - // finite element - // that is able - // to generate a - // space that all - // the other ones - // can be - // constrained to + // we first have to find the finite element that is + // able to generate a space that all the other ones can + // be constrained to const unsigned int dominating_fe_index = get_most_dominating_subface_fe_index (cell->face(face)); const FiniteElement &dominating_fe = dof_handler.get_fe()[dominating_fe_index]; - // check also - // that it is - // able to - // constrain the - // mother - // face. it - // should be, or - // we wouldn't - // have gotten - // into the - // branch for the - // 'complex' case + // check also that it is able to constrain the mother + // face. it should be, or we wouldn't have gotten into + // the branch for the 'complex' case Assert ((dominating_fe.compare_for_face_domination (cell->face(face)->get_fe(cell->face(face)->nth_active_fe_index(0))) == FiniteElementDomination::this_element_dominates) @@ -2862,10 +2462,8 @@ namespace DoFTools ExcInternalError()); - // first get the - // interpolation matrix - // from the mother to the - // virtual dofs + // first get the interpolation matrix from the mother + // to the virtual dofs Assert (dominating_fe.dofs_per_face <= cell->get_fe().dofs_per_face, ExcInternalError()); @@ -2876,10 +2474,8 @@ namespace DoFTools face_interpolation_matrices [dominating_fe_index][cell->active_fe_index()]); - // split this matrix into - // master and slave - // components. invert the - // master component + // split this matrix into master and slave components. + // invert the master component ensure_existence_of_master_dof_mask (cell->get_fe(), dominating_fe, @@ -2906,14 +2502,8 @@ namespace DoFTools = (split_face_interpolation_matrices [dominating_fe_index][cell->active_fe_index()]->second); - // now compute - // the constraint - // matrix as the - // product - // between the - // inverse matrix - // and the slave - // part + // now compute the constraint matrix as the product + // between the inverse matrix and the slave part constraint_matrix.reinit (cell->get_fe().dofs_per_face - dominating_fe.dofs_per_face, dominating_fe.dofs_per_face); @@ -2921,19 +2511,13 @@ namespace DoFTools .mmult (constraint_matrix, restrict_mother_to_virtual_master_inv); - // then figure - // out the global - // numbers of - // master and - // slave dofs and - // apply - // constraints + // then figure out the global numbers of master and + // slave dofs and apply constraints scratch_dofs.resize (cell->get_fe().dofs_per_face); cell->face(face)->get_dof_indices (scratch_dofs, cell->active_fe_index ()); - // split dofs into master - // and slave components + // split dofs into master and slave components master_dofs.clear (); slave_dofs.clear (); for (unsigned int i=0; iget_fe().dofs_per_face; ++i) @@ -2954,26 +2538,13 @@ namespace DoFTools - // next we have - // to deal with - // the - // subfaces. do - // as discussed - // in the hp - // paper + // next we have to deal with the subfaces. do as + // discussed in the hp paper for (unsigned int sf=0; sfface(face)->n_children(); ++sf) { - // ignore - // interfaces - // with - // artificial - // cells as - // well as - // interfaces - // between - // ghost - // cells in 2d + // ignore interfaces with artificial cells as well + // as interfaces between ghost cells in 2d if (cell->neighbor_child_on_subface (face, sf)->is_artificial() || (dim==2 && cell->is_ghost() @@ -2990,11 +2561,8 @@ namespace DoFTools const FiniteElement &subface_fe = dof_handler.get_fe()[subface_fe_index]; - // first get the - // interpolation - // matrix from the - // subface to the - // virtual dofs + // first get the interpolation matrix from the + // subface to the virtual dofs Assert (dominating_fe.dofs_per_face <= subface_fe.dofs_per_face, ExcInternalError()); @@ -3030,15 +2598,8 @@ namespace DoFTools } case FiniteElementDomination::no_requirements: - // there - // are no - // continuity - // requirements - // between - // the two - // elements. record - // no - // constraints + // there are no continuity requirements between the two + // elements. record no constraints break; default: @@ -3048,40 +2609,24 @@ namespace DoFTools } else { - // this face has no - // children, but it - // could still be that - // it is shared by two - // cells that use a - // different fe index + // this face has no children, but it could still be that it + // is shared by two cells that use a different fe index Assert (cell->face(face) ->fe_index_is_active(cell->active_fe_index()) == true, ExcInternalError()); - // see if there is a - // neighbor that is - // an artificial - // cell. in that - // case, we're not - // interested in this - // interface. we test - // this case first - // since artificial - // cells may not have - // an active_fe_index - // set, etc + // see if there is a neighbor that is an artificial cell. + // in that case, we're not interested in this interface. we + // test this case first since artificial cells may not have + // an active_fe_index set, etc if (!cell->at_boundary(face) && cell->neighbor(face)->is_artificial()) continue; - // Only if there is - // a neighbor with - // a different - // active_fe_index - // and the same h-level, - // some action has - // to be taken. + // Only if there is a neighbor with a different + // active_fe_index and the same h-level, some action has to + // be taken. if ((DoFHandlerSupportsDifferentFEs::value == true) && !cell->face(face)->at_boundary () @@ -3094,17 +2639,13 @@ namespace DoFTools { const typename DH::cell_iterator neighbor = cell->neighbor (face); - // see which side of the - // face we have to - // constrain + // see which side of the face we have to constrain switch (cell->get_fe().compare_for_face_domination (neighbor->get_fe ())) { case FiniteElementDomination::this_element_dominates: { - // Get DoFs on - // dominating and - // dominated side of - // the face + // Get DoFs on dominating and dominated side of the + // face master_dofs.resize (cell->get_fe().dofs_per_face); cell->face(face)->get_dof_indices (master_dofs, cell->active_fe_index ()); @@ -3113,26 +2654,20 @@ namespace DoFTools cell->face(face)->get_dof_indices (slave_dofs, neighbor->active_fe_index ()); - // break if the n_master_dofs == 0, - // because we are attempting to - // constrain to an element that has + // break if the n_master_dofs == 0, because we are + // attempting to constrain to an element that has // has no face dofs if (master_dofs.size() == 0) break; - // make sure - // the element - // constraints - // for this - // face are - // available + // make sure the element constraints for this face + // are available ensure_existence_of_face_matrix (cell->get_fe(), neighbor->get_fe(), face_interpolation_matrices [cell->active_fe_index()][neighbor->active_fe_index()]); - // Add constraints to global constraint - // matrix. + // Add constraints to global constraint matrix. filter_constraints (master_dofs, slave_dofs, *(face_interpolation_matrices @@ -3145,61 +2680,34 @@ namespace DoFTools case FiniteElementDomination::other_element_dominates: { - // we don't do anything - // here since we will - // come back to this - // face from the other - // cell, at which time - // we will fall into - // the first case - // clause above + // we don't do anything here since we will come + // back to this face from the other cell, at which + // time we will fall into the first case clause + // above break; } case FiniteElementDomination::either_element_can_dominate: { - // it appears as if - // neither element has - // any constraints on - // its neighbor. this - // may be because - // neither element has - // any DoFs on faces at - // all. or that the two - // elements are - // actually the same, - // although they happen - // to run under - // different fe_indices - // (this is what - // happens in - // hp/hp_hanging_nodes_01 - // for example). + // it appears as if neither element has any + // constraints on its neighbor. this may be because + // neither element has any DoFs on faces at all. or + // that the two elements are actually the same, + // although they happen to run under different + // fe_indices (this is what happens in + // hp/hp_hanging_nodes_01 for example). // - // another possibility - // is what happens in - // crash_13. there, we - // have - // FESystem(FE_Q(1),FE_DGQ(0)) - // vs. FESystem(FE_Q(1),FE_DGQ(1)). - // neither of them - // dominates the - // other. the point is - // that it doesn't - // matter, since - // hopefully for this - // case, both sides' - // dofs should have - // been unified. + // another possibility is what happens in crash_13. + // there, we have FESystem(FE_Q(1),FE_DGQ(0)) vs. + // FESystem(FE_Q(1),FE_DGQ(1)). neither of them + // dominates the other. the point is that it + // doesn't matter, since hopefully for this case, + // both sides' dofs should have been unified. // - // make sure this is - // actually true. this - // actually only - // matters, of course, - // if either of the two - // finite elements - // actually do have - // dofs on the face + // make sure this is actually true. this actually + // only matters, of course, if either of the two + // finite elements actually do have dofs on the + // face if ((cell->get_fe().dofs_per_face != 0) || (cell->neighbor(face)->get_fe().dofs_per_face != 0)) @@ -3209,10 +2717,8 @@ namespace DoFTools cell->neighbor(face)->get_fe().dofs_per_face, ExcNotImplemented()); - // (ab)use the master - // and slave dofs - // arrays for a - // moment here + // (ab)use the master and slave dofs arrays for + // a moment here master_dofs.resize (cell->get_fe().dofs_per_face); cell->face(face) ->get_dof_indices (master_dofs, @@ -3232,16 +2738,11 @@ namespace DoFTools case FiniteElementDomination::neither_element_dominates: { - // we don't presently - // know what exactly to - // do here. it isn't quite - // clear what exactly we - // would have to do - // here. sit tight until - // someone trips over the - // following statement - // and see what exactly - // is going on + // we don't presently know what exactly to do here. + // it isn't quite clear what exactly we would have + // to do here. sit tight until someone trips over + // the following statement and see what exactly is + // going on Assert (false, ExcNotImplemented()); break; } @@ -3253,8 +2754,7 @@ namespace DoFTools } default: - // we shouldn't get - // here + // we shouldn't get here Assert (false, ExcInternalError()); } } @@ -3270,14 +2770,10 @@ namespace DoFTools make_hanging_node_constraints (const DH &dof_handler, ConstraintMatrix &constraints) { - // Decide whether to use the - // new or old make_hanging_node_constraints - // function. If all the FiniteElement - // or all elements in a FECollection support - // the new face constraint matrix, the - // new code will be used. - // Otherwise, the old implementation is used - // for the moment. + // Decide whether to use the new or old make_hanging_node_constraints + // function. If all the FiniteElement or all elements in a FECollection + // support the new face constraint matrix, the new code will be used. + // Otherwise, the old implementation is used for the moment. if (dof_handler.get_fe().hp_constraints_are_implemented ()) internal:: make_hp_hanging_node_constraints (dof_handler, @@ -3329,28 +2825,30 @@ namespace DoFTools // A lookup table on how to go through the child faces depending on the // orientation: - static const int lookup_table_2d[2][2] - // flip - = { {0, 1}, // false - {1, 0}, // true + static const int lookup_table_2d[2][2] = + { + // flip: + {0, 1}, // false + {1, 0}, // true }; - static const int lookup_table_3d[2][2][2][4] - // orientation flip rotation - = { { { {0, 2, 1, 3}, // false false false - {2, 3, 0, 1}, - }, // false false true - { {3, 1, 2, 0}, // false true false - {1, 0, 3, 2}, + static const int lookup_table_3d[2][2][2][4] = + { + // orientation flip rotation + { { {0, 2, 1, 3}, // false false false + {2, 3, 0, 1}, // false false true + }, + { {3, 1, 2, 0}, // false true false + {1, 0, 3, 2}, // false true true + }, + }, + { { {0, 1, 2, 3}, // true false false + {1, 3, 0, 2}, // true false true }, - }, // false true true - { { {0, 1, 2, 3}, // true false false - {1, 3, 0, 2}, - }, // true false true - { {3, 2, 1, 0}, // true true false - {2, 0, 3, 1}, + { {3, 2, 1, 0}, // true true false + {2, 0, 3, 1}, // true true true }, - }, // true true true + }, }; // In the case that both faces have children, we loop over all @@ -3431,19 +2929,20 @@ namespace DoFTools // face_flip, // face_rotation) // But this will give us an index on a cell - something we cannot work - // with. But luckily we can match them back :-] + // with directly. But luckily we can match them back :-] std::map cell_to_rotated_face_index; // Build up a cell to face index for face_2: for (unsigned int i = 0; i < dofs_per_face; ++i) { - const unsigned int cell_index = fe. - face_to_cell_index(i, - 0, // It doesn't really matter, just assume - // we're on the first face... - true, false, false // default orientation - ); + const unsigned int cell_index = + fe.face_to_cell_index(i, + 0, /* It doesn't really matter, just assume + * we're on the first face... + */ + true, false, false // default orientation + ); cell_to_rotated_face_index[cell_index] = i; } @@ -3452,18 +2951,17 @@ namespace DoFTools // Query the correct face_index on face_2 respecting the given // orientation: const unsigned int j = - cell_to_rotated_face_index[fe. - face_to_cell_index(i, - 0, // It doesn't really matter, just assume - // we're on the first face... - face_orientation, - face_flip, - face_rotation)]; + cell_to_rotated_face_index[fe.face_to_cell_index(i, + 0, /* It doesn't really matter, just assume + * we're on the first face... + */ + face_orientation, + face_flip, + face_rotation)]; // And finally constrain the two DoFs respecting component_mask: - if ((component_mask.n_selected_components(fe.n_components()) == fe.n_components()) - || - (component_mask[fe.face_system_to_component_index(i).first] == true)) + if (component_mask.n_selected_components(fe.n_components()) == fe.n_components() + || component_mask[fe.face_system_to_component_index(i).first]) { if (!constraint_matrix.is_constrained(dofs_1[i])) { @@ -3511,9 +3009,7 @@ namespace DoFTools Assert (0<=direction && direction*> - (&dof_handler.get_tria()) - == - 0), + (&dof_handler.get_tria()) == 0), ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); @@ -3600,9 +3096,7 @@ namespace DoFTools ExcNotImplemented()); Assert ((dynamic_cast*> - (&dof_handler.get_tria()) - == - 0), + (&dof_handler.get_tria()) == 0), ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); @@ -3714,20 +3208,16 @@ namespace DoFTools } - // this internal function assigns to each dof - // the respective component of the vector - // system. + // this internal function assigns to each dof the respective component + // of the vector system. // - // the output array dofs_by_component - // lists for each dof the corresponding - // vector component. if the DoFHandler is - // based on a parallel distributed - // triangulation then the output array is - // index by + // the output array dofs_by_component lists for each dof the + // corresponding vector component. if the DoFHandler is based on a + // parallel distributed triangulation then the output array is index by // dof.locally_owned_dofs().index_within_set(indices[i]) // - // if an element is non-primitive then we assign to each degree of freedom - // the following component: + // if an element is non-primitive then we assign to each degree of + // freedom the following component: // - if the nonzero components that belong to a shape function are not // selected in the component_mask, then the shape function is assigned // to the first nonzero vector component that corresponds to this @@ -3747,18 +3237,13 @@ namespace DoFTools ExcDimensionMismatch(dofs_by_component.size(), dof.n_locally_owned_dofs())); - // next set up a table for the degrees - // of freedom on each of the cells - // (regardless of the fact whether it - // is listed in the component_select - // argument or not) + // next set up a table for the degrees of freedom on each of the + // cells (regardless of the fact whether it is listed in the + // component_select argument or not) // - // for each element 'f' of the - // FECollection, - // local_component_association[f][d] - // then returns the vector component - // that degree of freedom 'd' belongs - // to + // for each element 'f' of the FECollection, + // local_component_association[f][d] then returns the vector + // component that degree of freedom 'd' belongs to std::vector > local_component_association (fe_collection.size()); for (unsigned int f=0; f indices; for (typename DH::active_cell_iterator c=dof.begin_active(); c!=dof.end(); ++ c) @@ -3788,15 +3272,12 @@ namespace DoFTools } - // this is the function corresponding to the one above but working - // on blocks instead of components. + // this is the function corresponding to the one above but working on + // blocks instead of components. // - // the output array dofs_by_block - // lists for each dof the corresponding - // vector block. if the DoFHandler is - // based on a parallel distributed - // triangulation then the output array is - // index by + // the output array dofs_by_block lists for each dof the corresponding + // vector block. if the DoFHandler is based on a parallel distributed + // triangulation then the output array is index by // dof.locally_owned_dofs().index_within_set(indices[i]) template inline @@ -3811,18 +3292,13 @@ namespace DoFTools ExcDimensionMismatch(dofs_by_block.size(), dof.n_locally_owned_dofs())); - // next set up a table for the degrees - // of freedom on each of the cells - // (regardless of the fact whether it - // is listed in the component_select - // argument or not) + // next set up a table for the degrees of freedom on each of the + // cells (regardless of the fact whether it is listed in the + // component_select argument or not) // - // for each element 'f' of the - // FECollection, - // local_block_association[f][d] - // then returns the vector block - // that degree of freedom 'd' belongs - // to + // for each element 'f' of the FECollection, + // local_block_association[f][d] then returns the vector block that + // degree of freedom 'd' belongs to std::vector > local_block_association (fe_collection.size()); for (unsigned int f=0; f indices; for (typename DH::active_cell_iterator c=dof.begin_active(); c!=dof.end(); ++ c) @@ -3880,15 +3355,12 @@ namespace DoFTools Assert (fe_is_primitive(dof_handler) == true, typename FiniteElement::ExcFENotPrimitive()); - // store a flag whether we should care - // about different components. this is - // just a simplification, we could ask - // for this at every single place - // equally well + // store a flag whether we should care about different components. this + // is just a simplification, we could ask for this at every single + // place equally well const bool consider_components = (n_components(dof_handler) != 1); - // zero out the components that we - // will touch + // zero out the components that we will touch if (consider_components == false) dof_data = 0; else @@ -3904,8 +3376,7 @@ namespace DoFTools dof_data(i) = 0; } - // count how often we have added a value - // in the sum for each dof + // count how often we have added a value in the sum for each dof std::vector touch_count (dof_handler.n_dofs(), 0); typename DH::active_cell_iterator cell = dof_handler.begin_active(), @@ -3920,31 +3391,25 @@ namespace DoFTools cell->get_dof_indices (dof_indices); for (unsigned int i=0; iget_fe().system_to_component_index(i).first == component)) { - // sum up contribution of the - // present_cell to this dof + // sum up contribution of the present_cell to this dof dof_data(dof_indices[i]) += cell_data(present_cell); - // note that we added another - // summand + + // note that we added another summand ++touch_count[dof_indices[i]]; } } - // compute the mean value on all the - // dofs by dividing with the number + // compute the mean value on all the dofs by dividing with the number // of summands. for (unsigned int i=0; i dofs_by_component (dof.n_locally_owned_dofs()); internal::get_component_association (dof, component_mask, dofs_by_component); @@ -4000,7 +3462,7 @@ namespace DoFTools } -//TODO: Unify the following two functions with the non-hp case + // TODO: Unify the following two functions with the non-hp case template void @@ -4016,10 +3478,8 @@ namespace DoFTools Assert(selected_dofs.size() == dof.n_dofs(), ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs())); - // two special cases: no component - // is selected, and all components - // are selected; both rather - // stupid, but easy to catch + // two special cases: no component is selected, and all components are + // selected; both rather stupid, but easy to catch if (component_mask.n_selected_components(n_components(dof)) == 0) { std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); @@ -4035,9 +3495,8 @@ namespace DoFTools // preset all values by false std::fill_n (selected_dofs.begin(), dof.n_dofs(), false); - // get the component association of each DoF - // and then select the ones that match the - // given set of components + // get the component association of each DoF and then select the ones + // that match the given set of components std::vector dofs_by_component (dof.n_dofs()); internal::get_component_association (dof, component_mask, dofs_by_component); @@ -4090,10 +3549,8 @@ namespace DoFTools Assert(selected_dofs.size() == dof.n_dofs(level), ExcDimensionMismatch(selected_dofs.size(), dof.n_dofs(level))); - // two special cases: no component - // is selected, and all components - // are selected, both rather - // stupid, but easy to catch + // two special cases: no component is selected, and all components are + // selected, both rather stupid, but easy to catch if (component_mask.n_selected_components(n_components(dof)) == 0) { std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false); @@ -4108,18 +3565,15 @@ namespace DoFTools // preset all values by false std::fill_n (selected_dofs.begin(), dof.n_dofs(level), false); - // next set up a table for the - // degrees of freedom on each of - // the cells whether it is - // something interesting or not + // next set up a table for the degrees of freedom on each of the cells + // whether it is something interesting or not std::vector local_component_asssociation = internal::get_local_component_association (fe, component_mask); std::vector local_selected_dofs (fe.dofs_per_cell); for (unsigned int i=0; i indices(fe.dofs_per_cell); typename MGDoFHandler::cell_iterator c; for (c = dof.begin(level) ; c != dof.end(level) ; ++ c) @@ -4190,15 +3644,12 @@ namespace DoFTools selected_dofs.clear (); selected_dofs.set_size(dof_handler.n_dofs()); - // let's see whether we have to - // check for certain boundary - // indicators or whether we can - // accept all + // let's see whether we have to check for certain boundary indicators + // or whether we can accept all const bool check_boundary_indicator = (boundary_indicators.size() != 0); - // also see whether we have to - // check whether a certain vector - // component is selected, or all + // also see whether we have to check whether a certain vector component + // is selected, or all const bool check_vector_component = ((component_mask.represents_the_all_selected_mask() == false) || @@ -4208,23 +3659,17 @@ namespace DoFTools 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) @@ -4244,23 +3689,12 @@ namespace DoFTools if (!check_vector_component) selected_dofs.add_index (face_dof_indices[i]); else - // check for - // component is - // required. somewhat - // tricky as usual - // for the case that - // the shape function - // is non-primitive, - // but use usual - // convention (see - // docs) + // check for component is required. somewhat tricky as + // usual for the case that the shape function is + // non-primitive, but use usual convention (see docs) { - // first get at the - // cell-global - // number of a face - // dof, to ask the - // fe certain - // questions + // first get at the cell-global number of a face dof, + // to ask the fe certain questions const unsigned int cell_index = (dim == 1 ? i @@ -4312,38 +3746,27 @@ namespace DoFTools Assert (boundary_indicators.find (numbers::internal_face_boundary_id) == boundary_indicators.end(), ExcInvalidBoundaryIndicator()); - // let's see whether we have to - // check for certain boundary - // indicators or whether we can - // accept all + // let's see whether we have to check for certain boundary indicators + // or whether we can accept all const bool check_boundary_indicator = (boundary_indicators.size() != 0); - // also see whether we have to - // check whether a certain vector - // component is selected, or all + // also see whether we have to check whether a certain vector component + // is selected, or all 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; 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; @@ -4365,16 +3788,9 @@ namespace DoFTools if (!check_vector_component) selected_dofs[cell_dof_indices[i]] = true; else - // check for - // component is - // required. somewhat - // tricky as usual - // for the case that - // the shape function - // is non-primitive, - // but use usual - // convention (see - // docs) + // check for component is required. somewhat tricky + // as usual for the case that the shape function is + // non-primitive, but use usual convention (see docs) { if (fe.is_primitive (i)) selected_dofs[cell_dof_indices[i]] @@ -4429,8 +3845,7 @@ namespace DoFTools const FiniteElement &fe = dof_handler.get_fe(); - // this function is similar to the - // make_sparsity_pattern function, + // this function is similar to the make_sparsity_pattern function, // see there for more information typename dealii::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), @@ -4465,8 +3880,7 @@ namespace DoFTools const FiniteElement &fe = dof_handler.get_fe(); - // this function is similar to the - // make_sparsity_pattern function, + // this function is similar to the make_sparsity_pattern function, // see there for more information typename dealii::DoFHandler::active_cell_iterator @@ -4482,16 +3896,15 @@ namespace DoFTools for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) selected_dofs[face->child(0)->vertex_dof_index(2,dof)] = true; - // dof numbers on the centers of - // the lines bounding this face + // dof numbers on the centers of the lines bounding this + // face for (unsigned int line=0; line<4; ++line) for (unsigned int dof=0; dof!=fe.dofs_per_vertex; ++dof) selected_dofs[face->line(line)->child(0)->vertex_dof_index(1,dof)] = true; - // next the dofs on the lines interior - // to the face; the order of these - // lines is laid down in the - // FiniteElement class documentation + // next the dofs on the lines interior to the face; the + // order of these lines is laid down in the FiniteElement + // class documentation for (unsigned int dof=0; dofchild(0)->line(1)->dof_index(dof)] = true; for (unsigned int dof=0; dofline(line)->child(child)->dof_index(dof)] = true; - // finally, for the dofs interior - // to the four child faces + // finally, for the dofs interior to the four child faces for (unsigned int child=0; child<4; ++child) for (unsigned int dof=0; dof!=fe.dofs_per_quad; ++dof) selected_dofs[face->child(child)->dof_index(dof)] = true; @@ -4546,9 +3958,8 @@ namespace DoFTools 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(); @@ -4585,12 +3996,9 @@ namespace DoFTools // 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; @@ -4624,12 +4032,9 @@ namespace DoFTools // 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; @@ -4681,8 +4086,7 @@ namespace DoFTools n_selected_dofs += std::count (dofs_by_component.begin(), dofs_by_component.end(), i); - // First count the number of dofs - // in the current component. + // First count the number of dofs in the current component. constant_modes.resize (n_components_selected, std::vector(n_selected_dofs, false)); std::vector component_list (n_components, 0); @@ -4721,9 +4125,8 @@ namespace DoFTools get_subdomain_association (const DH &dof_handler, std::vector &subdomain_association) { - // if the Triangulation is distributed, the - // only thing we can usefully ask is for - // its locally owned subdomain + // if the Triangulation is distributed, the only thing we can usefully + // ask is for its locally owned subdomain Assert ((dynamic_cast *> (&dof_handler.get_tria()) == 0), @@ -4742,16 +4145,12 @@ namespace DoFTools 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(); @@ -4766,16 +4165,10 @@ namespace DoFTools local_dof_indices.resize (dofs_per_cell); cell->get_dof_indices (local_dof_indices); - // set subdomain ids. if dofs - // already have their values - // set then they must be on - // partition interfaces. in - // that case randomly assign - // them to either the previous - // association or the current - // one, where we take "random" - // to be "once this way once - // that way" + // set subdomain ids. if dofs already have their values set then + // they must be on partition interfaces. in that case randomly + // assign them to either the previous association or the current + // one, where we take "random" to be "once this way once that way" for (unsigned int i=0; i 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; typename DH::active_cell_iterator @@ -4884,9 +4274,8 @@ namespace DoFTools dof_handler.get_fe().n_components())); std::fill (n_dofs_on_subdomain.begin(), n_dofs_on_subdomain.end(), 0); - // in debug mode, make sure that there are - // some cells at least with this subdomain - // id + // in debug mode, make sure that there are some cells at least with + // this subdomain id #ifdef DEBUG { bool found = false; @@ -4923,10 +4312,10 @@ namespace DoFTools namespace internal { -//TODO: why is this function so complicated? It would be nice to have -// comments that explain why we can't just loop over all components -// and count the entries in dofs_by_component that have this -// component's index + // TODO: why is this function so complicated? It would be nice to have + // comments that explain why we can't just loop over all components and + // count the entries in dofs_by_component that have this component's + // index template void resolve_components (const FiniteElement &fe, @@ -4955,9 +4344,9 @@ namespace DoFTools dofs_by_component.end(), component); - // if we have non-primitive FEs and want all - // components to show the number of dofs, need - // to copy the result to those components + // if we have non-primitive FEs and want all components + // to show the number of dofs, need to copy the result to + // those components if (!base.is_primitive() && !only_once) for (unsigned int dd=1; dd void resolve_components (const hp::FECollection &fe_collection, @@ -4976,9 +4366,9 @@ namespace DoFTools 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 dofs_by_component (dof_handler.n_locally_owned_dofs()); internal::get_component_association (dof_handler, ComponentMask(), dofs_by_component); @@ -5089,10 +4476,8 @@ namespace DoFTools only_once, dofs_per_component, component); Assert (n_components == component, ExcInternalError()); - // finally sanity check. this is - // only valid if the finite element - // is actually primitive, so - // exclude other elements from this + // finally sanity check. this is only valid if the finite element is + // actually primitive, so exclude other elements from this Assert ((internal::all_elements_are_primitive(dof_handler.get_fe()) == false) || (std::accumulate (dofs_per_component.begin(), @@ -5136,8 +4521,7 @@ namespace DoFTools const FiniteElement &fe = fe_collection[this_fe]; std::fill (dofs_per_block.begin(), dofs_per_block.end(), 0U); - // If the empty vector was given as - // default argument, set up this + // If the empty vector was given as default argument, set up this // vector as identity. if (target_block.size()==0) { @@ -5160,18 +4544,14 @@ namespace DoFTools AssertDimension (dofs_per_block.size(), n_target_blocks); - // special case for only one - // block. treat this first - // since it does not require any - // computations + // special case for only one block. treat this first since it does + // not require any computations if (n_blocks == 1) { dofs_per_block[0] = dof_handler.n_dofs(); return; } - // otherwise determine the number - // of dofs in each block - // separately. + // otherwise determine the number of dofs in each block separately. std::vector dofs_by_block (dof_handler.n_locally_owned_dofs()); internal::get_block_association (dof_handler, dofs_by_block); @@ -5183,10 +4563,8 @@ namespace DoFTools #ifdef DEAL_II_USE_P4EST #if DEAL_II_COMPILER_SUPPORTS_MPI - // if we are working on a parallel - // mesh, we now need to collect - // this information from all - // processors + // if we are working on a parallel mesh, we now need to collect + // this information from all processors if (const parallel::distributed::Triangulation *tria = (dynamic_cast*> (&dof_handler.get_tria()))) @@ -5220,14 +4598,10 @@ namespace DoFTools namespace { /** - * This is a function that is - * called by the _2 function and - * that operates on a range of - * cells only. It is used to - * split up the whole range of - * cells into chunks which are - * then worked on in parallel, if - * multithreading is available. + * This is a function that is called by the _2 function and that + * operates on a range of cells only. It is used to split up the + * whole range of cells into chunks which are then worked on in + * parallel, if multithreading is available. */ template void @@ -5241,36 +4615,23 @@ namespace DoFTools const typename dealii::DoFHandler::active_cell_iterator &begin, const typename dealii::DoFHandler::active_cell_iterator &end) { - // aliases to the finite elements - // used by the dof handlers: + // aliases to the finite elements used by the dof handlers: const FiniteElement &coarse_fe = coarse_grid.get_fe(); - // for each cell on the parameter grid: - // find out which degrees of freedom on the - // fine grid correspond in which way to - // the degrees of freedom on the parameter - // grid + // for each cell on the parameter grid: find out which degrees of + // freedom on the fine grid correspond in which way to the degrees + // of freedom on the parameter grid // - // since for continuous FEs some - // dofs exist on more than one - // cell, we have to track which - // ones were already visited. the - // problem is that if we visit a - // dof first on one cell and - // compute its weight with respect - // to some global dofs to be - // non-zero, and later visit the - // dof again on another cell and - // (since we are on another cell) - // recompute the weights with - // respect to the same dofs as - // above to be zero now, we have to - // preserve them. we therefore - // overwrite all weights if they - // are nonzero and do not enforce - // zero weights since that might be - // only due to the fact that we are - // on another cell. + // since for continuous FEs some dofs exist on more than one cell, + // we have to track which ones were already visited. the problem is + // that if we visit a dof first on one cell and compute its weight + // with respect to some global dofs to be non-zero, and later visit + // the dof again on another cell and (since we are on another cell) + // recompute the weights with respect to the same dofs as above to + // be zero now, we have to preserve them. we therefore overwrite + // all weights if they are nonzero and do not enforce zero weights + // since that might be only due to the fact that we are on another + // cell. // // example: // coarse grid @@ -5289,31 +4650,20 @@ namespace DoFTools // | | | | | // *--x--y--*--* // - // when on cell 1, we compute the - // weights of dof 'x' to be 1/2 - // from parameter dofs 0 and 1, - // respectively. however, when - // later we are on cell 2, we again - // compute the prolongation of - // shape function 1 restricted to - // cell 2 to the globla grid and - // find that the weight of global - // dof 'x' now is zero. however, we - // should not overwrite the old - // value. + // when on cell 1, we compute the weights of dof 'x' to be 1/2 from + // parameter dofs 0 and 1, respectively. however, when later we are + // on cell 2, we again compute the prolongation of shape function 1 + // restricted to cell 2 to the globla grid and find that the weight + // of global dof 'x' now is zero. however, we should not overwrite + // the old value. // - // we therefore always only set - // nonzero values. why adding up is - // not useful: dof 'y' would get - // weight 1 from parameter dof 1 on - // both cells 1 and 2, but the - // correct weight is nevertheless - // only 1. - - // vector to hold the representation of - // a single degree of freedom on the - // coarse grid (for the selected fe) - // on the fine grid + // we therefore always only set nonzero values. why adding up is + // not useful: dof 'y' would get weight 1 from parameter dof 1 on + // both cells 1 and 2, but the correct weight is nevertheless only + // 1. + + // vector to hold the representation of a single degree of freedom + // on the coarse grid (for the selected fe) on the fine grid const unsigned int n_fine_dofs = weight_mapping.size(); dealii::Vector global_parameter_representation (n_fine_dofs); @@ -5322,14 +4672,12 @@ namespace DoFTools for (cell=begin; cell!=end; ++cell) { - // get the global indices of the - // parameter dofs on this parameter - // grid cell + // get the global indices of the parameter dofs on this + // parameter grid cell cell->get_dof_indices (parameter_dof_indices); - // loop over all dofs on this - // cell and check whether they - // are interesting for us + // loop over all dofs on this cell and check whether they are + // interesting for us for (unsigned int local_dof=0; local_dof set_dof_values_by_interpolation (parameter_dofs[local_parameter_dof], global_parameter_representation); - // now that we've got the global - // representation of each parameter - // dof, we've only got to clobber the - // non-zero entries in that vector and - // store the result + // now that we've got the global representation of each + // parameter dof, we've only got to clobber the non-zero + // entries in that vector and store the result // - // what we have learned: if entry @p{i} - // of the global vector holds the value - // @p{v[i]}, then this is the weight with - // which the present dof contributes - // to @p{i}. there may be several such - // @p{i}s and their weights' sum should - // be one. Then, @p{v[i]} should be - // equal to @p{\sum_j w_{ij} p[j]} with - // @p{p[j]} be the values of the degrees - // of freedom on the coarse grid. we - // can thus compute constraints which - // link the degrees of freedom @p{v[i]} - // on the fine grid to those on the - // coarse grid, @p{p[j]}. Now to use - // these as real constraints, rather - // than as additional equations, we - // have to identify representants - // among the @p{i} for each @p{j}. this will - // be done by simply taking the first - // @p{i} for which @p{w_{ij}==1}. + // what we have learned: if entry @p{i} of the global + // vector holds the value @p{v[i]}, then this is the + // weight with which the present dof contributes to + // @p{i}. there may be several such @p{i}s and their + // weights' sum should be one. Then, @p{v[i]} should be + // equal to @p{\sum_j w_{ij} p[j]} with @p{p[j]} be the + // values of the degrees of freedom on the coarse grid. + // we can thus compute constraints which link the degrees + // of freedom @p{v[i]} on the fine grid to those on the + // coarse grid, @p{p[j]}. Now to use these as real + // constraints, rather than as additional equations, we + // have to identify representants among the @p{i} for + // each @p{j}. this will be done by simply taking the + // first @p{i} for which @p{w_{ij}==1}. // - // guard modification of - // the weights array by a - // Mutex. since it should - // happen rather rarely - // that there are several - // threads operating on - // different intergrid - // weights, have only one - // mutex for all of them + // guard modification of the weights array by a Mutex. + // since it should happen rather rarely that there are + // several threads operating on different intergrid + // weights, have only one mutex for all of them static Threads::ThreadMutex mutex; Threads::ThreadMutex::ScopedLock lock (mutex); for (unsigned int i=0; i void @@ -5430,15 +4757,15 @@ namespace DoFTools const std::vector &weight_mapping, std::vector > &weights) { - // simply distribute the range of - // cells to different threads + // simply distribute the range of cells to different threads typedef typename dealii::DoFHandler::active_cell_iterator active_cell_iterator; std::vector > cell_intervals = Threads::split_range (coarse_grid.begin_active(), coarse_grid.end(), multithread_info.n_default_threads); -//TODO: use WorkStream here + // TODO: use WorkStream here + Threads::TaskGroup<> tasks; void (*fun_ptr) (const dealii::DoFHandler &, const unsigned int , @@ -5464,11 +4791,9 @@ namespace DoFTools /** - * This is a helper function that - * is used in the computation of - * integrid constraints. See the - * function for a thorough - * description of how it works. + * This is a helper function that is used in the computation of + * integrid constraints. See the function for a thorough description + * of how it works. */ template unsigned int @@ -5481,8 +4806,7 @@ namespace DoFTools std::vector > &weights, std::vector &weight_mapping) { - // aliases to the finite elements - // used by the dof handlers: + // aliases to the finite elements used by the dof handlers: const FiniteElement &coarse_fe = coarse_grid.get_fe(), &fine_fe = fine_grid.get_fe(); @@ -5493,46 +4817,38 @@ namespace DoFTools // local numbers of dofs const unsigned int fine_dofs_per_cell = fine_fe.dofs_per_cell; - // alias the number of dofs per - // cell belonging to the - // coarse_component which is to be - // the restriction of the fine + // alias the number of dofs per cell belonging to the + // coarse_component which is to be the restriction of the fine // grid: const unsigned int coarse_dofs_per_cell_component = coarse_fe.base_element(coarse_fe.component_to_base_index(coarse_component).first).dofs_per_cell; - // Try to find out whether the - // grids stem from the same coarse - // grid. This is a rather crude - // test, but better than nothing + // Try to find out whether the grids stem from the same coarse + // grid. This is a rather crude test, but better than nothing Assert (coarse_grid.get_tria().n_cells(0) == fine_grid.get_tria().n_cells(0), ExcGridsDontMatch()); - // check whether the map correlates - // the right objects + // check whether the map correlates the right objects Assert (&coarse_to_fine_grid_map.get_source_grid() == &coarse_grid, ExcGridsDontMatch ()); Assert (&coarse_to_fine_grid_map.get_destination_grid() == &fine_grid, ExcGridsDontMatch ()); - // check whether component numbers - // are valid + // check whether component numbers are valid AssertIndexRange (coarse_component,coarse_fe.n_components()); AssertIndexRange (fine_component, fine_fe.n_components()); - // check whether respective finite - // elements are equal + + // check whether respective finite elements are equal Assert (coarse_fe.base_element (coarse_fe.component_to_base_index(coarse_component).first) == fine_fe.base_element (fine_fe.component_to_base_index(fine_component).first), ExcFiniteElementsDontMatch()); #ifdef DEBUG - // if in debug mode, check whether - // the coarse grid is indeed - // coarser everywhere than the fine - // grid + // if in debug mode, check whether the coarse grid is indeed + // coarser everywhere than the fine grid for (typename dealii::DoFHandler::active_cell_iterator cell=coarse_grid.begin_active(); cell != coarse_grid.end(); ++cell) @@ -5540,40 +4856,35 @@ namespace DoFTools ExcGridNotCoarser()); #endif - - /* - * From here on: the term `parameter' refers to the selected component - * on the coarse grid and its analogon on the fine grid. The naming of - * variables containing this term is due to the fact that - * `selected_component' is longer, but also due to the fact that the - * code of this function was initially written for a program where the - * component which we wanted to match between grids was actually the - * `parameter' variable. + * From here on: the term `parameter' refers to the selected + * component on the coarse grid and its analogon on the fine grid. + * The naming of variables containing this term is due to the fact + * that `selected_component' is longer, but also due to the fact + * that the code of this function was initially written for a + * program where the component which we wanted to match between + * grids was actually the `parameter' variable. * - * Likewise, the terms `parameter grid' and `state grid' refer to the - * coarse and fine grids, respectively. + * Likewise, the terms `parameter grid' and `state grid' refer to + * the coarse and fine grids, respectively. * - * Changing the names of variables would in principle be a good idea, - * but would not make things simpler and would be another source of - * errors. If anyone feels like doing so: patches would be welcome! + * Changing the names of variables would in principle be a good + * idea, but would not make things simpler and would be another + * source of errors. If anyone feels like doing so: patches would + * be welcome! */ - // set up vectors of cell-local - // data; each vector represents one - // degree of freedom of the - // coarse-grid variable in the - // fine-grid element + // set up vectors of cell-local data; each vector represents one + // degree of freedom of the coarse-grid variable in the fine-grid + // element std::vector > parameter_dofs (coarse_dofs_per_cell_component, dealii::Vector(fine_dofs_per_cell)); - // for each coarse dof: find its - // position within the fine element - // and set this value to one in the - // respective vector (all other values - // are zero by construction) + // for each coarse dof: find its position within the fine element + // and set this value to one in the respective vector (all other + // values are zero by construction) for (unsigned int local_coarse_dof=0; local_coarse_dof dof_is_interesting (fine_grid.n_dofs(), false); std::vector local_dof_indices (fine_fe.dofs_per_cell); @@ -5634,10 +4942,8 @@ namespace DoFTools { cell->get_dof_indices (local_dof_indices); for (unsigned int i=0; i > &coarse_to_fine_grid_map, ConstraintMatrix &constraints) { - // store the weights with which a dof - // on the parameter grid contributes - // to a dof on the fine grid. see the - // long doc below for more info + // store the weights with which a dof on the parameter grid contributes + // to a dof on the fine grid. see the long doc below for more info // - // allocate as many rows as there are - // parameter dofs on the coarse grid - // and as many columns as there are - // parameter dofs on the fine grid. + // allocate as many rows as there are parameter dofs on the coarse grid + // and as many columns as there are parameter dofs on the fine grid. // - // weight_mapping is used to map the - // global (fine grid) parameter dof + // weight_mapping is used to map the global (fine grid) parameter dof // indices to the columns // - // in the original implementation, - // the weights array was actually - // of FullMatrix type. this - // wasted huge amounts of memory, - // but was fast. nonetheless, since - // the memory consumption was - // quadratic in the number of - // degrees of freedom, this was not - // very practical, so we now use a - // vector of rows of the matrix, - // and in each row a vector of - // pairs (colnum,value). this seems - // like the best tradeoff between - // memory and speed, as it is now - // linear in memory and still fast + // in the original implementation, the weights array was actually of + // FullMatrix type. this wasted huge amounts of memory, but was + // fast. nonetheless, since the memory consumption was quadratic in the + // number of degrees of freedom, this was not very practical, so we now + // use a vector of rows of the matrix, and in each row a vector of + // pairs (colnum,value). this seems like the best tradeoff between + // memory and speed, as it is now linear in memory and still fast // enough. // - // to save some memory and since - // the weights are usually - // (negative) powers of 2, we - // choose the value type of the - // matrix to be @p{float} rather - // than @p{double}. + // to save some memory and since the weights are usually (negative) + // powers of 2, we choose the value type of the matrix to be @p{float} + // rather than @p{double}. std::vector > weights; - // this is this mapping. there is one - // entry for each dof on the fine grid; - // if it is a parameter dof, then its - // value is the column in weights for - // that parameter dof, if it is any - // other dof, then its value is -1, - // indicating an error + // this is this mapping. there is one entry for each dof on the fine + // grid; if it is a parameter dof, then its value is the column in + // weights for that parameter dof, if it is any other dof, then its + // value is -1, indicating an error std::vector weight_mapping; const unsigned int n_parameters_on_fine_grid @@ -5780,9 +5054,8 @@ namespace DoFTools n_fine_dofs = fine_grid.n_dofs(); - // get an array in which we store - // which dof on the coarse grid is - // a parameter and which is not + // get an array in which we store which dof on the coarse grid is a + // parameter and which is not std::vector coarse_dof_is_parameter (coarse_grid.n_dofs()); if (true) { @@ -5792,40 +5065,26 @@ namespace DoFTools extract_dofs (coarse_grid, ComponentMask(mask), coarse_dof_is_parameter); } - // now we know that the weights in - // each row constitute a - // constraint. enter this into the - // constraints object + // now we know that the weights in each row constitute a constraint. + // enter this into the constraints object // - // first task: for each parameter - // dof on the parameter grid, find - // a representant on the fine, - // global grid. this is possible - // since we use conforming finite - // element. we take this - // representant to be the first - // element in this row with weight - // identical to one. the - // representant will become an - // unconstrained degree of freedom, - // while all others will be - // constrained to this dof (and - // possibly others) + // first task: for each parameter dof on the parameter grid, find a + // representant on the fine, global grid. this is possible since we use + // conforming finite element. we take this representant to be the first + // element in this row with weight identical to one. the representant + // will become an unconstrained degree of freedom, while all others + // will be constrained to this dof (and possibly others) std::vector representants(n_coarse_dofs, -1); for (unsigned int parameter_dof=0; parameter_dof 0, ExcInternalError()); - // find the column where the - // representant is mentioned + // find the column where the representant is mentioned std::map::const_iterator i = weights[parameter_dof].begin(); for (; i!=weights[parameter_dof].end(); ++i) if (i->second == 1) @@ -5833,59 +5092,42 @@ namespace DoFTools Assert (i!=weights[parameter_dof].end(), ExcInternalError()); const unsigned int column = i->first; - // now we know in which column of - // weights the representant is, but - // we don't know its global index. get - // it using the inverse operation of - // the weight_mapping + // now we know in which column of weights the representant is, + // but we don't know its global index. get it using the inverse + // operation of the weight_mapping unsigned int global_dof=0; for (; global_dof(column)) break; Assert (global_dof < weight_mapping.size(), ExcInternalError()); - // now enter the representants global - // index into our list + // now enter the representants global index into our list representants[parameter_dof] = global_dof; } else { - // consistency check: if this - // is no parameter dof on the - // coarse grid, then the - // respective row must be - // empty! + // consistency check: if this is no parameter dof on the coarse + // grid, then the respective row must be empty! Assert (weights[parameter_dof].size() == 0, ExcInternalError()); }; - // note for people that want to - // optimize this function: the - // largest part of the computing - // time is spent in the following, - // rather innocent block of - // code. basically, it must be the - // ConstraintMatrix::add_entry call - // which takes the bulk of the - // time, but it is not known to the - // author how to make it faster... + // note for people that want to optimize this function: the largest + // part of the computing time is spent in the following, rather + // innocent block of code. basically, it must be the + // ConstraintMatrix::add_entry call which takes the bulk of the time, + // but it is not known to the author how to make it faster... std::vector > constraint_line; for (unsigned int global_dof=0; global_dofsecond == 1) && (representants[first_used_row] == static_cast(global_dof))) - // dof unconstrained or - // constrained to itself - // (in case this cell is - // mapped to itself, rather - // than to children of + // dof unconstrained or constrained to itself (in case this + // cell is mapped to itself, rather than to children of // itself) continue; } @@ -5946,52 +5185,33 @@ namespace DoFTools const InterGridMap > &coarse_to_fine_grid_map, std::vector > &transfer_representation) { - // store the weights with which a dof - // on the parameter grid contributes - // to a dof on the fine grid. see the - // long doc below for more info + // store the weights with which a dof on the parameter grid contributes + // to a dof on the fine grid. see the long doc below for more info // - // allocate as many rows as there are - // parameter dofs on the coarse grid - // and as many columns as there are - // parameter dofs on the fine grid. + // allocate as many rows as there are parameter dofs on the coarse grid + // and as many columns as there are parameter dofs on the fine grid. // - // weight_mapping is used to map the - // global (fine grid) parameter dof + // weight_mapping is used to map the global (fine grid) parameter dof // indices to the columns // - // in the original implementation, - // the weights array was actually - // of FullMatrix type. this - // wasted huge amounts of memory, - // but was fast. nonetheless, since - // the memory consumption was - // quadratic in the number of - // degrees of freedom, this was not - // very practical, so we now use a - // vector of rows of the matrix, - // and in each row a vector of - // pairs (colnum,value). this seems - // like the best tradeoff between - // memory and speed, as it is now - // linear in memory and still fast + // in the original implementation, the weights array was actually of + // FullMatrix type. this wasted huge amounts of memory, but was + // fast. nonetheless, since the memory consumption was quadratic in the + // number of degrees of freedom, this was not very practical, so we now + // use a vector of rows of the matrix, and in each row a vector of + // pairs (colnum,value). this seems like the best tradeoff between + // memory and speed, as it is now linear in memory and still fast // enough. // - // to save some memory and since - // the weights are usually - // (negative) powers of 2, we - // choose the value type of the - // matrix to be @p{float} rather - // than @p{double}. + // to save some memory and since the weights are usually (negative) + // powers of 2, we choose the value type of the matrix to be @p{float} + // rather than @p{double}. std::vector > weights; - // this is this mapping. there is one - // entry for each dof on the fine grid; - // if it is a parameter dof, then its - // value is the column in weights for - // that parameter dof, if it is any - // other dof, then its value is -1, - // indicating an error + // this is this mapping. there is one entry for each dof on the fine + // grid; if it is a parameter dof, then its value is the column in + // weights for that parameter dof, if it is any other dof, then its + // value is -1, indicating an error std::vector weight_mapping; internal::compute_intergrid_weights_1 (coarse_grid, coarse_component, @@ -5999,21 +5219,18 @@ namespace DoFTools coarse_to_fine_grid_map, weights, weight_mapping); - // now compute the requested - // representation + // now compute the requested representation const unsigned int n_global_parm_dofs = std::count_if (weight_mapping.begin(), weight_mapping.end(), std::bind2nd (std::not_equal_to (), -1)); - // first construct the inverse - // mapping of weight_mapping + // first construct the inverse mapping of weight_mapping std::vector inverse_weight_mapping (n_global_parm_dofs, DoFHandler::invalid_dof_index); for (unsigned int i=0; ihas_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 typename DH::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) @@ -6155,8 +5363,7 @@ namespace DoFTools for (unsigned int fe_index = 0; fe_index < fe_collection.size(); ++fe_index) { - // check whether every fe in the collection - // has support points + // check whether every fe in the collection has support points Assert(fe_collection[fe_index].has_support_points(), typename FiniteElement::ExcFEHasNoSupportPoints()); q_coll_dummy.push_back( @@ -6164,18 +5371,13 @@ namespace DoFTools fe_collection[fe_index].get_unit_support_points())); } - // now loop over all cells and - // enquire the support points on - // each of these. we use dummy - // quadrature formulas where the - // quadrature points are located at - // the unit support points to - // enquire the location of the - // support points in real space + // Now loop over all cells and enquire the support points on each + // of these. we use dummy quadrature formulas where the quadrature + // points are located at the unit support points to enquire the + // location of the support points in real space. // - // the weights of the quadrature - // rule have been set to invalid values - // by the used constructor. + // The weights of the quadrature rule have been set to invalid + // values by the used constructor. hp::FEValues hp_fe_values(mapping, fe_collection, q_coll_dummy, update_quadrature_points); typename DH::active_cell_iterator cell = @@ -6237,8 +5439,8 @@ namespace DoFTools ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + // Let the internal function do all the work, just make sure that it + // gets a MappingCollection const hp::MappingCollection mapping_collection(mapping); internal::map_dofs_to_support_points (mapping_collection, @@ -6261,8 +5463,8 @@ namespace DoFTools ExcMessage ("This function can not be used with distributed triangulations." "See the documentation for more information.")); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + // Let the internal function do all the work, just make sure that it + // gets a MappingCollection internal::map_dofs_to_support_points (mapping, dof_handler, support_points); @@ -6277,8 +5479,8 @@ namespace DoFTools { support_points.clear(); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + // Let the internal function do all the work, just make sure that it + // gets a MappingCollection const hp::MappingCollection mapping_collection(mapping); internal::map_dofs_to_support_points (mapping_collection, @@ -6295,8 +5497,8 @@ namespace DoFTools { support_points.clear(); - //Let the internal function do all the work, - //just make sure that it gets a MappingCollection + // Let the internal function do all the work, just make sure that it + // gets a MappingCollection internal::map_dofs_to_support_points (mapping, dof_handler, support_points); @@ -6401,21 +5603,16 @@ namespace DoFTools || (face->boundary_indicator() == boundary_indicator))) { - // get indices and physical - // location on this face + // get indices and physical location on this face face_dofs.resize (fe.dofs_per_face); face->get_dof_indices (face_dofs, cell->active_fe_index()); - // enter those dofs into the list - // that match the component + // enter those dofs into the list that match the component // signature. for (unsigned int i=0; iget_fe().get_nonzero_components (i); bool nonzero = false; @@ -6503,9 +5700,7 @@ namespace DoFTools if (interior_only) { - // Exclude degrees of - // freedom on faces - // opposite to the vertex + // Exclude degrees of freedom on faces opposite to the vertex exclude.resize(fe.dofs_per_cell); std::fill(exclude.begin(), exclude.end(), false); const unsigned int dpf = fe.dofs_per_face; @@ -6551,9 +5746,8 @@ namespace DoFTools for (unsigned int child=0; childn_children(); ++child) { cell = pcell->child(child); - // For hp, only this line - // here would have to be - // replaced. + + // For hp, only this line here would have to be replaced. const FiniteElement &fe = dof_handler.get_fe(); const unsigned int n_dofs = fe.dofs_per_cell; indices.resize(n_dofs); @@ -6563,9 +5757,7 @@ namespace DoFTools if (interior_dofs_only) { - // Eliminate dofs on - // faces of the child - // which are on faces + // Eliminate dofs on faces of the child which are on faces // of the parent const unsigned int dpf = fe.dofs_per_face; @@ -6575,11 +5767,9 @@ namespace DoFTools for (unsigned int i=0; i::faces_per_cell; ++face) if (cell->at_boundary(face)) @@ -6610,23 +5800,21 @@ namespace DoFTools typename DH::cell_iterator cell; typename DH::cell_iterator endc = dof_handler.end(level); - // Vector mapping from vertex - // index in the triangulation to - // consecutive block indices on - // this level - // The number of cells at a vertex + // Vector mapping from vertex index in the triangulation to consecutive + // block indices on this level The number of cells at a vertex std::vector vertex_cell_count(dof_handler.get_tria().n_vertices(), 0); + // Is a vertex at the boundary? std::vector vertex_boundary(dof_handler.get_tria().n_vertices(), false); + std::vector vertex_mapping(dof_handler.get_tria().n_vertices(), numbers::invalid_unsigned_int); - // Estimate for the number of - // dofs at this point + + // Estimate for the number of dofs at this point std::vector vertex_dof_count(dof_handler.get_tria().n_vertices(), 0); - // Identify all vertices active - // on this level and remember - // some data about them + // Identify all vertices active on this level and remember some data + // about them for (cell=dof_handler.begin(level); cell != endc; ++cell) for (unsigned int v=0; v::vertices_per_cell; ++v) { @@ -6643,20 +5831,16 @@ namespace DoFTools vertex_boundary[vg] = true; } } - // From now on, only vertices - // with positive dof count are - // "in". + // From now on, only vertices with positive dof count are "in". - // Remove vertices at boundaries - // or in corners + // Remove vertices at boundaries or in corners for (unsigned int vg=0; vg indices; @@ -6696,9 +5877,8 @@ namespace DoFTools if (interior_only) { - // Exclude degrees of - // freedom on faces - // opposite to the vertex + // Exclude degrees of freedom on faces opposite to the + // vertex exclude.resize(fe.dofs_per_cell); std::fill(exclude.begin(), exclude.end(), false); const unsigned int dpf = fe.dofs_per_face;