Assert (vec.has_ghost_elements() == false,
ExcMessage ("This operation can only be performed on vectors without ghost elements."));
- typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
- ConstraintLine index_comparison;
- index_comparison.line = vec.local_range().first;
- const constraint_iterator begin_my_constraints =
- Utilities::lower_bound (lines.begin(),lines.end(),index_comparison);
-
- index_comparison.line = vec.local_range().second;
- const constraint_iterator end_my_constraints
- = Utilities::lower_bound(lines.begin(),lines.end(),index_comparison);
-
- // all indices we need to read from
- IndexSet my_indices (vec.size());
-
- const std::pair<unsigned int, unsigned int>
- local_range = vec.local_range();
+ IndexSet vec_owned_elements (vec.size());
+ vec_owned_elements.add_range(vec.local_range().first,
+ vec.local_range().second);
- my_indices.add_range (local_range.first, local_range.second);
+ typedef std::vector<ConstraintLine>::const_iterator constraint_iterator;
- std::set<unsigned int> individual_indices;
- for (constraint_iterator it = begin_my_constraints;
- it != end_my_constraints; ++it)
+#ifdef DEBUG
+ // the algorithm below pushes contributions c_ij x_j to a constrained
+ // DoF x_j *from the side of the owner of x_j*, as opposed to pulling
+ // it from the owner of the target side x_i. this relies on the
+ // assumption that both target and source know about the constraint
+ // on x_i.
+ //
+ // the disadvantage is that it is incredibly difficult to find out what
+ // is happening if the assumption is not satisfied. to help debug
+ // problems of this kind, we do the following in a first step if
+ // debugging is enabled:
+ // - every processor who owns an x_j where c_ij!=0 sends a one
+ // to the owner of x_i to add up
+ // - the owner of x_i knows how many nonzero entries exist, so can
+ // verify that the correct number of ones has been added
+ // - if the sum is not correct, then apparently one of the owners
+ // of the x_j's did not know to send its one and, consequently,
+ // will also not know to send the correct c_ij*x_j later on,
+ // leading to a bug
+ set_zero (vec);
+ for (constraint_iterator it = lines.begin();
+ it != lines.end(); ++it)
for (unsigned int i=0; i<it->entries.size(); ++i)
- if ((it->entries[i].first < local_range.first)
- ||
- (it->entries[i].first >= local_range.second))
- individual_indices.insert (it->entries[i].first);
-
- my_indices.add_indices (individual_indices.begin(),
- individual_indices.end());
-
- IndexSet local_range_is (vec.size());
- local_range_is.add_range(local_range.first, local_range.second);
-
+ if (vec_owned_elements.is_element (it->entries[i].first))
+ vec(it->line) += 1;
+ vec.compress (VectorOperation::add);
+
+ DeclException2 (ExcIncorrectConstraint, int, int,
+ << "While distributing the constraint for DoF "
+ << arg1 << ", it turns out that one of the processors "
+ << "who own the " << arg2
+ << " degrees of freedom that x_" << arg1
+ << " is constrained against does not know about "
+ << "the constraint on x_" << arg1
+ << ". Did you not initialize the ConstraintMatrix "
+ << "with the appropriate locally_relevant set so "
+ << "that every processor who owns a DoF that constrains "
+ << "another DoF also knows about this constraint?");
+ for (constraint_iterator it = lines.begin();
+ it != lines.end(); ++it)
+ if (vec_owned_elements.is_element(it->line))
+ Assert (vec(it->line) == it->entries.size(),
+ ExcIncorrectConstraint(it->line, it->entries.size()));
+#endif
- // create a vector and import those indices
- PETScWrappers::MPI::Vector ghost_vec (vec.get_mpi_communicator(),
- local_range_is,
- my_indices);
- ghost_vec = vec;
- ghost_vec.update_ghost_values();
+ // set the values of all constrained DoFs to zero, then add the
+ // contributions for each constrained DoF of all the ones we
+ // own locally
+ set_zero (vec);
- // finally do the distribution on own constraints
- for (constraint_iterator it = begin_my_constraints;
- it != end_my_constraints; ++it)
- {
- // fill entry in line next_constraint.line by adding the different
- // contributions
- PetscScalar new_value = it->inhomogeneity;
- for (unsigned int i=0; i<it->entries.size(); ++i)
- new_value += (PetscScalar(ghost_vec(it->entries[i].first)) *
- it->entries[i].second);
- vec(it->line) = new_value;
- }
+ for (constraint_iterator it = lines.begin();
+ it != lines.end(); ++it)
+ for (unsigned int i=0; i<it->entries.size(); ++i)
+ if (vec_owned_elements.is_element(it->entries[i].first))
+ vec(it->line) += it->entries[i].second *
+ vec(it->entries[i].first);
- vec.compress (VectorOperation::insert);
+ // now compress to communicate the entries that we added to
+ // and that weren't to local processors to the owner
+ vec.compress (VectorOperation::add);
}
//
// compared to the _01 test, here the ConstraintMatrix object acts on an index
// set that only includes the locally owned vector elements, without any
-// overlap. this verifies that we really only need to know about the *targets*
-// of constraints locally, not the *sources*. (at the time of writing this
-// test, ConstraintMatrix::distribute for PETSc distributed vectors actively
-// imported the source vector elements for those target elements we locally
-// own.)
+// overlap. this means that the owners of a source DoF in a constraint
+// do not know about the constraint if the constrained DoF is owned by
+// another processor. This is not allowed, and should trigger an
+// exception.
#include "../tests.h"
}
cm.close ();
- // now distribute these constraints
- cm.distribute (vec);
-
- // verify correctness
- vec.compress ();
-
- if (myid != 0)
- Assert (vec(vec.local_range().first+10) == vec.local_range().first-25,
- ExcInternalError());
-
- if (myid != n_processes-1)
- Assert (vec(vec.local_range().first+90) == vec.local_range().first+105,
- ExcInternalError());
-
- for (unsigned int i=vec.local_range().first; i<vec.local_range().second; ++i)
+ // now distribute these constraints. this should trigger an exception
+ bool exception_caught = false;
+ deal_II_exceptions::disable_abort_on_exception();
+ try
{
- if ((i != vec.local_range().first+10)
- &&
- (i != vec.local_range().first+90))
- {
- double val = vec(i);
- Assert (std::fabs(val - i) <= 1e-6, ExcInternalError());
- }
+ cm.distribute (vec);
}
+ catch (...)
+ {
+ exception_caught = true;
+ }
+ if (n_processes > 1)
+ Assert (exception_caught == true, ExcInternalError());
- {
- double exact_l1 = 0;
-
- // add up original values of vector entries
- for (unsigned int i=0; i<vec.size(); ++i)
- exact_l1 += i;
-
- // but then correct for the constrained values
- for (unsigned int p=0; p<n_processes; ++p)
- {
- if (p != 0)
- exact_l1 = exact_l1 - (p*100+10) + (p*100-25);
- if (p != n_processes-1)
- exact_l1 = exact_l1 - (p*100+90) + (p*100+105);
- }
-
- Assert (vec.l1_norm() == exact_l1, ExcInternalError());
-
- if (myid == 0)
- deallog << "Norm = " << vec.l1_norm() << std::endl;
- }
+ if (myid == 0)
+ deallog << "OK" << std::endl;
}
projects / dealii-svn.git / commitdiff