]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement local_apply_boundary_values.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Mar 2004 20:27:37 +0000 (20:27 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 16 Mar 2004 20:27:37 +0000 (20:27 +0000)
git-svn-id: https://svn.dealii.org/trunk@8792 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/matrices.all_dimensions.cc

index e898ed0ab89716eff370c88753460210007d7688..5b218d01d60128648d924b1c85fc36c84ed764d0 100644 (file)
@@ -714,10 +714,35 @@ class MatrixTools : public MatrixCreator
                                      * the classes that are used to wrap
                                      * PETSc objects.
                                      *
-                                     * For a replacement function,
-                                     * see the documentation of the
-                                     * @ref{FilteredMatrix} class in
-                                     * the @p{LAC} sublibrary.
+                                     * Note that this function is not very
+                                     * efficient: it needs to alternatingly
+                                     * read and write into the matrix, a
+                                     * situation that PETSc does not handle
+                                     * too well. In addition, we only get rid
+                                     * of rows corresponding to boundary
+                                     * nodes, but the corresponding case of
+                                     * deleting the respective columns
+                                     * (i.e. if @arg eliminate_columns is @p
+                                     * true) is not presently implemented,
+                                     * and probably will never because it is
+                                     * too expensive without direct access to
+                                     * the PETSc data structures. A third
+                                     * reason against this function is that
+                                     * it doesn't handle the case where the
+                                     * matrix is distributed across an MPI
+                                     * system.
+                                     *
+                                     * In order to still be able to eliminate
+                                     * boundary values, it is better to get
+                                     * rid of them before the local matrices
+                                     * and vectors are distributed to the
+                                     * global ones, because then we don't
+                                     * have to mess with the sparse data
+                                     * structures. The
+                                     * local_apply_boundary_values() function
+                                     * does that, and is recommended for use
+                                     * instead of the global one for PETSc
+                                     * matrices and vectors.
                                      */
 #ifdef DEAL_II_USE_PETSC
     static void
@@ -727,6 +752,36 @@ class MatrixTools : public MatrixCreator
                           PETScWrappers::VectorBase  &right_hand_side,
                           const bool             eliminate_columns = true);
 #endif
+
+                                     /**
+                                      * Rather than applying boundary values
+                                      * to the global matrix and vector after
+                                      * assembly, this function does so @em
+                                      * before assembling, by modifying the
+                                      * local matrix and vector
+                                      * contributions. If you call this
+                                      * function on all local contributions,
+                                      * the resulting matrix will have the
+                                      * same entries, and the final call to
+                                      * apply_boundary_values() on the global
+                                      * system will not be necessary.
+                                      *
+                                      * Since this function does not have to
+                                      * work on the complicated data
+                                      * structures of sparse matrices, it is
+                                      * relatively cheap. It may therefore be
+                                      * a win if you have many fixed degrees
+                                      * of freedom (e.g. boundary nodes), or
+                                      * if access to the sparse matrix is
+                                      * expensive (e.g. for block sparse
+                                      * matrices, or for PETSc matrices).
+                                      */
+    static void
+    local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                                 const std::vector<unsigned int> &local_dof_indices,
+                                 FullMatrix<double> &local_matrix,
+                                 Vector<double>     &local_rhs,
+                                 const bool          eliminate_columns);
     
                                     /**
                                      * Exception
index 7834575749d7b4075cc1f600867fc5a98eacad2a..3301c70663de4095385a2263deccad9c03560d71 100644 (file)
@@ -13,6 +13,7 @@
 
 
 #include <numerics/matrices.h>
+#include <lac/full_matrix.h>
 #include <lac/vector.h>
 #include <lac/block_vector.h>
 #include <lac/sparse_matrix.h>
@@ -642,6 +643,74 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
 #endif
 
 
+void
+MatrixTools::
+local_apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                             const std::vector<unsigned int> &local_dof_indices,
+                             FullMatrix<double> &local_matrix,
+                             Vector<double>     &local_rhs,
+                             const bool          eliminate_columns)
+{
+  Assert (local_dof_indices.size() == local_matrix.m(),
+          ExcDimensionMismatch(local_dof_indices.size(),
+                               local_matrix.m()));
+  Assert (local_dof_indices.size() == local_matrix.n(),
+          ExcDimensionMismatch(local_dof_indices.size(),
+                               local_matrix.n()));
+  Assert (local_dof_indices.size() == local_rhs.size(),
+          ExcDimensionMismatch(local_dof_indices.size(),
+                               local_rhs.size()));
+
+                                   // if there is nothing to do, then exit
+                                   // right away
+  if (boundary_values.size() == 0)
+    return;
+
+                                   // otherwise traverse all the dofs used in
+                                   // the local matrices and vectors and see
+                                   // what's there to do
+  const unsigned int n_local_dofs = local_dof_indices.size();
+  for (unsigned int i=0; i<n_local_dofs; ++i)
+    {
+      const std::map<unsigned int, double>::const_iterator
+        boundary_value = boundary_values.find (local_dof_indices[i]);
+      if (boundary_value != boundary_values.end())
+        {
+                                           // remove this row, except for the
+                                           // diagonal element
+          for (unsigned j=0; j<n_local_dofs; ++j)
+            if (i != j)
+              local_matrix(i,j) = 0;
+
+                                           // replace diagonal entry by its
+                                           // absolute value to make sure that
+                                           // everything remains positive, or
+                                           // by one if zero
+          if (local_matrix(i,i) == 0.)
+            local_matrix(i,i) = 1.;
+          else
+            local_matrix(i,i) = std::fabs(local_matrix(i,i));
+          
+                                           // and replace rhs entry by correct
+                                           // value
+          local_rhs(i) = local_matrix(i,i) * boundary_value->second;
+
+                                           // finally do the elimination step
+                                           // if requested
+          if (eliminate_columns == true)
+            {
+              for (unsigned int row=0; row<n_local_dofs; ++row)
+                if (row != i)
+                  {
+                    local_rhs(row) -= local_matrix(row,i) * boundary_value->second;
+                    local_matrix(row,i) = 0;
+                  }
+            }
+        }
+    }
+}
+
+
 
 
 

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.