From 6a27a9d28ba69ce410818ae6e0213e050aa78398 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Mon, 23 Feb 2009 17:18:49 +0000 Subject: [PATCH] Now the ConstraintMatrix seems to work also for inhomogeneous constraints, at least when using the distribute_local_to_global function. Still need to check the correctness of the other functions, which I will do some other time. git-svn-id: https://svn.dealii.org/trunk@18417 0785d39b-7218-0410-832d-ea1e28bc413d --- .../include/dofs/dof_constraints.templates.h | 64 +++++++---- .../deal.II/source/dofs/dof_constraints.cc | 104 ++++++++---------- 2 files changed, 91 insertions(+), 77 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_constraints.templates.h b/deal.II/deal.II/include/dofs/dof_constraints.templates.h index 5dccbde9d9..a8913a81f7 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -458,14 +458,21 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed, // diagonal, zero otherwise { for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p) - for (unsigned int q=0; - q!=lines[distribute[column]].entries.size(); ++q) - uncondensed.add (lines[distribute[row]].entries[p].first, - lines[distribute[column]].entries[q].first, - entry->value() * - lines[distribute[row]].entries[p].second * - lines[distribute[column]].entries[q].second); - + { + for (unsigned int q=0; + q!=lines[distribute[column]].entries.size(); ++q) + uncondensed.add (lines[distribute[row]].entries[p].first, + lines[distribute[column]].entries[q].first, + entry->value() * + lines[distribute[row]].entries[p].second * + lines[distribute[column]].entries[q].second); + + if (use_vectors == true) + vec(lines[distribute[row]].entries[p].first) -= + entry->value() * lines[distribute[row]].entries[p].second * + lines[distribute[row]].inhomogeneity; + } + // set old entry to correct // value entry->value() = (row == column ? average_diagonal : 0. ); @@ -657,12 +664,19 @@ ConstraintMatrix::condense (BlockSparseMatrix &uncondensed, const double old_value = entry->value (); for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p) - for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q) - uncondensed.add (lines[distribute[row]].entries[p].first, - lines[distribute[global_col]].entries[q].first, - old_value * - lines[distribute[row]].entries[p].second * - lines[distribute[global_col]].entries[q].second); + { + for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q) + uncondensed.add (lines[distribute[row]].entries[p].first, + lines[distribute[global_col]].entries[q].first, + old_value * + lines[distribute[row]].entries[p].second * + lines[distribute[global_col]].entries[q].second); + + if (use_vectors == true) + vec(lines[distribute[row]].entries[p].first) -= + old_value * lines[distribute[row]].entries[p].second * + lines[distribute[global_col]].inhomogeneity; + } entry->value() = (row == global_col ? average_diagonal : 0. ); } @@ -980,12 +994,19 @@ distribute_local_to_global (const FullMatrix &local_matrix, // the values again in // the matrix directly for (unsigned int p=0; pentries.size(); ++p) - for (unsigned int q=0; qentries.size(); ++q) - global_matrix.add (position_i->entries[p].first, - position_j->entries[q].first, - local_matrix(i,j) * - position_i->entries[p].second * - position_j->entries[q].second); + { + for (unsigned int q=0; qentries.size(); ++q) + global_matrix.add (position_i->entries[p].first, + position_j->entries[q].first, + local_matrix(i,j) * + position_i->entries[p].second * + position_j->entries[q].second); + + if (use_vectors == true) + global_vector (position_i->entries[p].first) -= + local_matrix(j,i) * position_i->entries[p].second * + position_j->inhomogeneity; + } // to make sure that the // global matrix remains @@ -1031,7 +1052,8 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (i == j) { column_indices[col_counter] = local_dof_indices[j]; - column_values[col_counter] = local_matrix(i,j); + column_values[col_counter] = + std::max(std::fabs(local_matrix(i,j)), 1e-10); col_counter++; } } diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index aea1a7ae1f..823d11a4bc 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -289,50 +289,41 @@ void ConstraintMatrix::close () } else // the DoF that we - // encountered is - // not constrained - // by a linear - // combination of - // other dofs but - // is equal to zero - // (i.e. its chain - // of entries is - // empty). in that - // case, we can't - // just overwrite - // the current - // entry, but we + // encountered is not + // constrained by a + // linear combination of + // other dofs but is + // equal to zero + // (i.e. its chain of + // entries is empty). in + // that case, we can't + // just overwrite the + // current entry, but we // have to actually // eliminate it { line->entries.erase (line->entries.begin()+entry); } - - //line->inhomogeneity += constrained_line->inhomogeneity * - // weight; - //constrained_line->inhomogeneity = 0; + line->inhomogeneity += constrained_line->inhomogeneity * + weight; - // now that we're - // here, do not - // increase index by - // one but rather - // make another pass - // for the present - // entry because we + // now that we're here, do + // not increase index by + // one but rather make + // another pass for the + // present entry because we // have replaced the - // present entry by - // another one, or - // because we have - // deleted it and - // shifted all - // following ones one + // present entry by another + // one, or because we have + // deleted it and shifted + // all following ones one // forward } else // entry not further - // constrained. just - // move ahead by one + // constrained. just move + // ahead by one ++entry; } @@ -342,41 +333,38 @@ void ConstraintMatrix::close () if (chained_constraint_replaced == false) break; - // increase iteration - // count. note that we should - // not iterate more times than - // there are constraints, since - // this puts a natural upper - // bound on the length of - // constraint chains + // increase iteration count. note + // that we should not iterate more + // times than there are constraints, + // since this puts a natural upper + // bound on the length of constraint + // chains ++iteration; Assert (iteration <= lines.size(), ExcInternalError()); } - // finally sort the entries and - // re-scale them if necessary. in - // this step, we also throw out - // duplicates as mentioned above + // finally sort the entries and re-scale + // them if necessary. in this step, we + // also throw out duplicates as mentioned + // above for (std::vector::iterator line = lines.begin(); line!=lines.end(); ++line) { std::sort (line->entries.begin(), line->entries.end()); - // loop over the now sorted - // list and see whether any of - // the entries references the - // same dofs more than once + // loop over the now sorted list and + // see whether any of the entries + // references the same dofs more than + // once for (unsigned int i=1; ientries.size(); ++i) if (line->entries[i].first == line->entries[i-1].first) { // ok, we've found a - // duplicate. go on to - // count how many - // duplicates there are - // so that we can - // allocate the right - // amount of memory + // duplicate. go on to count + // how many duplicates there + // are so that we can allocate + // the right amount of memory unsigned int duplicates = 1; for (unsigned int j=i+1; jentries.size(); ++j) if (line->entries[j].first == line->entries[j-1].first) @@ -461,8 +449,11 @@ void ConstraintMatrix::close () for (unsigned int i=0; ientries.size(); ++i) sum += line->entries[i].second; if ((sum != 1.0) && (std::fabs (sum-1.) < 1.e-13)) - for (unsigned int i=0; ientries.size(); ++i) - line->entries[i].second /= sum; + { + for (unsigned int i=0; ientries.size(); ++i) + line->entries[i].second /= sum; + line->inhomogeneity /= sum; + } } #ifdef DEBUG @@ -692,8 +683,9 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) j!=tmp_other_lines[i]->entries.end(); ++j) tmp.push_back (std::make_pair(j->first, j->second*weight)); - //line->inhomogeneity += tmp_other_lines[i]->inhomogeneity * weight; }; + line->inhomogeneity += tmp_other_lines[i]->inhomogeneity * + line->entries[i].second; }; // finally exchange old and // newly resolved line -- 2.39.5