Add inhomogeneity to the ConstraintMatrix class. Still need to write some new functio...

git-svn-id: https://svn.dealii.org/trunk@18390 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h
index 176092f490e95a5188e8fb9762a368e086eaaf74..ea3ec36ac8b8dd742ca73a24126f94bc6adf80cd 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_constraints.h
+++ b/deal.II/deal.II/include/dofs/dof_constraints.h
@@ -270,7 +270,7 @@ class BlockIndices;
  * a second step.
  *
  * @ingroup dofs
- * @author Wolfgang Bangerth, 1998, 2004, 2008
+ * @author Wolfgang Bangerth, Martin Kronbichler, 1998, 2004, 2008, 2009
  */
 class ConstraintMatrix : public Subscriptor
 {
@@ -331,6 +331,15 @@ class ConstraintMatrix : public Subscriptor
     void add_entries (const unsigned int                                  line,
                      const std::vector<std::pair<unsigned int,double> > &col_val_pairs);
 
+                                    /**
+                                     * Set an imhomogeneity to the
+                                     * constraint line <i>i</i>,
+                                     * according to the discussion in
+                                     * the general class description.
+                                     */
+    void set_inhomogeneity (const unsigned int line,
+                           const double       value);
+
                                     /**
                                      * Close the filling of
                                      * entries. Since the lines of a
@@ -1073,7 +1082,7 @@ class ConstraintMatrix : public Subscriptor
 
                                     /**
                                      * Delete hanging nodes in a
-                                     * vector.  Sets all hanging node
+                                     * vector. Sets all hanging node
                                      * values to zero. The @p
                                      * VectorType may be a
                                      * Vector<float>, Vector<double>,
@@ -1230,6 +1239,11 @@ class ConstraintMatrix : public Subscriptor
                                          */
        std::vector<std::pair<unsigned int,double> > entries;
 
+                                        /**
+                                         * Value of the inhomogeneity.
+                                         */
+        double inhomogeneity;
+
                                         /**
                                          * This operator is a bit
                                          * weird and unintuitive: it
@@ -1422,6 +1436,7 @@ ConstraintMatrix::add_line (const unsigned int line)
                                   // list
   lines.push_back (ConstraintLine());
   lines.back().line = line;
+  lines.back().inhomogeneity = 0.;
 }
 
 
@@ -1482,6 +1497,38 @@ ConstraintMatrix::add_entry (const unsigned int line,
 
 
 
+inline
+void
+ConstraintMatrix::set_inhomogeneity (const unsigned int line,
+                                    const double       value)
+{
+  std::vector<ConstraintLine>::iterator line_ptr;
+  const std::vector<ConstraintLine>::const_iterator start=lines.begin();
+
+                                  // the usual case is that the line where
+                                  // a value is entered is the one we
+                                  // added last, so we search backward
+  for (line_ptr=(lines.end()-1); line_ptr!=start; --line_ptr)
+    if (line_ptr->line == line)
+      break;
+
+                                  // if the loop didn't break, then
+                                  // line_ptr must be begin().
+                                  // we have an error if that doesn't
+                                  // point to 'line' then
+  Assert (line_ptr->line==line, ExcLineInexistant(line));
+                                  // second check: the respective
+                                  // flag in the
+                                  // constraint_line_exists field
+                                  // must exist
+  Assert (line < constraint_line_exists.size(), ExcInternalError());
+  Assert (constraint_line_exists[line] == true, ExcInternalError());
+
+  line_ptr->inhomogeneity = value;
+}
+
+
+
 inline
 bool
 ConstraintMatrix::is_constrained (const unsigned int index) const 

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 befc30a16bd22f9abcfd370a2d6b0e937f07e6df..69c40bb38c2c050f85d2cd985e404df35a5bc7da 100644 (file)
--- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h
+++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h
@@ -1102,7 +1102,7 @@ ConstraintMatrix::distribute (const VectorType &condensed,
                                         // line was condensed away,
                                         // create it newly. first set
                                         // it to zero
-       uncondensed(line) = 0.;
+       uncondensed(line) = next_constraint->inhomogeneity;
                                         // then add the different
                                         // contributions
        for (unsigned int i=0; i<next_constraint->entries.size(); ++i)
@@ -1126,7 +1126,7 @@ ConstraintMatrix::distribute (VectorType &vec) const
                                       // fill entry in line
                                       // next_constraint.line by adding the
                                       // different contributions
-      double new_value = 0;
+      double new_value = next_constraint->inhomogeneity;
       for (unsigned int i=0; i<next_constraint->entries.size(); ++i)
        new_value += (vec(next_constraint->entries[i].first) *
                       next_constraint->entries[i].second);

diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc
index 80c251d5eec7ca9f6a8203790a12551f567ef146..846b9df307d230e4d02b4ca3c5f864f337b53de8 100644 (file)
--- a/deal.II/deal.II/source/dofs/dof_constraints.cc
+++ b/deal.II/deal.II/source/dofs/dof_constraints.cc
@@ -306,6 +306,9 @@ void ConstraintMatrix::close ()
                    line->entries.erase (line->entries.begin()+entry);
                  }
                    
+
+               line->inhomogeneity += constrained_line->inhomogeneity *
+                                      weight;
                
                                                 // now that we're
                                                 // here, do not
@@ -685,6 +688,8 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
                     j=tmp_other_lines[i]->entries.begin();
                   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;
            };
        };
                                       // finally exchange old and