]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Add apply_boundary_values for block matrices and vectors.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 May 2000 14:30:39 +0000 (14:30 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 17 May 2000 14:30:39 +0000 (14:30 +0000)
At present not implemented, but I check it in to work at it at home.

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

deal.II/deal.II/source/numerics/matrices.cc

index 4d26689880844606e8d4fba5a9f1b9fe6cc1de41..168e7f2ee73b9b60cacd227e5c9cab707080e99c 100644 (file)
@@ -675,6 +675,175 @@ MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundar
 };
 
 
+
+/*
+
+template <int dim>
+template <int blocks>
+void
+MatrixTools<dim>::apply_boundary_values (const map<unsigned int,double> &boundary_values,
+                                        BlockSparseMatrix<double,blocks,blocks>  &matrix,
+                                        BlockVector<blocks,double>   &solution,
+                                        BlockVector<blocks,double>   &right_hand_side,
+                                        const bool                    preserve_symmetry)
+{
+  Assert (matrix.n() == matrix.m(),
+         ExcDimensionsDontMatch(matrix.n(), matrix.m()));
+  Assert (matrix.n() == right_hand_side.size(),
+         ExcDimensionsDontMatch(matrix.n(), right_hand_side.size()));
+  Assert (matrix.n() == solution.size(),
+         ExcDimensionsDontMatch(matrix.n(), solution.size()));
+  Assert (matrix.get_row_indices() == matrix.get_column_indices(),
+         ExcMatrixNotBlockSquare());
+  
+                                  // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+  map<unsigned int,double>::const_iterator  dof  = boundary_values.begin(),
+                                           endd = boundary_values.end();
+  const unsigned int n_dofs             = matrix.m();
+
+                                  // if a diagonal entry is zero
+                                  // later, then we use another
+                                  // number instead. take it to be
+                                  // the first nonzero diagonal
+                                  // element of the matrix, or 1 if
+                                  // there is no such thing
+  double first_nonzero_diagonal_entry = 1;
+  for (unsigned int i=0; i<n_dofs; ++i)
+    if (matrix.diag_element(i) != 0)
+      {
+       first_nonzero_diagonal_entry = matrix.diag_element(i);
+       break;
+      };
+
+  
+  for (; dof != endd; ++dof)
+    {
+      Assert (dof->first < n_dofs, ExcInternalError());
+      
+      const unsigned int dof_number = dof->first;
+                                      // for each boundary dof:
+      
+                                      // set entries of this line
+                                      // to zero except for the diagonal
+                                      // entry. Note that the diagonal
+                                      // entry is always the first one
+                                      // for square matrices, i.e.
+                                      // we shall not set
+                                      // matrix.global_entry(
+                                      //     sparsity_rowstart[dof.first])
+      const unsigned int last = sparsity_rowstart[dof_number+1];
+      for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+       matrix.global_entry(j) = 0.;
+
+
+                                      // set right hand side to
+                                      // wanted value: if main diagonal
+                                      // entry nonzero, don't touch it
+                                      // and scale rhs accordingly. If
+                                      // zero, take the first main
+                                      // diagonal entry we can find, or
+                                      // one if no nonzero main diagonal
+                                      // element exists. Normally, however,
+                                      // the main diagonal entry should
+                                      // not be zero.
+                                      //
+                                      // store the new rhs entry to make
+                                      // the gauss step more efficient
+      double new_rhs;
+      if (matrix.diag_element(dof_number) != 0.0)
+       new_rhs = right_hand_side(dof_number)
+               = dof->second * matrix.diag_element(dof_number);
+      else
+       {
+         matrix.set (dof_number, dof_number,
+                     first_nonzero_diagonal_entry);
+         new_rhs = right_hand_side(dof_number)
+                 = dof->second * first_nonzero_diagonal_entry;
+       };
+
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row
+      if (preserve_symmetry)
+       {
+                                          // store the only nonzero entry
+                                          // of this line for the Gauss
+                                          // elimination step
+         const double diagonal_entry = matrix.diag_element(dof_number);
+         
+                                          // we have to loop over all
+                                          // rows of the matrix which
+                                          // have a nonzero entry in
+                                          // the column which we work
+                                          // in presently. if the
+                                          // sparsity pattern is
+                                          // symmetric, then we can
+                                          // get the positions of
+                                          // these rows cheaply by
+                                          // looking at the nonzero
+                                          // column numbers of the
+                                          // present row. we need not
+                                          // look at the first entry,
+                                          // since that is the
+                                          // diagonal element and
+                                          // thus the present row
+         for (unsigned int j=sparsity_rowstart[dof_number]+1; j<last; ++j)
+           {
+             const unsigned int row = sparsity_colnums[j];
+
+                                              // find the position of
+                                              // element
+                                              // (row,dof_number)
+             const unsigned int *
+               p = lower_bound(&sparsity_colnums[sparsity_rowstart[row]+1],
+                               &sparsity_colnums[sparsity_rowstart[row+1]],
+                               dof_number);
+
+                                              // check whether this line has
+                                              // an entry in the regarding column
+                                              // (check for ==dof_number and
+                                              // != next_row, since if
+                                              // row==dof_number-1, *p is a
+                                              // past-the-end pointer but points
+                                              // to dof_number anyway...)
+                                              //
+                                              // there should be such an entry!
+             Assert ((*p == dof_number) &&
+                     (p != &sparsity_colnums[sparsity_rowstart[row+1]]),
+                     ExcInternalError());
+
+             const unsigned int global_entry
+               = (p - &sparsity_colnums[sparsity_rowstart[0]]);
+             
+                                              // correct right hand side
+             right_hand_side(row) -= matrix.global_entry(global_entry) /
+                                     diagonal_entry * new_rhs;
+             
+                                              // set matrix entry to zero
+             matrix.global_entry(global_entry) = 0.;
+           };
+       };
+
+                                      // preset solution vector
+      solution(dof_number) = dof->second;
+    };
+};
+
+
+*/
+
+
+
 template <int dim>
 MassMatrix<dim>::MassMatrix (const Function<dim> * const rhs,
                             const Function<dim> * const a) :

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.