Dear all,
in my opinion there is a small bug in the implementation of the
ConstraintMatrix member function add_entries in file constraint_matrix.cc
(release Version 9.0, lines 224-261). The current piece of code is:
224 void
225 ConstraintMatrix::add_entries
226 (const size_type line,
227 const std::vector<std::pair<size_type,double> > &col_val_pairs)
228 {
229 Assert (sorted==false, ExcMatrixIsClosed());
230 Assert (is_constrained(line), ExcLineInexistant(line));
231
232 ConstraintLine *line_ptr = &lines[lines_cache[calculate_line_index(
line)]];
233 Assert (line_ptr->index == line, ExcInternalError());
234
235 // if in debug mode, check whether an entry for this column already
236 // exists and if its the same as the one entered at present
237 //
238 // in any case: skip this entry if an entry for this column already
239 // exists, since we don't want to enter it twice
240 for (std::vector<std::pair<size_type,double> >::const_iterator
241 col_val_pair = col_val_pairs.begin();
242 col_val_pair!=col_val_pairs.end(); ++col_val_pair)
243 {
244 Assert (line != col_val_pair->first,
245 ExcMessage ("Can't constrain a degree of freedom to
itself"));
246
247 for (ConstraintLine::Entries::const_iterator
248 p=line_ptr->entries.begin();
249 p != line_ptr->entries.end(); ++p)
250 if (p->first == col_val_pair->first)
251 {
252 // entry exists, break innermost loop
253 Assert (p->second == col_val_pair->second,
254 ExcEntryAlreadyExists(line, col_val_pair->first,
255 p->second, col_val_pair->
second));
256 break;
257 }
258
259 line_ptr->entries.push_back (*col_val_pair);
260 }
261 }
Although, it is stated in the comments (lines 235-239) there is no skip of
the entry col_val_pair in the case this entry already exist. The break
command terminates the innermost loop but the push_back operation is
performed nevertheless, i.e. for every element of col_val_pairs. This leads
to duplicates in the constraint matrix.
I solved the problem by just addind the three lines 246a, 255a, 258a.
224 void
225 ConstraintMatrix::add_entries
226 (const size_type line,
227 const std::vector<std::pair<size_type,double> > &col_val_pairs)
228 {
229 Assert (sorted==false, ExcMatrixIsClosed());
230 Assert (is_constrained(line), ExcLineInexistant(line));
231
232 ConstraintLine *line_ptr = &lines[lines_cache[calculate_line_index(
line)]];
233 Assert (line_ptr->index == line, ExcInternalError());
234
235 // if in debug mode, check whether an entry for this column already
236 // exists and if its the same as the one entered at present
237 //
238 // in any case: skip this entry if an entry for this column already
239 // exists, since we don't want to enter it twice
240 for (std::vector<std::pair<size_type,double> >::const_iterator
241 col_val_pair = col_val_pairs.begin();
242 col_val_pair!=col_val_pairs.end(); ++col_val_pair)
243 {
244 Assert (line != col_val_pair->first,
245 ExcMessage ("Can't constrain a degree of freedom to
itself"));
246
246a bool entry_exists = false;
247 for (ConstraintLine::Entries::const_iterator
248 p=line_ptr->entries.begin();
249 p != line_ptr->entries.end(); ++p)
250 if (p->first == col_val_pair->first)
251 {
252 // entry exists, break innermost loop
253 Assert (p->second == col_val_pair->second,
254 ExcEntryAlreadyExists(line, col_val_pair->first,
255 p->second,
col_val_pair->second));
255a entry_exists = true;
256 break;
257 }
258
258a if (entry_exists == fase)
259 line_ptr->entries.push_back (*col_val_pair);
260 }
261 }
May be there is a more elegant way to do this.
Best regards
Dustin
--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see
https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups
"deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
For more options, visit https://groups.google.com/d/optout.
