* VectorTools::interpolate_boundary_values() function that returns
* a std::map and use it as input for MatrixTools::apply_boundary_values()
* to set boundary nodes to their correct value.
+ *
+ * Either behavior can also be achieved by building two separate
+ * ConstraintMatrix objects and calling ConstraintMatrix::merge function with
+ * a particular second argument.
*/
<h3>lac</h3>
<ol>
+ <li><p>New: The ConstraintMatrix::merge function now takes a second
+ argument that indicates what should happen if the two objects to be
+ merged have constraints on the very same degree of freedom.
+ <br>
+ (WB 2010/09/18)
+ </p></li>
+
<li><p>Updated: Changes to the <code>SLEPcWrappers</code> to preserve
compatibility with SLEPc version 3.1. Main new updated features are
new solver classes: (i) Power; and (ii) Generalized Davidson.
class ConstraintMatrix : public Subscriptor
{
public:
+ /**
+ * An enum that describes what should
+ * happen if the two ConstraintMatrix
+ * objects involved in a call to the
+ * merge() function happen to have
+ * constraints on the same degrees of
+ * freedom.
+ */
+ enum MergeConflictBehavior
+ {
+ /**
+ * Throw an exception if the two
+ * objects concerned have
+ * conflicting constraints on the
+ * same degree of freedom.
+ */
+ no_conflicts_allowed,
+
+ /**
+ * In an operation
+ * <code>cm1.merge(cm2)</code>, if
+ * <code>cm1</code> and
+ * <code>cm2</code> have
+ * constraints on the same degree
+ * of freedom, take the one from
+ * <code>cm1</code>.
+ */
+ left_object_wins,
+
+ /**
+ * In an operation
+ * <code>cm1.merge(cm2)</code>, if
+ * <code>cm1</code> and
+ * <code>cm2</code> have
+ * constraints on the same degree
+ * of freedom, take the one from
+ * <code>cm2</code>.
+ */
+ right_object_wins
+ };
+
/**
* Constructor. The supplied IndexSet
* defines which indices might be
* other argument is closed, then
* merging may be significantly faster.
*
- * Note that the constraints in each of
+ * Using the default value of the second
+ * arguments, the constraints in each of
* the two objects (the old one
* represented by this object and the
* argument) may not refer to the same
* freedom that is constrained in one
* object may not be constrained in the
* second. If this is nevertheless the
- * case, an exception is thrown.
+ * case, an exception is thrown. However,
+ * this behavior can be changed by
+ * providing a different value for the
+ * second argument.
*/
- void merge (const ConstraintMatrix &other_constraints);
+ void merge (const ConstraintMatrix &other_constraints,
+ const MergeConflictBehavior merge_conflict_behavior = no_conflicts_allowed);
/**
* Shift all entries of this matrix
* constrained.
*/
const std::vector<std::pair<unsigned int,double> >*
- get_constraint_entries (unsigned int line) const;
+ get_constraint_entries (const unsigned int line) const;
/**
* Returns the value of the inhomogeneity
-void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
+void
+ConstraintMatrix::merge (const ConstraintMatrix &other_constraints,
+ const MergeConflictBehavior merge_conflict_behavior)
{
//TODO: this doesn't work with IndexSets yet. [TH]
AssertThrow(local_lines.size()==0, ExcNotImplemented());
AssertThrow(other_constraints.local_lines.size()==0, ExcNotImplemented());
- // first check whether the
- // constraints in the two objects
- // are for different degrees of
- // freedom
-#ifdef DEBUG
- if (true)
- {
- // first insert all dofs in
- // this object into a list...
- std::set<unsigned int> this_dofs;
- for (std::vector<ConstraintLine>::const_iterator line=lines.begin();
- line!=lines.end(); ++line)
- this_dofs.insert (line->line);
-
- // ...then check whether it
- // appears in the other object
- // as well (in which case we
- // would need to issue an
- // error). note that we have to
- // do this in a somewhat
- // complicated style since the
- // two objects may not be
- // sorted
- for (std::vector<ConstraintLine>::const_iterator
- line=other_constraints.lines.begin();
- 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));
- }
-#endif
-
// store the previous state with
// respect to sorting
const bool object_was_sorted = sorted;
lines_cache.resize(other_constraints.lines_cache.size(),
numbers::invalid_unsigned_int);
- // first action is to fold into the
- // present object possible
- // constraints in the second
- // object. for this, loop over all
+ // first action is to fold into the present
+ // object possible constraints in the
+ // second object. we don't strictly need to
+ // do this any more since the
+ // ConstraintMatrix has learned to deal
+ // with chains of constraints in the
+ // close() function, but we have
+ // traditionally done this and it's not
+ // overly hard to do.
+ //
+ // 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
+ // if necessary.
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.clear ();
- tmp_other_lines.clear ();
- tmp_other_lines.reserve (line->entries.size());
-
-
- bool entries_to_resolve = false;
-
- for (unsigned int i=0; i<line->entries.size(); ++i)
- {
- if (other_constraints.sorted == true)
- {
- // as the array is
- // sorted, use a
- // bindary find to
- // check for the
- // existence of the
- // element. if it does
- // not exist, then the
- // pointer may still
- // point into tha
- // array, but to an
- // element of which the
- // indices do not
- // match, so return the
- // end iterator
- ConstraintLine index_comparison;
- index_comparison.line = line->entries[i].first;
-
- std::vector<ConstraintLine>::const_iterator
- it = std::lower_bound (other_constraints.lines.begin (),
- other_constraints.lines.end (),
- index_comparison);
- if ((it != other_constraints.lines.end ()) &&
- (it->line != index_comparison.line))
- it = other_constraints.lines.end ();
-
- tmp_other_lines.push_back (it);
- }
- else
- {
- std::vector<ConstraintLine>::const_iterator
- it = 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)
- {
- it = p;
- break;
- }
-
- tmp_other_lines.push_back (it);
- };
-
- if (tmp_other_lines.back() != other_constraints.lines.end ())
- entries_to_resolve = true;
- };
- Assert (tmp_other_lines.size() == line->entries.size(),
- ExcInternalError());
-
- // 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())
+ // constrained, or if we won't take
+ // the constraint from the other
+ // object, then simply copy it over
+ if (!other_constraints.is_constrained(line->entries[i].first)
+ ||
+ ((merge_conflict_behavior != right_object_wins)
+ &&
+ other_constraints.is_constrained(line->entries[i].first)
+ &&
+ this->is_constrained(line->entries[i].first)))
tmp.push_back(line->entries[i]);
else
// otherwise resolve
// with multiplied
// weights
{
- Assert (tmp_other_lines[i]->line == line->entries[i].first,
+ const std::vector<std::pair<unsigned int,double> >*
+ other_line
+ = other_constraints.get_constraint_entries (line->entries[i].first);
+ Assert (other_line != 0,
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)
+ j=other_line->begin();
+ j!=other_line->end(); ++j)
tmp.push_back (std::make_pair(j->first, j->second*weight));
- line->inhomogeneity += tmp_other_lines[i]->inhomogeneity *
- line->entries[i].second;
- };
- };
+ line->inhomogeneity += other_constraints.get_inhomogeneity(line->entries[i].first) *
+ line->entries[i].second;
+ }
+ }
// 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());
+ // next action: append those lines at the
+ // end that we want to add
+ for (std::vector<ConstraintLine>::const_iterator
+ line=other_constraints.lines.begin();
+ line!=other_constraints.lines.end(); ++line)
+ if (!is_constrained(line->line))
+ lines.push_back (*line);
+ else
+ {
+ // the constrained dof we want to
+ // copy from the other object is also
+ // constrained here. let's see what
+ // we should do with that
+ switch (merge_conflict_behavior)
+ {
+ case no_conflicts_allowed:
+ AssertThrow (false,
+ ExcDoFIsConstrainedFromBothObjects (line->line));
+ break;
+
+ case left_object_wins:
+ // ignore this constraint
+ break;
+
+ case right_object_wins:
+ // we need to replace the
+ // existing constraint by
+ // the one from the other
+ // object
+ lines[lines_cache[calculate_line_index(line->line)]].entries
+ = line->entries;
+ lines[lines_cache[calculate_line_index(line->line)]].inhomogeneity
+ = line->inhomogeneity;
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ }
+
// if the object was sorted before,
// then make sure it is so
// afterwards as well. otherwise
--- /dev/null
+//---------------------------- constraints_merge_03.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2003, 2004, 2005, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- constraints_merge_03.cc ---------------------------
+
+
+// merge and print a bunch of ConstrainMatrices. test the case where we have
+// conflicting constraints and don't allow any conflicts
+
+#include "../tests.h"
+#include <lac/constraint_matrix.h>
+#include <base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+std::ofstream logfile("constraints_merge_03/output");
+
+
+void merge_check ()
+{
+ deallog << "Checking ConstraintMatrix::merge" << std::endl;
+
+ // check twice, once with closed
+ // objects, once with open ones
+ for (unsigned int run=0; run<2; ++run)
+ {
+ deallog << "Checking with " << (run == 0 ? "open" : "closed")
+ << " objects" << std::endl;
+
+ // check that the `merge' function
+ // works correctly
+ ConstraintMatrix c1, c2;
+
+ // enter simple line
+ c1.add_line (0);
+ c1.add_entry (0, 11, 1.);
+ c1.set_inhomogeneity (0, 42);
+
+ // fill second constraints
+ // object that has a conflict
+ c2.add_line (0);
+ c2.add_entry (0, 13, 2.);
+ c2.set_inhomogeneity (0, 142);
+ // in one of the two runs,
+ // close the objects
+ if (run == 1)
+ {
+ c1.close ();
+ c2.close ();
+ };
+
+ // now merge the two and print the
+ // results
+ try
+ {
+ c1.merge (c2, ConstraintMatrix::no_conflicts_allowed);
+ }
+ catch (...)
+ {
+ deallog << "OK" << std::endl;
+ }
+
+ c1.print (logfile);
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ merge_check ();
+}
+
--- /dev/null
+
+DEAL::Checking ConstraintMatrix::merge
+DEAL::Checking with open objects
+DEAL::OK
+ 0 11: 1.0
+ 0: 42.
+DEAL::Checking with closed objects
+DEAL::OK
+ 0 11: 1.0
+ 0: 42.
--- /dev/null
+//---------------------------- constraints_merge_04.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2003, 2004, 2005, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- constraints_merge_04.cc ---------------------------
+
+
+// merge and print a bunch of ConstrainMatrices. test the case where we have
+// conflicting constraints and the left object wins
+
+#include "../tests.h"
+#include <lac/constraint_matrix.h>
+#include <base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+std::ofstream logfile("constraints_merge_04/output");
+
+
+void merge_check ()
+{
+ deallog << "Checking ConstraintMatrix::merge" << std::endl;
+
+ // check twice, once with closed
+ // objects, once with open ones
+ for (unsigned int run=0; run<2; ++run)
+ {
+ deallog << "Checking with " << (run == 0 ? "open" : "closed")
+ << " objects" << std::endl;
+
+ // check that the `merge' function
+ // works correctly
+ ConstraintMatrix c1, c2;
+
+ // enter simple line
+ c1.add_line (0);
+ c1.add_entry (0, 11, 1.);
+ c1.set_inhomogeneity (0, 42);
+
+ // fill second constraints
+ // object that has a conflict
+ c2.add_line (0);
+ c2.add_entry (0, 13, 2.);
+ c2.set_inhomogeneity (0, 142);
+ // in one of the two runs,
+ // close the objects
+ if (run == 1)
+ {
+ c1.close ();
+ c2.close ();
+ };
+
+ // now merge the two and print the
+ // results
+ try
+ {
+ c1.merge (c2, ConstraintMatrix::left_object_wins);
+ }
+ catch (...)
+ {
+ Assert (false, ExcInternalError());
+ }
+
+ c1.print (logfile);
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ merge_check ();
+}
+
--- /dev/null
+
+DEAL::Checking ConstraintMatrix::merge
+DEAL::Checking with open objects
+ 0 11: 1.0
+ 0: 42.
+DEAL::Checking with closed objects
+ 0 11: 1.0
+ 0: 42.
--- /dev/null
+//---------------------------- constraints_merge_05.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2003, 2004, 2005, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- constraints_merge_05.cc ---------------------------
+
+
+// merge and print a bunch of ConstrainMatrices. test the case where we have
+// conflicting constraints and the right object wins
+
+#include "../tests.h"
+#include <lac/constraint_matrix.h>
+#include <base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+std::ofstream logfile("constraints_merge_05/output");
+
+
+void merge_check ()
+{
+ deallog << "Checking ConstraintMatrix::merge" << std::endl;
+
+ // check twice, once with closed
+ // objects, once with open ones
+ for (unsigned int run=0; run<2; ++run)
+ {
+ deallog << "Checking with " << (run == 0 ? "open" : "closed")
+ << " objects" << std::endl;
+
+ // check that the `merge' function
+ // works correctly
+ ConstraintMatrix c1, c2;
+
+ // enter simple line
+ c1.add_line (0);
+ c1.add_entry (0, 11, 1.);
+ c1.set_inhomogeneity (0, 42);
+
+ // fill second constraints
+ // object that has a conflict
+ c2.add_line (0);
+ c2.add_entry (0, 13, 2.);
+ c2.set_inhomogeneity (0, 142);
+ // in one of the two runs,
+ // close the objects
+ if (run == 1)
+ {
+ c1.close ();
+ c2.close ();
+ };
+
+ // now merge the two and print the
+ // results
+ try
+ {
+ c1.merge (c2, ConstraintMatrix::right_object_wins);
+ }
+ catch (...)
+ {
+ Assert (false, ExcInternalError());
+ }
+
+ c1.print (logfile);
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ merge_check ();
+}
+
--- /dev/null
+
+DEAL::Checking ConstraintMatrix::merge
+DEAL::Checking with open objects
+ 0 13: 2.0
+ 0: 1.4e+02
+DEAL::Checking with closed objects
+ 0 13: 2.0
+ 0: 1.4e+02
--- /dev/null
+//---------------------------- constraints_merge_06.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2003, 2004, 2005, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- constraints_merge_06.cc ---------------------------
+
+
+// merge and print a bunch of ConstrainMatrices. test the case where we have
+// conflicting constraints and the right object wins
+//
+// like _05, but the winning constraint introduces another chain
+
+#include "../tests.h"
+#include <lac/constraint_matrix.h>
+#include <base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+std::ofstream logfile("constraints_merge_06/output");
+
+
+void merge_check ()
+{
+ deallog << "Checking ConstraintMatrix::merge" << std::endl;
+
+ // check twice, once with closed
+ // objects, once with open ones
+ for (unsigned int run=0; run<2; ++run)
+ {
+ deallog << "Checking with " << (run == 0 ? "open" : "closed")
+ << " objects" << std::endl;
+
+ // check that the `merge' function
+ // works correctly
+ ConstraintMatrix c1, c2;
+
+ // enter simple line
+ c1.add_line (0);
+ c1.add_entry (0, 11, 1.);
+ c1.set_inhomogeneity (0, 42);
+
+ c1.add_line (13);
+ c1.add_entry (13, 2, 0.5);
+ c1.set_inhomogeneity (13, 2);
+
+ // fill second constraints
+ // object that has a conflict
+ c2.add_line (0);
+ c2.add_entry (0, 13, 2.);
+ c2.set_inhomogeneity (0, 142);
+ // in one of the two runs,
+ // close the objects
+ if (run == 1)
+ {
+ c1.close ();
+ c2.close ();
+ };
+
+ // now merge the two and print the
+ // results
+ try
+ {
+ c1.merge (c2, ConstraintMatrix::right_object_wins);
+ }
+ catch (...)
+ {
+ Assert (false, ExcInternalError());
+ }
+
+ c1.print (logfile);
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ merge_check ();
+}
+
--- /dev/null
+
+DEAL::Checking ConstraintMatrix::merge
+DEAL::Checking with open objects
+ 0 13: 2.0
+ 0: 1.4e+02
+ 13 2: 0.50
+ 13: 2.0
+DEAL::Checking with closed objects
+ 0 2: 1.0
+ 0: 1.5e+02
+ 13 2: 0.50
+ 13: 2.0
--- /dev/null
+//---------------------------- constraints_merge_07.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 1998, 1999, 2000, 2001, 2003, 2004, 2005, 2008, 2010 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- constraints_merge_07.cc ---------------------------
+
+
+// merge and print a bunch of ConstrainMatrices. test the case where we have
+// conflicting constraints and the right object wins
+//
+// like _06, but the winning constraint introduces now introduces a cycle in
+// the graph of constraints
+
+#include "../tests.h"
+#include <lac/constraint_matrix.h>
+#include <base/logstream.h>
+
+#include <fstream>
+#include <iomanip>
+
+
+std::ofstream logfile("constraints_merge_07/output");
+
+
+void merge_check ()
+{
+ deallog << "Checking ConstraintMatrix::merge" << std::endl;
+
+ // check twice, once with closed
+ // objects, once with open ones
+ for (unsigned int run=0; run<2; ++run)
+ {
+ deallog << "Checking with " << (run == 0 ? "open" : "closed")
+ << " objects" << std::endl;
+
+ // check that the `merge' function
+ // works correctly
+ ConstraintMatrix c1, c2;
+
+ // enter simple line
+ c1.add_line (0);
+ c1.add_entry (0, 11, 1.);
+ c1.set_inhomogeneity (0, 42);
+
+ c1.add_line (13);
+ c1.add_entry (13, 2, 0.5);
+ c1.set_inhomogeneity (13, 2);
+
+ c1.add_line (2);
+ c1.add_entry (2, 0, 0.25);
+ c1.set_inhomogeneity (2, 22);
+
+ // fill second constraints
+ // object that has a conflict
+ c2.add_line (0);
+ c2.add_entry (0, 13, 2.);
+ c2.set_inhomogeneity (0, 142);
+ // in one of the two runs,
+ // close the objects
+ if (run == 1)
+ {
+ c1.close ();
+ c2.close ();
+ };
+
+ // now merge the two and print the
+ // results
+ try
+ {
+ c1.merge (c2, ConstraintMatrix::right_object_wins);
+ }
+ catch (...)
+ {
+ Assert (false, ExcInternalError());
+ }
+
+ c1.print (logfile);
+ }
+}
+
+
+int main ()
+{
+ deallog << std::setprecision (2);
+ logfile << std::setprecision (2);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ merge_check ();
+}
+
--- /dev/null
+
+DEAL::Checking ConstraintMatrix::merge
+DEAL::Checking with open objects
+ 0 13: 2.0
+ 0: 1.4e+02
+ 13 2: 0.50
+ 13: 2.0
+ 2 13: 0.50
+ 2: 58.
+DEAL::Checking with closed objects
+ 0 11: 0.25
+ 0: 1.8e+02
+ 2 11: 0.25
+ 2: 32.
+ 13 11: 0.12
+ 13: 18.