From 6fc5a74416ae093e7af4f29bac34fbc57a8ec5b9 Mon Sep 17 00:00:00 2001 From: Joscha Gedicke Date: Tue, 8 Mar 2016 10:09:35 +0100 Subject: [PATCH] solved row vs column issue --- include/deal.II/lac/constraint_matrix.h | 12 ++++- .../deal.II/lac/constraint_matrix.templates.h | 45 +++++++++++++------ include/deal.II/meshworker/simple.h | 27 +++++------ source/lac/constraint_matrix.cc | 7 +-- source/lac/constraint_matrix.inst.in | 6 +-- 5 files changed, 60 insertions(+), 37 deletions(-) diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h index b40514fd4c..4fe2034994 100644 --- a/include/deal.II/lac/constraint_matrix.h +++ b/include/deal.II/lac/constraint_matrix.h @@ -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. * * local_dof_indices_row is the set of row indices and * local_dof_indices_col is the set of column indices of the local matrix. + * diagonal=false says whether the two index sets are equal or not. + * + * If both index sets are equal, diagonal must be set to true or we + * simply use the previous function. If both index sets are different (diagonal=false) + * the global_vector is modified to handle inhomogeneties but no + * entries from local_vector are added. Note that the edge integrals for inner + * edged for DG do not contribute any values to the right hand side. */ template void @@ -756,7 +763,8 @@ public: const std::vector &local_dof_indices_row, const std::vector &local_dof_indices_col, VectorType &global_vector, - const FullMatrix &local_matrix) const; + const FullMatrix &local_matrix, + bool diagonal = false) const; /** * Enter a single value into a result vector, obeying constraints. diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h index a87fffb115..b990132b5b 100644 --- a/include/deal.II/lac/constraint_matrix.templates.h +++ b/include/deal.II/lac/constraint_matrix.templates.h @@ -623,7 +623,7 @@ distribute_local_to_global (const Vector &local_vector, VectorType &global_vector, const FullMatrix &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 &local_vector, const std::vector &local_dof_indices_row, const std::vector &local_dof_indices_col, VectorType &global_vector, - const FullMatrix &local_matrix) const + const FullMatrix &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 &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 &local_vector, // constrained. If so, distribute the constraint const double val = position->inhomogeneity; if (val != 0) - for (size_type j=0; j &local_vector, // now distribute the constraint, // but make sure we don't touch // the entries of fixed dofs - for (size_type j=0; jentries.size(); ++j) + if (diagonal) { - Assert (!(!local_lines.size() + for (size_type j=0; jentries.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; + } } } } diff --git a/include/deal.II/meshworker/simple.h b/include/deal.II/meshworker/simple.h index a08e9b00be..b04373c43d 100644 --- a/include/deal.II/meshworker/simple.h +++ b/include/deal.II/meshworker/simple.h @@ -1198,9 +1198,6 @@ namespace MeshWorker if (ResidualSimple::constraints == 0) { - for (unsigned int i=0; i= MatrixSimple::threshold) @@ -1208,7 +1205,7 @@ namespace MeshWorker } else { - ResidualSimple::constraints->distribute_local_to_global(vector, i1, i2, *v, M); + ResidualSimple::constraints->distribute_local_to_global(vector, i1, i2, *v, M, false); ResidualSimple::constraints->distribute_local_to_global(M, i1, i2, *MatrixSimple::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]); } } diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc index 45247f8192..0391d54e33 100644 --- a/source/lac/constraint_matrix.cc +++ b/source/lac/constraint_matrix.cc @@ -1231,7 +1231,8 @@ ConstraintMatrix::resolve_indices (std::vector &indices const std::vector &, \ const std::vector &, \ VectorType &, \ - const FullMatrix &) const + const FullMatrix &, \ + bool) const #define PARALLEL_VECTOR_FUNCTIONS(VectorType) \ template void ConstraintMatrix:: \ @@ -1244,8 +1245,8 @@ ConstraintMatrix::resolve_indices (std::vector &indices const std::vector &, \ const std::vector &, \ VectorType &, \ - const FullMatrix &) const - + const FullMatrix &, \ + bool) const #ifdef DEAL_II_WITH_PETSC VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector); diff --git a/source/lac/constraint_matrix.inst.in b/source/lac/constraint_matrix.inst.in index 8e3e3486be..5010a6f75c 100644 --- a/source/lac/constraint_matrix.inst.in +++ b/source/lac/constraint_matrix.inst.in @@ -20,7 +20,7 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES) template void ConstraintMatrix::distribute_local_to_global > ( const Vector&, const std::vector &, T &, const FullMatrix&) const; template void ConstraintMatrix::distribute_local_to_global > ( - const Vector&, const std::vector &, const std::vector &, T &, const FullMatrix&) const; + const Vector&, const std::vector &, const std::vector &, T &, const FullMatrix&, bool) const; template void ConstraintMatrix::set_zero >(T &) const; } @@ -32,7 +32,7 @@ for (S: REAL_SCALARS; T : DEAL_II_VEC_TEMPLATES) template void ConstraintMatrix::distribute_local_to_global > ( const Vector&, const std::vector &, parallel::distributed::T &, const FullMatrix&) const; template void ConstraintMatrix::distribute_local_to_global > ( - const Vector&, const std::vector &, const std::vector &, parallel::distributed::T &, const FullMatrix&) const; + const Vector&, const std::vector &, const std::vector &, parallel::distributed::T &, const FullMatrix&, bool) const; template void ConstraintMatrix::set_zero >(parallel::distributed::T &) const; } @@ -44,7 +44,7 @@ for (V: EXTERNAL_SEQUENTIAL_VECTORS) template void ConstraintMatrix::distribute_local_to_global ( const Vector&, const std::vector &, V&, const FullMatrix&) const; template void ConstraintMatrix::distribute_local_to_global ( - const Vector&, const std::vector &, const std::vector &, V&, const FullMatrix&) const; + const Vector&, const std::vector &, const std::vector &, V&, const FullMatrix&, bool) const; template void ConstraintMatrix::set_zero(V&) const; } -- 2.39.5