Dustin,
good observation. I think you have found a bug indeed. I've put this
here for now:
https://github.com/dealii/dealii/issues/7101
Your fix seems reasonable to me. Do you want to put that into a pull
request? Let us know if you need help with that (and with writing a test).
(For reference, post-9.0, the class is now called AffineConstraints.)
Best and thanks for pointing this out!
Wolfgang
On 08/22/2018 08:30 AM, Dustin Kumor wrote:
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:
|
224void
225ConstraintMatrix::add_entries
226(constsize_type line,
227conststd::vector<std::pair<size_type,double>>&col_val_pairs)
228{
229Assert(sorted==false,ExcMatrixIsClosed());
230Assert(is_constrained(line),ExcLineInexistant(line));
231
232ConstraintLine*line_ptr =&lines[lines_cache[calculate_line_index(line)]];
233Assert(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
240for(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{
244Assert(line !=col_val_pair->first,
245ExcMessage("Can't constrain a degree of freedom to itself"));
246
247for(ConstraintLine::Entries::const_iterator
248 p=line_ptr->entries.begin();
249 p !=line_ptr->entries.end();++p)
250if(p->first ==col_val_pair->first)
251{
252// entry exists, break innermost loop
253Assert(p->second ==col_val_pair->second,
254ExcEntryAlreadyExists(line,col_val_pair->first,
255
p->second,col_val_pair->second));
256break;
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_pairin the case this entry already exist. The
breakcommand terminates the innermost loop but the push_backoperation 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.
|
224void
225ConstraintMatrix::add_entries
226(constsize_type line,
227conststd::vector<std::pair<size_type,double>>&col_val_pairs)
228{
229Assert(sorted==false,ExcMatrixIsClosed());
230Assert(is_constrained(line),ExcLineInexistant(line));
231
232ConstraintLine*line_ptr =&lines[lines_cache[calculate_line_index(line)]];
233Assert(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
240for(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 {
244Assert(line != col_val_pair->first,
245ExcMessage("Can't constrain a degree of freedom to itself"));
246
246aboolentry_exists =false;
247for(ConstraintLine::Entries::const_iterator
248 p=line_ptr->entries.begin();
249 p !=line_ptr->entries.end();++p)
250if(p->first ==col_val_pair->first)
251{
252// entry exists, break innermost loop
253Assert(p->second == col_val_pair->second,
254ExcEntryAlreadyExists(line, col_val_pair->first,
255 p->second,
col_val_pair->second));
255a entry_exists =true;
256break;
257}
258
258aif(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 dealii+unsubscr...@googlegroups.com
<mailto:dealii+unsubscr...@googlegroups.com>.
For more options, visit https://groups.google.com/d/optout.
--
------------------------------------------------------------------------
Wolfgang Bangerth email: bange...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/
--
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 dealii+unsubscr...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.
