From e4c01806391e94450a065786487ef55caffc39f5 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 30 Jul 2001 21:18:02 +0000 Subject: [PATCH] Extend ConstraintMatrix::merge. git-svn-id: https://svn.dealii.org/trunk@4864 0785d39b-7218-0410-832d-ea1e28bc413d --- .../deal.II/include/dofs/dof_constraints.h | 94 ++++++++++- .../deal.II/source/dofs/dof_constraints.cc | 146 +++++++++++++++++- 2 files changed, 232 insertions(+), 8 deletions(-) diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index cd293b1a04..cc1c2ab683 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -189,18 +189,80 @@ class ConstraintMatrix : public Subscriptor * @p{close} called before), if * this object was closed before, * then it will be closed - * afterwards as well. + * afterwards as well. Note, + * however, that if the other + * argument is closed, then + * merging may be significantly + * faster. * * Note that the constraints in * each of the two objects (the * old one represented by this * object and the argument) may * not refer to the same degree - * of freedom, since this may - * result in incompatible - * constraints. If this is - * nevertheless the case, an - * exception is thrown. + * of freedom, i.e. a degree of + * freedom that is constrained in + * one object may not be + * constrained in the second. If + * this is nevertheless the case, + * an exception is thrown. + * + * However, the following is + * possible: if DoF @p{x} is + * constrained to dofs @p{x_i} + * for some set of indices @p{i}, + * then the DoFs @p{x_i} may be + * further constrained by the + * constraints object given as + * argument, although not to + * other DoFs that are + * constrained in either of the + * two objects. Note that it is + * not possible that the DoFs + * @p{x_i} are constrained within + * the present object. + * + * Because of simplicity of + * implementation, and also to + * avoid cycles, this operation + * is not symmetric: degrees of + * freedom that are constrained + * in the given argument object + * may not be constrained to DoFs + * that are themselves + * constrained within the present + * object. + * + * The aim of these merging + * operations is that if, for + * example, you have hanging + * nodes that are constrained to + * the degrees of freedom + * adjacent to them, you cannot + * originally, i.e. within one + * object, constrain these + * adjacent nodes + * further. However, that may be + * desirable in some cases, for + * example if they belong to a + * symmetry boundary for which + * the nodes on one side of the + * domain should have the same + * values as those on the other + * side. In that case, you would + * first construct a costraints + * object holding the hanging + * nodes constraints, and a + * second one that contains the + * constraints due to the + * symmetry boundary. You would + * then finally merge this second + * one into the first, possibly + * eliminating constraints of + * hanging nodes to adjacent + * boundary nodes by constraints + * to nodes at the opposite + * boundary. */ void merge (const ConstraintMatrix &other_constraints); @@ -498,6 +560,15 @@ class ConstraintMatrix : public Subscriptor int, << "Degree of freedom " << arg1 << " is constrained from both object in a merge operation."); + /** + * Exception + */ + DeclException1 (ExcDoFIsConstrainedToConstrainedDoF, + int, + << "In the given argument a degree of freedom is constrained " + << "to another DoF with number " << arg1 + << ", which however is constrained by this object. This is not" + << " allowed."); private: @@ -542,6 +613,17 @@ class ConstraintMatrix : public Subscriptor */ bool operator < (const ConstraintLine &) const; + /** + * This operator is likewise + * weird: it checks whether + * the line indices of the + * two operands are equal, + * irrespective of the fact + * that the contents of the + * line may be different. + */ + bool operator == (const ConstraintLine &) const; + /** * Determine an estimate for the * memory consumption (in bytes) diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index dc827038d9..6271db259f 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -47,6 +47,15 @@ ConstraintMatrix::ConstraintLine::operator < (const ConstraintLine &a) const +inline +bool +ConstraintMatrix::ConstraintLine::operator == (const ConstraintLine &a) const +{ + return line == a.line; +}; + + + unsigned int ConstraintMatrix::ConstraintLine::memory_consumption () const { @@ -260,6 +269,25 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) line!=other_constraints.lines.end(); ++line) AssertThrow (this_dofs.find (line->line) == this_dofs.end(), ExcDoFIsConstrainedFromBothObjects (line->line)); + + // finally check the following: + // while we allow that in this + // object nodes are constrained + // to other nodes that are + // constrained in the given + // argument, we do not allow + // the reverse, i.e. the nodes + // to which the constraints in + // the other object hold may + // not be constrained here + for (std::vector::const_iterator + line=other_constraints.lines.begin(); + line!=other_constraints.lines.end(); ++line) + for (std::vector >::const_iterator + e=line->entries.begin(); + e!=line->entries.end(); ++e) + AssertThrow (this_dofs.find (e->first) == this_dofs.end(), + ExcDoFIsConstrainedToConstrainedDoF (e->first)); }; // store the previous state with @@ -267,7 +295,122 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints) const bool object_was_sorted = sorted; sorted = false; - // append new lines at the end + + // first action is to fold into the + // present object possible + // constraints in the second + // object. for this, loop over all + // constraints and replace the + // constraint lines with a new one + // where constraints are replaced + // if necessary. use the same tmp + // object over again to avoid + // excessive memory allocation + std::vector > tmp; + std::vector::const_iterator> tmp_other_lines; + for (std::vector::iterator line=lines.begin(); + line!=lines.end(); ++line) + { + // copy the line of old object + // modulo dofs constrained in + // the second object. for this + // purpose, first search the + // respective constraint line + // (if any, otherwise a null + // pointer) in the other object + // for each of the entries in + // this line + // + // store whether we have to + // resolve entries, since if + // not there is no need to copy + // the line one-to-one + tmp.resize (0); + tmp_other_lines.resize (0); + tmp_other_lines.resize (line->entries.size()); + + bool entries_to_resolve = false; + + for (unsigned int i=0; ientries.size(); ++i) + { + if (other_constraints.sorted == true) + { + ConstraintLine index_comparison; + index_comparison.line = line->entries[i].first; +//TODO:[WB] we should use a binary search, since we are sure that the array is sorted + tmp_other_lines[i] = + std::find (other_constraints.lines.begin (), + other_constraints.lines.end (), + index_comparison); + } + else + { + tmp_other_lines[i] = other_constraints.lines.end (); + for (std::vector::const_iterator + p=other_constraints.lines.begin(); + p!=other_constraints.lines.end(); ++p) + if (p->line == line->entries[i].first) + { + tmp_other_lines[i] = p; + break; + }; + }; + + if (tmp_other_lines[i] != other_constraints.lines.end ()) + entries_to_resolve = true; + }; + + + // now we have for each entry + // in the present line whether + // it needs to be resolved + // using the new object, and if + // so which constraint line to + // use. first check whether we + // have to resolve anything at + // all, otherwise leave the + // line as is + if (entries_to_resolve == false) + continue; + + // something to resolve, so go + // about it + for (unsigned int i=0; ientries.size(); ++i) + { + // if the present dof is not + // constrained, then simply + // copy it over + if (tmp_other_lines[i] == other_constraints.lines.end()) + tmp.push_back(line->entries[i]); + else + // otherwise resolve + // further constraints by + // replacing the old + // entry by a sequence of + // new entries taken from + // the other object, but + // with multiplied + // weights + { + Assert (tmp_other_lines[i]->line == line->entries[i].first, + ExcInternalError()); + + const double weight = line->entries[i].second; + for (std::vector >::const_iterator + j=tmp_other_lines[i]->entries.begin(); + j!=tmp_other_lines[i]->entries.end(); ++j) + tmp.push_back (make_pair(j->first, j->second*weight)); + }; + }; + // finally exchange old and + // newly resolved line + line->entries.swap (tmp); + }; + + + + // next action: append new lines at + // the end lines.insert (lines.end(), other_constraints.lines.begin(), other_constraints.lines.end()); @@ -935,7 +1078,6 @@ bool ConstraintMatrix::is_constrained (const unsigned int index) const { if (sorted == true) { - ConstraintLine index_comparison; index_comparison.line = index; -- 2.39.5