dolfin team mailing list archive

[Joachim Berdal Haga <jobh@xxxxxxxxx>]

And another one. I realized DirichletBC::apply is not symmetric, so
this small modification to zero_columns makes it possible to
symmetricise the diagonal blocks too.
I know this is not the preferred way to get a symmetric matrix, but
making it possible to do it on the LA level opens up some nice
possibilities for block matrices in python.
-j.

commit 0ae486596a7067305e1c4740eeacf45ab7202d64
Author: Joachim B Haga <jobh@xxxxxxxxx>
Date:   Thu Jan 20 13:33:52 2011 +0100

    Modify DirichletBC::zero_columns so that it optionally leaves the diagonal entries alone.

diff --git a/dolfin/fem/DirichletBC.cpp b/dolfin/fem/DirichletBC.cpp
index 75b7b17..9da5358 100644
--- a/dolfin/fem/DirichletBC.cpp
+++ b/dolfin/fem/DirichletBC.cpp
@@ -209,7 +209,7 @@ void DirichletBC::zero(GenericMatrix& A) const
   A.apply("insert");
 }
 //-----------------------------------------------------------------------------
-void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b) const
+void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b, bool keep_diagonal) const
 {
   Map bv_map;
   get_boundary_values(bv_map, _method);
@@ -241,17 +241,19 @@ void DirichletBC::zero_columns(GenericMatrix& A, GenericVector& b) const
     for (uint j=0; j<cols.size(); j++)
     {
       const uint col = cols[j];
-      if (is_dof[col] && vals[j] != 0.0)
+      if (keep_diagonal && col==row)
+        continue;
+      if (!is_dof[col] || vals[j] == 0.0)
+        continue;
+
+      if (!row_changed)
       {
-        if (!row_changed)
-        {
-          row_changed = true;
-          b_rows.push_back(row);
-          b_vals.push_back(0.0);
-        }
-        b_vals.back() -= dof_val[col]*vals[j];
-        vals[j] = 0.0;
+        row_changed = true;
+        b_rows.push_back(row);
+        b_vals.push_back(0.0);
       }
+      b_vals.back() -= dof_val[col]*vals[j];
+      vals[j] = 0.0;
     }
     if (row_changed) {
       A.setrow(row, cols, vals);
diff --git a/dolfin/fem/DirichletBC.h b/dolfin/fem/DirichletBC.h
index cf15127..2a10d10 100644
--- a/dolfin/fem/DirichletBC.h
+++ b/dolfin/fem/DirichletBC.h
@@ -171,7 +171,7 @@ namespace dolfin
 
     /// Make columns associated with boundary conditions zero, and
     /// update the RHS to reflect the changes. Useful for non-diagonals.
-    void zero_columns(GenericMatrix& A, GenericVector& b) const;
+    void zero_columns(GenericMatrix& A, GenericVector& b, bool keep_diagonal=false) const;
 
     /// Return boundary markers (facets stored as pairs of cells and local
     /// facet numbers)