From: kronbichler Date: Mon, 6 Apr 2009 15:48:15 +0000 (+0000) Subject: Make the local_to_global functions a bit easier to read. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d85b0e858e86f5c005cf4adcb77758ef8b433132;p=dealii-svn.git Make the local_to_global functions a bit easier to read. git-svn-id: https://svn.dealii.org/trunk@18559 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/vectors.templates.h b/deal.II/deal.II/include/numerics/vectors.templates.h index 350ad1504f..9ac535b381 100644 --- a/deal.II/deal.II/include/numerics/vectors.templates.h +++ b/deal.II/deal.II/include/numerics/vectors.templates.h @@ -467,8 +467,7 @@ void VectorTools::project (const Mapping &mapping, // still be fast { CompressedSimpleSparsityPattern csp (dof.n_dofs(), dof.n_dofs()); - DoFTools::make_sparsity_pattern (dof, csp); - constraints.condense(csp); + DoFTools::make_sparsity_pattern (dof, csp, constraints); sparsity.copy_from (csp); } diff --git a/deal.II/deal.II/source/dofs/dof_tools.cc b/deal.II/deal.II/source/dofs/dof_tools.cc index b015033dd4..974267dbbb 100644 --- a/deal.II/deal.II/source/dofs/dof_tools.cc +++ b/deal.II/deal.II/source/dofs/dof_tools.cc @@ -132,7 +132,7 @@ DoFTools::make_sparsity_pattern ( ExcDimensionMismatch(couplings.n_cols(), dof.get_fe().n_components())); 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 @@ -141,43 +141,55 @@ DoFTools::make_sparsity_pattern ( // functions, which takes some additional // thought std::vector > dof_mask(fe_collection.size()); - for (unsigned int f=0; f dofs_on_this_cell(fe_collection.max_dofs_per_cell()); diff --git a/deal.II/doc/authors.html b/deal.II/doc/authors.html index 13139b9977..6e2d573a0f 100644 --- a/deal.II/doc/authors.html +++ b/deal.II/doc/authors.html @@ -103,6 +103,7 @@ contact with them, send email to authors at dealii.org.
  • Martin Kronbichler: Parts of step-22 and step-31 tutorial programs, interfaces to Trilinos, + fast implementation of local_to_global functions in ConstraintMatrix, and many enhancements in random places.
  • Tobias Leicht: diff --git a/deal.II/lac/include/lac/constraint_matrix.templates.h b/deal.II/lac/include/lac/constraint_matrix.templates.h index 307417bac8..711d96334e 100644 --- a/deal.II/lac/include/lac/constraint_matrix.templates.h +++ b/deal.II/lac/include/lac/constraint_matrix.templates.h @@ -963,11 +963,15 @@ namespace internals { struct distributing { - distributing (const unsigned int local_row = deal_II_numbers::invalid_unsigned_int); + distributing (const unsigned int global_row = deal_II_numbers::invalid_unsigned_int, + const unsigned int local_row = deal_II_numbers::invalid_unsigned_int); + unsigned int global_row; unsigned int local_row; std::vector > constraints; }; - distributing::distributing (const unsigned int local_row) : + distributing::distributing (const unsigned int global_row, + const unsigned int local_row) : + global_row (global_row), local_row (local_row) {} // a version of std::lower_bound for a @@ -975,20 +979,20 @@ namespace internals // without taking into account the second // argument of the pair. inline - std::vector >::iterator - lower_bound (std::vector >::iterator begin, - std::vector >::iterator end, + std::vector::iterator + lower_bound (std::vector::iterator begin, + std::vector::iterator end, unsigned int val) { unsigned int length = end - begin; unsigned int half; - std::vector >::iterator middle; + std::vector::iterator middle; while (length > 0) { half = length >> 1; middle = begin + half; - if (middle->first < val) + if (middle->global_row < val) { begin = middle; ++begin; @@ -1009,19 +1013,18 @@ namespace internals // and much faster. inline void - insert_index (std::vector > &my_indices, + insert_index (std::vector &my_indices, const unsigned int row, const unsigned int local_row, const std::pair constraint = std::make_pair(0,0)) { - typedef std::vector >::iterator - index_iterator; + typedef std::vector::iterator index_iterator; index_iterator pos, pos1; - if (my_indices.size() == 0 || my_indices.back().first < row) + if (my_indices.size() == 0 || my_indices.back().global_row < row) { - my_indices.push_back(std::make_pair(row,distributing())); + my_indices.push_back(distributing(row)); pos1 = my_indices.end()-1; } else @@ -1030,30 +1033,44 @@ namespace internals my_indices.end(), row); - if (pos->first == row) + if (pos->global_row == row) pos1 = pos; else - pos1 = my_indices.insert(pos,std::make_pair - (row,distributing())); + pos1 = my_indices.insert(pos, distributing(row)); } if (local_row == deal_II_numbers::invalid_unsigned_int) - pos1->second.constraints.push_back (constraint); + pos1->constraints.push_back (constraint); else - pos1->second.local_row = local_row; + pos1->local_row = local_row; } + + // this sort algorithm sorts a vector of + // distributing elements, but does not + // take the constraints into + // account. this means that in case that + // constraints are already inserted, this + // function does not work as + // expected. shellsort is very fast in + // case the indices are already sorted + // (which is the usual case with DG + // elements), and not too slow in other + // cases inline void - list_shellsort (std::vector > &my_indices) + list_shellsort (std::vector &my_indices) { - // now sort the actual dofs using a shell - // sort which is very fast in case the - // indices are already sorted (which is - // the usual case with DG elements) unsigned int i, j, j2, temp, templ, istep; unsigned step; + // in debug mode, check whether the + // constraints are really empty. +#ifdef DEBUG + for (unsigned int i=0; i 0) @@ -1063,19 +1080,19 @@ namespace internals istep = step; j = i; j2 = j-istep; - temp = my_indices[i].first; - templ = my_indices[i].second.local_row; - if (my_indices[j2].first > temp) + temp = my_indices[i].global_row; + templ = my_indices[i].local_row; + if (my_indices[j2].global_row > temp) { - while ((j >= istep) && (my_indices[j2].first > temp)) + while ((j >= istep) && (my_indices[j2].global_row > temp)) { - my_indices[j].first = my_indices[j2].first; - my_indices[j].second.local_row = my_indices[j2].second.local_row; + my_indices[j].global_row = my_indices[j2].global_row; + my_indices[j].local_row = my_indices[j2].local_row; j = j2; j2 -= istep; } - my_indices[j].first = temp; - my_indices[j].second.local_row = templ; + my_indices[j].global_row = temp; + my_indices[j].local_row = templ; } } step = step>>1; @@ -1152,7 +1169,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // constraints). Choosing an STL map or // anything else I know of would be much // more expensive here! - std::vector > my_indices (n_local_dofs); + std::vector my_indices (n_local_dofs); std::vector > constraint_lines; constraint_lines.reserve(n_local_dofs); @@ -1172,8 +1189,8 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (constraint_line_exists.size() <= local_dof_indices[i] || constraint_line_exists[local_dof_indices[i]] == false) { - my_indices[added_rows].first = local_dof_indices[i]; - my_indices[added_rows].second.local_row = i; + my_indices[added_rows].global_row = local_dof_indices[i]; + my_indices[added_rows].local_row = i; ++added_rows; continue; } @@ -1257,8 +1274,8 @@ distribute_local_to_global (const FullMatrix &local_matrix, for (unsigned int i=0; i &local_matrix, for (unsigned int j=0; j < n_actual_dofs; ++j) { - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; Assert(loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); if (local_matrix(loc_row,loc_col) != 0) { vals[col_counter] = local_matrix(loc_row,loc_col); - cols[col_counter] = my_indices[j].first; + cols[col_counter] = my_indices[j].global_row; col_counter++; } } @@ -1309,7 +1326,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, for (unsigned int j=0; j < n_actual_dofs; ++j) { double col_val; - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; // case 1: row has direct contribution in // local matrix @@ -1333,11 +1350,11 @@ distribute_local_to_global (const FullMatrix &local_matrix, // account for indirect contributions by // constraints - for (unsigned int p=0; p &local_matrix, // constraints in row, going trough the // direct and indirect references in the // given column. - for (unsigned int q=0; q= 0 && loc_col < n_local_dofs, ExcInternalError()); - add_this = local_matrix(my_indices[i].second.constraints[q].first, + add_this = local_matrix(my_indices[i].constraints[q].first, loc_col); } else add_this = 0; - for (unsigned int p=0; p &local_matrix, // append it to the array of values. if (col_val != 0) { - cols[col_counter] = my_indices[j].first; + cols[col_counter] = my_indices[j].global_row; vals[col_counter] = col_val; col_counter++; } @@ -1402,14 +1419,14 @@ distribute_local_to_global (const FullMatrix &local_matrix, local_matrix(loc_row,constraint_lines[i].first); } - for (unsigned int q=0; q < my_indices[i].second.constraints.size(); ++q) + for (unsigned int q=0; q < my_indices[i].constraints.size(); ++q) { - const unsigned int loc_row_q = my_indices[i].second.constraints[q].first; + const unsigned int loc_row_q = my_indices[i].constraints[q].first; double add_this = local_vector (loc_row_q); for (unsigned int k=0; kinhomogeneity * local_matrix(loc_row_q,constraint_lines[k].first); - val += add_this * my_indices[i].second.constraints[q].second; + val += add_this * my_indices[i].constraints[q].second; } } } @@ -1471,7 +1488,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, average_diagonal += std::fabs (local_matrix(i,i)); average_diagonal /= n_local_dofs; - std::vector > my_indices (n_local_dofs); + std::vector my_indices (n_local_dofs); std::vector > constraint_lines; constraint_lines.reserve(n_local_dofs); @@ -1483,8 +1500,8 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (constraint_line_exists.size() <= local_dof_indices[i] || constraint_line_exists[local_dof_indices[i]] == false) { - my_indices[added_rows].first = local_dof_indices[i]; - my_indices[added_rows].second.local_row = i; + my_indices[added_rows].global_row = local_dof_indices[i]; + my_indices[added_rows].local_row = i; ++added_rows; continue; } @@ -1533,7 +1550,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, std::vector localized_indices (n_actual_dofs); for (unsigned int i=0; i &local_matrix, { unsigned int col_counter = 0; std::vector::iterator block_it = block_ends.begin(); - const unsigned int row = my_indices[i].first; - const unsigned int loc_row = my_indices[i].second.local_row; + const unsigned int row = my_indices[i].global_row; + const unsigned int loc_row = my_indices[i].local_row; double val = 0; if (have_indirect_rows == false) @@ -1593,7 +1610,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, ExcInternalError()); } - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; Assert(loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); @@ -1607,7 +1624,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, // now comes the hack that sets the // correct end of block in case we work // with block matrices. - while (n_actual_dofs == *block_it && block_it != block_ends.end()) + while (block_it != block_ends.end() && n_actual_dofs == *block_it) { actual_block_ends[block_it-block_ends.begin()] = col_counter; ++block_it; @@ -1638,7 +1655,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, } double col_val; - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; if (loc_row != deal_II_numbers::invalid_unsigned_int) { @@ -1654,35 +1671,35 @@ distribute_local_to_global (const FullMatrix &local_matrix, else col_val = 0; - for (unsigned int p=0; p= 0 && loc_col < n_local_dofs, ExcInternalError()); - add_this = local_matrix(my_indices[i].second.constraints[q].first, + add_this = local_matrix(my_indices[i].constraints[q].first, loc_col); } else add_this = 0; - for (unsigned int p=0; p &local_matrix, local_matrix(loc_row,constraint_lines[i].first); } - for (unsigned int q=0; q < my_indices[i].second.constraints.size(); ++q) + for (unsigned int q=0; q < my_indices[i].constraints.size(); ++q) { - const unsigned int loc_row_q = my_indices[i].second.constraints[q].first; + const unsigned int loc_row_q = my_indices[i].constraints[q].first; double add_this = local_vector (loc_row_q); for (unsigned int k=0; kinhomogeneity * local_matrix(loc_row_q,constraint_lines[k].first); - val += add_this * my_indices[i].second.constraints[q].second; + val += add_this * my_indices[i].constraints[q].second; } } } @@ -1778,7 +1795,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, ExcDimensionMismatch(dof_mask.n_cols(), n_local_dofs)); } - std::vector > my_indices (n_local_dofs); + std::vector my_indices (n_local_dofs); std::vector > constraint_lines; constraint_lines.reserve(n_local_dofs); @@ -1798,8 +1815,8 @@ add_entries_local_to_global (const std::vector &local_dof_indices, if (constraint_line_exists.size() <= local_dof_indices[i] || constraint_line_exists[local_dof_indices[i]] == false) { - my_indices[added_rows].first = local_dof_indices[i]; - my_indices[added_rows].second.local_row = i; + my_indices[added_rows].global_row = local_dof_indices[i]; + my_indices[added_rows].local_row = i; ++added_rows; continue; } @@ -1881,12 +1898,33 @@ add_entries_local_to_global (const std::vector &local_dof_indices, // sparsity pattern. std::vector cols (n_actual_dofs); - // now do the actual job. + // easy case - we add all indices (i.e., + // the dof_mask is not active). so to + // each global row we add all the + // indices. + if (dof_mask_is_active == false) + { + for (unsigned int i=0; i &local_dof_indices, for (unsigned int j=0; j < n_actual_dofs; ++j) { - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; Assert(loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); - bool add_this = true; - if (dof_mask_is_active == true) - if (dof_mask[loc_row][loc_col] == false) - add_this = false; - - if (add_this == true) - cols[col_counter++] = my_indices[j].first; - + if (dof_mask[loc_row][loc_col] == true) + cols[col_counter++] = my_indices[j].global_row; } } @@ -1920,7 +1952,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { for (unsigned int j=0; j < n_actual_dofs; ++j) { - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; bool add_this = false; @@ -1937,56 +1969,34 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { Assert (loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); - if (dof_mask_is_active == true) - { - if (dof_mask[loc_row][loc_col] == true) - goto add_this_index; - } - else + if (dof_mask[loc_row][loc_col] == true) goto add_this_index; } // account for indirect contributions by // constraints - for (unsigned int p=0; p= 0 && loc_col < n_local_dofs, ExcInternalError()); - if (dof_mask_is_active == true) - { - if (dof_mask[my_indices[i].second.constraints[q].first][loc_col] == true) - goto add_this_index; - } - else + if (dof_mask[my_indices[i].constraints[q].first][loc_col] == true) goto add_this_index; } - for (unsigned int p=0; p &local_dof_indices, if (add_this == true) { add_this_index: - cols[col_counter++] = my_indices[j].first; + cols[col_counter++] = my_indices[j].global_row; } } } @@ -2043,7 +2053,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, ExcDimensionMismatch(dof_mask.n_cols(), n_local_dofs)); } - std::vector > my_indices (n_local_dofs); + std::vector my_indices (n_local_dofs); std::vector > constraint_lines; constraint_lines.reserve(n_local_dofs); @@ -2063,8 +2073,8 @@ add_entries_local_to_global (const std::vector &local_dof_indices, if (constraint_line_exists.size() <= local_dof_indices[i] || constraint_line_exists[local_dof_indices[i]] == false) { - my_indices[added_rows].first = local_dof_indices[i]; - my_indices[added_rows].second.local_row = i; + my_indices[added_rows].global_row = local_dof_indices[i]; + my_indices[added_rows].local_row = i; ++added_rows; continue; } @@ -2144,7 +2154,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, std::vector localized_indices (n_actual_dofs); for (unsigned int i=0; i &local_dof_indices, localized_indices[i] = sparsity_pattern.get_row_indices(). global_to_local(localized_indices[i]).second; } + + // easy case when we add all the indices + // unconditionally + if (dof_mask_is_active == false) + { + unsigned int block_start = 0; + Threads::ThreadMutex::ScopedLock lock(mutex); + for (unsigned int block=0; block::iterator index_it = localized_indices.begin(); + for (unsigned int block_col = 0; block_col actual_block_ends (num_blocks); std::vector cols (n_actual_dofs); @@ -2180,8 +2222,8 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { unsigned int col_counter = 0; std::vector::iterator block_it = block_ends.begin(); - const unsigned int row = my_indices[i].first; - const unsigned int loc_row = my_indices[i].second.local_row; + const unsigned int row = my_indices[i].global_row; + const unsigned int loc_row = my_indices[i].local_row; if (have_indirect_rows == false) { @@ -2201,16 +2243,12 @@ add_entries_local_to_global (const std::vector &local_dof_indices, ExcInternalError()); } - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; Assert(loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); bool add_this = true; - if (dof_mask_is_active == true) - if (dof_mask[loc_row][loc_col] == false) - add_this = false; - - if (add_this == true) + if (dof_mask[loc_row][loc_col] == true) cols[col_counter++] = localized_indices[j]; } @@ -2244,7 +2282,7 @@ add_entries_local_to_global (const std::vector &local_dof_indices, ExcInternalError()); } - const unsigned int loc_col = my_indices[j].second.local_row; + const unsigned int loc_col = my_indices[j].local_row; bool add_this = false; @@ -2257,50 +2295,28 @@ add_entries_local_to_global (const std::vector &local_dof_indices, { Assert (loc_col >= 0 && loc_col < n_local_dofs, ExcInternalError()); - if (dof_mask_is_active == true) - { - if (dof_mask[loc_row][loc_col] == true) - goto add_this_index; - } - else + if (dof_mask[loc_row][loc_col] == true) goto add_this_index; } - for (unsigned int p=0; p= 0 && loc_col < n_local_dofs, ExcInternalError()); - if (dof_mask_is_active == true) - { - if (dof_mask[my_indices[i].second.constraints[q].first][loc_col] == true) - goto add_this_index; - } - else + if (dof_mask[my_indices[i].constraints[q].first][loc_col] == true) goto add_this_index; } - for (unsigned int p=0; p &local_dof_indices, // now comes the hack that sets the // correct end of block in case we work // with block matrices. - while (n_actual_dofs == *block_it && block_it != block_ends.end()) + while (block_it != block_ends.end() && n_actual_dofs == *block_it) { actual_block_ends[block_it-block_ends.begin()] = col_counter; ++block_it;