]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Extend ConstraintMatrix::merge.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Jul 2001 21:18:02 +0000 (21:18 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 30 Jul 2001 21:18:02 +0000 (21:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@4864 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.h
deal.II/deal.II/source/dofs/dof_constraints.cc

index cd293b1a04445eabca06ad142d6872a39e6d3e6b..cc1c2ab6834b990afc1f6475009ca1929d8741d8 100644 (file)
@@ -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)
index dc827038d91bcd804f10fcdefb5aa1b680e06232..6271db259fe5b5fe567882f624bc7c660445d639 100644 (file)
@@ -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<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
@@ -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<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());
@@ -935,7 +1078,6 @@ bool ConstraintMatrix::is_constrained (const unsigned int index) const
 {
   if (sorted == true)
     {
-      
       ConstraintLine index_comparison;
       index_comparison.line = index;
       

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.