]> https://gitweb.dealii.org/ - dealii.git/commitdiff
solved row vs column issue
authorJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Tue, 8 Mar 2016 09:09:35 +0000 (10:09 +0100)
committerJoscha Gedicke <joscha.gedicke@iwr.uni-heidelberg.de>
Tue, 8 Mar 2016 10:38:14 +0000 (11:38 +0100)
include/deal.II/lac/constraint_matrix.h
include/deal.II/lac/constraint_matrix.templates.h
include/deal.II/meshworker/simple.h
source/lac/constraint_matrix.cc
source/lac/constraint_matrix.inst.in

index b40514fd4c6c0d397795061bb6d5331a6714fd83..4fe203499472d098888cac58cfa4dcfaca273382 100644 (file)
@@ -744,11 +744,18 @@ public:
    * Same as the previous function, except that it uses two (possibly) different
    * index sets to correctly handle inhomogeneties when the local matrix is
    * computed from a combination of two neighboring elements, for example for an
-   * edge integral term in dG. Note that in the case that these two elements have
+   * edge integral term in DG. Note that in the case that these two elements have
    * different polynomial degree, the local matrix is rectangular.
    *
    * <tt>local_dof_indices_row<tt> is the set of row indices and
    * <tt>local_dof_indices_col<tt> is the set of column indices of the local matrix.
+   * <tt>diagonal=false<tt> says whether the two index sets are equal or not.
+   *
+   * If both index sets are equal, <tt>diagonal<tt> must be set to true or we
+   * simply use the previous function. If both index sets are different (diagonal=false)
+   * the <tt>global_vector<tt> is modified to handle inhomogeneties but no
+   * entries from <tt>local_vector<tt> are added. Note that the edge integrals for inner
+   * edged for DG do not contribute any values to the right hand side.
    */
   template <typename VectorType, typename LocalType>
   void
@@ -756,7 +763,8 @@ public:
                               const std::vector<size_type> &local_dof_indices_row,
                               const std::vector<size_type> &local_dof_indices_col,
                               VectorType                   &global_vector,
-                              const FullMatrix<LocalType>  &local_matrix) const;
+                              const FullMatrix<LocalType>  &local_matrix,
+                              bool diagonal = false) const;
 
   /**
    * Enter a single value into a result vector, obeying constraints.
index a87fffb1158db71290f741166be9588c29ae6a57..b990132b5b0499bed850300d4c6c083d8749ebe2 100644 (file)
@@ -623,7 +623,7 @@ distribute_local_to_global (const Vector<LocalType>      &local_vector,
                             VectorType                   &global_vector,
                             const FullMatrix<LocalType>  &local_matrix) const
 {
-  distribute_local_to_global(local_vector,local_dof_indices,local_dof_indices, global_vector, local_matrix);
+  distribute_local_to_global(local_vector,local_dof_indices,local_dof_indices, global_vector, local_matrix, true);
 }
 
 
@@ -635,16 +635,29 @@ distribute_local_to_global (const Vector<LocalType>      &local_vector,
                             const std::vector<size_type> &local_dof_indices_row,
                             const std::vector<size_type> &local_dof_indices_col,
                             VectorType                   &global_vector,
-                            const FullMatrix<LocalType>  &local_matrix) const
+                            const FullMatrix<LocalType>  &local_matrix,
+                            bool diagonal) const
 {
   Assert (sorted == true, ExcMatrixNotClosed());
-  AssertDimension (local_vector.size(), local_dof_indices_col.size());
+  AssertDimension (local_vector.size(), local_dof_indices_row.size());
   AssertDimension (local_matrix.m(), local_dof_indices_row.size());
   AssertDimension (local_matrix.n(), local_dof_indices_col.size());
-
-  const size_type n_local_dofs = local_vector.size();
+  
+  // diagonal checks if we have only one index set (if both are equal
+  // diagonal should be set to true).
+  // If true we do both, assemply of the right hand side (next lines)
+  // and (see further below) modifications of the right hand side
+  // according to the inhomogeneous constraints.
+  // Otherwise we only modify the right hand side according to
+  // local_matrix and the inhomogeneos constraints, and omit the vector add.
+
+  const size_type m_local_dofs = local_dof_indices_row.size();
+  const size_type n_local_dofs = local_dof_indices_col.size();
   if (lines.empty())
-    global_vector.add(local_dof_indices_col, local_vector);
+    {
+      if (diagonal)
+        global_vector.add(local_dof_indices_row, local_vector);
+    }
   else
     for (size_type i=0; i<n_local_dofs; ++i)
       {
@@ -654,7 +667,8 @@ distribute_local_to_global (const Vector<LocalType>      &local_vector,
         // the constraint
         if (is_constrained(local_dof_indices_col[i]) == false)
           {
-            global_vector(local_dof_indices_col[i]) += local_vector(i);
+            if (diagonal)
+              global_vector(local_dof_indices_row[i]) += local_vector(i);
             continue;
           }
 
@@ -669,7 +683,7 @@ distribute_local_to_global (const Vector<LocalType>      &local_vector,
         // constrained. If so, distribute the constraint
         const double val = position->inhomogeneity;
         if (val != 0)
-          for (size_type j=0; j<n_local_dofs; ++j)
+          for (size_type j=0; j<m_local_dofs; ++j)
             {
               if (is_constrained(local_dof_indices_row[j]) == false)
                 {
@@ -699,14 +713,17 @@ distribute_local_to_global (const Vector<LocalType>      &local_vector,
         // now distribute the constraint,
         // but make sure we don't touch
         // the entries of fixed dofs
-        for (size_type j=0; j<position->entries.size(); ++j)
+        if (diagonal)
           {
-            Assert (!(!local_lines.size()
+            for (size_type j=0; j<position->entries.size(); ++j)
+              {
+                Assert (!(!local_lines.size()
                       || local_lines.is_element(position->entries[j].first))
-                    || is_constrained(position->entries[j].first) == false,
-                    ExcMessage ("Tried to distribute to a fixed dof."));
-            global_vector(position->entries[j].first)
-            += local_vector(i) * position->entries[j].second;
+                      || is_constrained(position->entries[j].first) == false,
+                      ExcMessage ("Tried to distribute to a fixed dof."));
+                global_vector(position->entries[j].first)
+                += local_vector(i) * position->entries[j].second;
+              }
           }
       }
 }
index a08e9b00bea29fde06229c9a6b97c25e6503073c..b04373c43d18bb3d6194f5f991c626fc135f1ca2 100644 (file)
@@ -1198,9 +1198,6 @@ namespace MeshWorker
 
       if (ResidualSimple<VectorType>::constraints == 0)
         {
-          for (unsigned int i=0; i<i1.size(); ++i)
-            (*v)(i1[i]) += vector(i);
-
           for (unsigned int j=0; j<i1.size(); ++j)
             for (unsigned int k=0; k<i2.size(); ++k)
               if (std::fabs(M(j,k)) >= MatrixSimple<MatrixType>::threshold)
@@ -1208,7 +1205,7 @@ namespace MeshWorker
         }
       else
         {
-          ResidualSimple<VectorType>::constraints->distribute_local_to_global(vector, i1, i2, *v, M);
+          ResidualSimple<VectorType>::constraints->distribute_local_to_global(vector, i1, i2, *v, M, false);
           ResidualSimple<VectorType>::constraints->distribute_local_to_global(M, i1, i2, *MatrixSimple<MatrixType>::matrix[index]);
         }
     }
@@ -1238,11 +1235,11 @@ namespace MeshWorker
 
                 if (row == column)
                   assemble(info.matrix(k+m*n*n,false).matrix,
-                           info.vector(m).block(column), m,
-                           info.indices_by_block[column]);
+                           info.vector(m).block(row), m,
+                           info.indices_by_block[row]);
                 else
                   assemble(info.matrix(k+m*n*n,false).matrix,
-                           info.vector(m).block(column), m,
+                           info.vector(m).block(row), m,
                            info.indices_by_block[row],
                            info.indices_by_block[column]);
               }
@@ -1283,21 +1280,21 @@ namespace MeshWorker
 
                 if (row == column)
                   {
-                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
-                             info1.indices_by_block[column]);
-                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
-                             info2.indices_by_block[column]);
+                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(row), m,
+                             info1.indices_by_block[row]);
+                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(row), m,
+                             info2.indices_by_block[row]);
                   }
                 else
                   {
-                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(column), m,
+                    assemble(info1.matrix(k+m*n*n,false).matrix, info1.vector(m).block(row), m,
                              info1.indices_by_block[row], info1.indices_by_block[column]);
-                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(column), m,
+                    assemble(info2.matrix(k+m*n*n,false).matrix, info2.vector(m).block(row), m,
                              info2.indices_by_block[row], info2.indices_by_block[column]);
                   }
-                assemble(info1.matrix(k+m*n*n,true).matrix, info1.vector(m).block(column), m,
+                assemble(info1.matrix(k+m*n*n,true).matrix, info1.vector(m).block(row), m,
                          info1.indices_by_block[row], info2.indices_by_block[column]);
-                assemble(info2.matrix(k+m*n*n,true).matrix, info2.vector(m).block(column), m,
+                assemble(info2.matrix(k+m*n*n,true).matrix, info2.vector(m).block(row), m,
                          info2.indices_by_block[row], info1.indices_by_block[column]);
               }
         }
index 45247f81924cc564eddda76f492fb4eae83de43a..0391d54e33403e0f88c299742f9b5e333d163731 100644 (file)
@@ -1231,7 +1231,8 @@ ConstraintMatrix::resolve_indices (std::vector<types::global_dof_index> &indices
                                            const std::vector<ConstraintMatrix::size_type>  &, \
                                            const std::vector<ConstraintMatrix::size_type>  &, \
                                            VectorType                      &, \
-                                           const FullMatrix<VectorType::value_type>        &) const
+                                           const FullMatrix<VectorType::value_type> &, \
+                                           bool) const
 
 #define PARALLEL_VECTOR_FUNCTIONS(VectorType) \
   template void ConstraintMatrix:: \
@@ -1244,8 +1245,8 @@ ConstraintMatrix::resolve_indices (std::vector<types::global_dof_index> &indices
                                            const std::vector<ConstraintMatrix::size_type>  &, \
                                            const std::vector<ConstraintMatrix::size_type>  &, \
                                            VectorType                      &, \
-                                           const FullMatrix<VectorType::value_type>        &) const
-
+                                           const FullMatrix<VectorType::value_type> &, \
+                                           bool) const
 
 #ifdef DEAL_II_WITH_PETSC
 VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector);
index 8e3e3486be2dfeeca2243d106ca0b15ad4dcb2e0..5010a6f75c68b71807f998e0926fd9940fcc7990 100644 (file)
@@ -20,7 +20,7 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES)
     template void ConstraintMatrix::distribute_local_to_global<T<S> > (
       const Vector<S>&, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&) const;  
     template void ConstraintMatrix::distribute_local_to_global<T<S> > (
-      const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&) const;
+      const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, T<S> &, const FullMatrix<S>&, bool) const;
     template void ConstraintMatrix::set_zero<T<S> >(T<S> &) const;
   }
 
@@ -32,7 +32,7 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES)
     template void ConstraintMatrix::distribute_local_to_global<parallel::distributed::T<S> > (
       const Vector<S>&, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&) const;
     template void ConstraintMatrix::distribute_local_to_global<parallel::distributed::T<S> > (
-      const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&) const;
+      const Vector<S>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, parallel::distributed::T<S> &, const FullMatrix<S>&, bool) const;
     template void ConstraintMatrix::set_zero<parallel::distributed::T<S> >(parallel::distributed::T<S> &) const;
   }
 
@@ -44,7 +44,7 @@ for (V: EXTERNAL_SEQUENTIAL_VECTORS)
     template void ConstraintMatrix::distribute_local_to_global<V > (
       const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&) const;
     template void ConstraintMatrix::distribute_local_to_global<V > (
-      const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&) const;
+      const Vector<V::value_type>&, const std::vector<types::global_dof_index> &, const std::vector<types::global_dof_index> &, V&, const FullMatrix<V::value_type>&, bool) const;
     template void ConstraintMatrix::set_zero<V >(V&) const;
   }
 

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.