* @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);
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:
*/
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)
+inline
+bool
+ConstraintMatrix::ConstraintLine::operator == (const ConstraintLine &a) const
+{
+ return line == a.line;
+};
+
+
+
unsigned int
ConstraintMatrix::ConstraintLine::memory_consumption () const
{
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<ConstraintLine>::const_iterator
+ line=other_constraints.lines.begin();
+ line!=other_constraints.lines.end(); ++line)
+ for (std::vector<std::pair<unsigned int,double> >::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
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<std::pair<unsigned int,double> > tmp;
+ std::vector<std::vector<ConstraintLine>::const_iterator> tmp_other_lines;
+ for (std::vector<ConstraintLine>::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; i<line->entries.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<ConstraintLine>::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; i<line->entries.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<std::pair<unsigned int, double> >::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());
{
if (sorted == true)
{
-
ConstraintLine index_comparison;
index_comparison.line = index;
projects / dealii-svn.git / commitdiff