From 1d56ca72ec50c1c36944abadf9966762a937094d Mon Sep 17 00:00:00 2001 From: wolf Date: Wed, 9 May 2001 16:11:24 +0000 Subject: [PATCH] Introduce new function DoFTools::compute_intergrid_transfer_representation. Factor out common code with existing functions. git-svn-id: https://svn.dealii.org/trunk@4583 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/deal.II/include/dofs/dof_tools.h | 46 +- deal.II/deal.II/source/dofs/dof_tools.cc | 616 +++++++++++++---------- 2 files changed, 394 insertions(+), 268 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_tools.h b/deal.II/deal.II/include/dofs/dof_tools.h index 8d5222cfea..ebc867ed5e 100644 --- a/deal.II/deal.II/include/dofs/dof_tools.h +++ b/deal.II/deal.II/include/dofs/dof_tools.h @@ -807,7 +807,6 @@ class DoFTools * only need one object of this * type. */ -//TODO:[WB] document last argument template static void compute_intergrid_constraints (const DoFHandler &coarse_grid, @@ -815,9 +814,18 @@ class DoFTools const DoFHandler &fine_grid, const unsigned int fine_component, const InterGridMap &coarse_to_fine_grid_map, - ConstraintMatrix &constraints, - std::vector > *transfer_representation = 0); + ConstraintMatrix &constraints); +//TODO:[WB] document this function. put it in the news file + template + static void + compute_intergrid_transfer_representation (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap &coarse_to_fine_grid_map, + std::vector > &transfer_representation); + /** * Create a mapping from degree * of freedom indices to the @@ -966,6 +974,24 @@ class DoFTools DeclException0 (ExcInvalidBoundaryIndicator); private: + + /** + * 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 + static unsigned int + compute_intergrid_weights_1 (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap &coarse_to_fine_grid_map, + std::vector > &weights, + std::vector &weight_mapping); + /** * This is a helper function that * is used in the computation of @@ -975,12 +1001,12 @@ class DoFTools */ template static void - compute_intergrid_weights (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights); + compute_intergrid_weights_2 (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights); /** * This is a function that is @@ -994,7 +1020,7 @@ class DoFTools */ template static void - compute_intergrid_weights_1 (const DoFHandler &coarse_grid, + compute_intergrid_weights_3 (const DoFHandler &coarse_grid, const unsigned int coarse_component, const InterGridMap &coarse_to_fine_grid_map, const std::vector > ¶meter_dofs, diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index 7133f80b2d..97590cf689 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -1170,8 +1170,331 @@ DoFTools::compute_intergrid_constraints (const DoFHandler &coa const DoFHandler &fine_grid, const unsigned int fine_component, const InterGridMap &coarse_to_fine_grid_map, - ConstraintMatrix &constraints, - std::vector > *transfer_representation) + 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 + // + // 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 + // 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 + // 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}. + 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 + std::vector weight_mapping; + + const unsigned int n_parameters_on_fine_grid + = compute_intergrid_weights_1 (coarse_grid, coarse_component, fine_grid, fine_component, + coarse_to_fine_grid_map, weights, weight_mapping); + + // global numbers of dofs + const unsigned int n_coarse_dofs = coarse_grid.n_dofs(), + 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 + std::vector coarse_dof_is_parameter (coarse_grid.n_dofs()); + if (true) + { + std::vector mask (coarse_grid.get_fe().n_components(), + false); + mask[coarse_component] = true; + extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); + }; + + // 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) + 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 + std::map::const_iterator i = weights[parameter_dof].begin(); + for (; i!=weights[parameter_dof].end(); ++i) + if (i->second == 1) + break; + 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 + 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 + 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! + 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... + std::vector > constraint_line; + for (unsigned int global_dof=0; global_dof::const_iterator col_entry; + for (; first_used_rowsecond == 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 + // itself) + continue; + }; + + + // otherwise enter all constraints + constraints.add_line (global_dof); + + constraint_line.clear (); + for (unsigned int row=first_used_row; row::const_iterator + j = weights[row].find(col); + if ((j != weights[row].end()) && (j->second != 0)) + constraint_line.push_back (std::make_pair(representants[row], + j->second)); + }; + + constraints.add_entries (global_dof, constraint_line); + }; +}; + + + +template +void +DoFTools:: +compute_intergrid_transfer_representation (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + 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 + // + // 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 + // 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 + // 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}. + 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 + std::vector weight_mapping; + + compute_intergrid_weights_1 (coarse_grid, coarse_component, fine_grid, fine_component, + coarse_to_fine_grid_map, weights, weight_mapping); + + // now compute the requested + // representation + transfer_representation.clear (); + transfer_representation.resize (weights.size()); + + 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 + std::vector inverse_weight_mapping (n_global_parm_dofs, + DoFHandler::invalid_dof_index); + for (unsigned int i=0; i(-1)) + { + Assert (parameter_dof < n_global_parm_dofs, ExcInternalError()); + Assert (inverse_weight_mapping[parameter_dof] == DoFHandler::invalid_dof_index, + ExcInternalError()); + + inverse_weight_mapping[parameter_dof] = i; + }; + }; + + // next copy over weights array + // and replace respective + // numbers + const unsigned int n_coarse_dofs = coarse_grid.n_dofs(); + for (unsigned int i=0; i::const_iterator j = weights[i].begin(); + for (; j!=weights[i].end(); ++j) + transfer_representation[i][inverse_weight_mapping[j->first]] = j->second; + }; +}; + + + +template +unsigned int +DoFTools::compute_intergrid_weights_1 (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const DoFHandler &fine_grid, + const unsigned int fine_component, + const InterGridMap &coarse_to_fine_grid_map, + std::vector > &weights, + std::vector &weight_mapping) { // aliases to the finite elements // used by the dof handlers: @@ -1297,74 +1620,19 @@ DoFTools::compute_intergrid_constraints (const DoFHandler &coa dof_is_interesting[local_dof_indices[i]] = true; }; - n_parameters_on_fine_grid = count (dof_is_interesting.begin(), - dof_is_interesting.end(), - true); - }; - - // 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) - { - std::vector mask (coarse_grid.get_fe().n_components(), - false); - mask[coarse_component] = true; - extract_dofs (coarse_grid, mask, coarse_dof_is_parameter); - }; + n_parameters_on_fine_grid = std::count (dof_is_interesting.begin(), + dof_is_interesting.end(), + true); + }; - - - // 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. - // - // 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 - // 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}. - std::vector > weights(n_coarse_dofs); - // 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 (n_fine_dofs, -1); + // set up the weights mapping + weights.clear (); + weights.resize (n_coarse_dofs); - // set up this mapping + weight_mapping.clear (); + weight_mapping.resize (n_fine_dofs, -1); + if (true) { std::vector local_dof_indices(fine_fe.dofs_per_cell); @@ -1403,9 +1671,10 @@ DoFTools::compute_intergrid_constraints (const DoFHandler &coa // you want to read more // information on the algorithm // used. - compute_intergrid_weights (coarse_grid, coarse_component, - coarse_to_fine_grid_map, parameter_dofs, - weight_mapping, weights); + compute_intergrid_weights_2 (coarse_grid, coarse_component, + coarse_to_fine_grid_map, parameter_dofs, + weight_mapping, weights); + // ok, now we have all weights for each // dof on the fine grid. if in debug @@ -1440,197 +1709,21 @@ DoFTools::compute_intergrid_constraints (const DoFHandler &coa }; #endif - // if the user wants to have a - // representation of the transfer - // matrix, the provide it - if (transfer_representation != 0) - { - transfer_representation->clear (); - transfer_representation->resize (weights.size()); - - 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 - std::vector inverse_weight_mapping (n_global_parm_dofs, - DoFHandler::invalid_dof_index); - for (unsigned int i=0; i(-1)) - { - Assert (parameter_dof < n_global_parm_dofs, ExcInternalError()); - Assert (inverse_weight_mapping[parameter_dof] == DoFHandler::invalid_dof_index, - ExcInternalError()); - - inverse_weight_mapping[parameter_dof] = i; - }; - }; - - // next copy over weights array - // and replace respective - // numbers - for (unsigned int i=0; i::const_iterator j = weights[i].begin(); - for (; j!=weights[i].end(); ++j) - (*transfer_representation)[i][inverse_weight_mapping[j->first]] = j->second; - }; - }; - - // 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) - 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 - std::map::const_iterator i = weights[parameter_dof].begin(); - for (; i!=weights[parameter_dof].end(); ++i) - if (i->second == 1) - break; - 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 - 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 - 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! - 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... - std::vector > constraint_line; - for (unsigned int global_dof=0; global_dof::const_iterator col_entry; - for (; first_used_rowsecond == 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 - // itself) - continue; - }; - - - // otherwise enter all constraints - constraints.add_line (global_dof); - - constraint_line.clear (); - for (unsigned int row=first_used_row; row::const_iterator - j = weights[row].find(col); - if ((j != weights[row].end()) && (j->second != 0)) - constraint_line.push_back (std::make_pair(representants[row], - j->second)); - }; - - constraints.add_entries (global_dof, constraint_line); - }; + return n_parameters_on_fine_grid; }; + template void -DoFTools::compute_intergrid_weights (const DoFHandler &coarse_grid, - const unsigned int coarse_component, - const InterGridMap &coarse_to_fine_grid_map, - const std::vector > ¶meter_dofs, - const std::vector &weight_mapping, - std::vector > &weights) +DoFTools::compute_intergrid_weights_2 (const DoFHandler &coarse_grid, + const unsigned int coarse_component, + const InterGridMap &coarse_to_fine_grid_map, + const std::vector > ¶meter_dofs, + const std::vector &weight_mapping, + std::vector > &weights) { // simply distribute the range of // cells to different threads @@ -1643,7 +1736,7 @@ DoFTools::compute_intergrid_weights (const DoFHandler &coarse_ Threads::ThreadManager thread_manager; for (unsigned int i=0; i) + Threads::encapsulate (&DoFTools::template compute_intergrid_weights_3) .collect_args (coarse_grid, coarse_component, coarse_to_fine_grid_map, parameter_dofs, weight_mapping, weights, @@ -1658,7 +1751,7 @@ DoFTools::compute_intergrid_weights (const DoFHandler &coarse_ template void -DoFTools::compute_intergrid_weights_1 (const DoFHandler &coarse_grid, +DoFTools::compute_intergrid_weights_3 (const DoFHandler &coarse_grid, const unsigned int coarse_component, const InterGridMap &coarse_to_fine_grid_map, const std::vector > ¶meter_dofs, @@ -2121,8 +2214,15 @@ DoFTools::compute_intergrid_constraints (const DoFHandler &, const DoFHandler &, const unsigned int , const InterGridMap &, - ConstraintMatrix &, - std::vector > *); + ConstraintMatrix &); +template +void +DoFTools::compute_intergrid_transfer_representation (const DoFHandler &, + const unsigned int , + const DoFHandler &, + const unsigned int , + const InterGridMap &, + std::vector > &); -- 2.39.5