Small bug in ConstraintMatrix::add_entries () ?!

[Dustin Kumor]

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

[Wolfgang Bangerth]

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+un...@googlegroups.com
> <mailto:dealii+un...@googlegroups.com>.
> For more options, visit https://groups.google.com/d/optout.


--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Dustin Kumor]

Dear Wolfgang

I can put that into a pull request but I've never done this before so I'm going to need someone guiding me through the whole process.
Same for the test writing issue.

Best regards
Dustin

[Daniel Arndt]

Dustin,

I can put that into a pull request but I've never done this before so I'm going to need someone guiding me through the whole process.
Same for the test writing issue.

Have a look at https://github.com/dealii/dealii/wiki/Contributing and don't hesitate to ask if you have questions.

Best,

Daniel

[Dustin Kumor]

Daniel,

thanks for the hint. I will have a look at this on the weekend.

Best,
Dustin