]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Now the ConstraintMatrix seems to work also for inhomogeneous constraints, at least...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Feb 2009 17:18:49 +0000 (17:18 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 23 Feb 2009 17:18:49 +0000 (17:18 +0000)
git-svn-id: https://svn.dealii.org/trunk@18417 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/dofs/dof_constraints.templates.h
deal.II/deal.II/source/dofs/dof_constraints.cc

index 5dccbde9d9f83ecd96e20fe7b3ba39eea0d0d4eb..a8913a81f7a3080b6d4b69c2f5b7377a6d5688ae 100644 (file)
@@ -458,14 +458,21 @@ ConstraintMatrix::condense (SparseMatrix<number> &uncondensed,
                                                  // diagonal, zero otherwise
                 {
                   for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
-                    for (unsigned int q=0;
-                         q!=lines[distribute[column]].entries.size(); ++q)
-                      uncondensed.add (lines[distribute[row]].entries[p].first,
-                                       lines[distribute[column]].entries[q].first,
-                                       entry->value() *
-                                       lines[distribute[row]].entries[p].second *
-                                       lines[distribute[column]].entries[q].second);
-               
+                   {
+                     for (unsigned int q=0;
+                          q!=lines[distribute[column]].entries.size(); ++q)
+                       uncondensed.add (lines[distribute[row]].entries[p].first,
+                                        lines[distribute[column]].entries[q].first,
+                                        entry->value() *
+                                        lines[distribute[row]].entries[p].second *
+                                        lines[distribute[column]].entries[q].second);
+
+                     if (use_vectors == true)
+                       vec(lines[distribute[row]].entries[p].first) -= 
+                         entry->value() * lines[distribute[row]].entries[p].second *
+                         lines[distribute[row]].inhomogeneity;
+                   }
+
                                                    // set old entry to correct
                                                    // value
                   entry->value() = (row == column ? average_diagonal : 0. );
@@ -657,12 +664,19 @@ ConstraintMatrix::condense (BlockSparseMatrix<number> &uncondensed,
                       const double old_value = entry->value ();
                          
                       for (unsigned int p=0; p!=lines[distribute[row]].entries.size(); ++p)
-                        for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
-                          uncondensed.add (lines[distribute[row]].entries[p].first,
-                                           lines[distribute[global_col]].entries[q].first,
-                                           old_value *
-                                           lines[distribute[row]].entries[p].second *
-                                           lines[distribute[global_col]].entries[q].second);
+                       {
+                         for (unsigned int q=0; q!=lines[distribute[global_col]].entries.size(); ++q)
+                           uncondensed.add (lines[distribute[row]].entries[p].first,
+                                            lines[distribute[global_col]].entries[q].first,
+                                            old_value *
+                                            lines[distribute[row]].entries[p].second *
+                                            lines[distribute[global_col]].entries[q].second);
+
+                         if (use_vectors == true)
+                           vec(lines[distribute[row]].entries[p].first) -= 
+                             old_value * lines[distribute[row]].entries[p].second *
+                             lines[distribute[global_col]].inhomogeneity;
+                       }
 
                       entry->value() = (row == global_col ? average_diagonal : 0. );
                     }
@@ -980,12 +994,19 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                                                    // the values again in
                                                    // the matrix directly
                   for (unsigned int p=0; p<position_i->entries.size(); ++p)
-                    for (unsigned int q=0; q<position_j->entries.size(); ++q)
-                      global_matrix.add (position_i->entries[p].first,
-                                         position_j->entries[q].first,
-                                         local_matrix(i,j) *
-                                         position_i->entries[p].second *
-                                         position_j->entries[q].second);
+                   {
+                     for (unsigned int q=0; q<position_j->entries.size(); ++q)
+                       global_matrix.add (position_i->entries[p].first,
+                                          position_j->entries[q].first,
+                                          local_matrix(i,j) *
+                                          position_i->entries[p].second *
+                                          position_j->entries[q].second);
+
+                     if (use_vectors == true)
+                       global_vector (position_i->entries[p].first) -=
+                         local_matrix(j,i) * position_i->entries[p].second *
+                         position_j->inhomogeneity;
+                   }
 
                                                    // to make sure that the
                                                    // global matrix remains
@@ -1031,7 +1052,8 @@ distribute_local_to_global (const FullMatrix<double>        &local_matrix,
                   if (i == j)
                    {
                      column_indices[col_counter] = local_dof_indices[j];
-                     column_values[col_counter] = local_matrix(i,j);
+                     column_values[col_counter] = 
+                       std::max(std::fabs(local_matrix(i,j)), 1e-10);
                      col_counter++;
                    }
                 }
index aea1a7ae1f04cb8522568f9a564b889529fe10c9..823d11a4bc3ed2da7b06068603c5f4f5b6758758 100644 (file)
@@ -289,50 +289,41 @@ void ConstraintMatrix::close ()
                  }
                else
                                                   // the DoF that we
-                                                  // encountered is
-                                                  // not constrained
-                                                  // by a linear
-                                                  // combination of
-                                                  // other dofs but
-                                                  // is equal to zero
-                                                  // (i.e. its chain
-                                                  // of entries is
-                                                  // empty). in that
-                                                  // case, we can't
-                                                  // just overwrite
-                                                  // the current
-                                                  // entry, but we
+                                                  // encountered is not
+                                                  // constrained by a
+                                                  // linear combination of
+                                                  // other dofs but is
+                                                  // equal to zero
+                                                  // (i.e. its chain of
+                                                  // entries is empty). in
+                                                  // that case, we can't
+                                                  // just overwrite the
+                                                  // current entry, but we
                                                   // have to actually
                                                   // eliminate it
                  {
                    line->entries.erase (line->entries.begin()+entry);
                  }
-                   
 
-               //line->inhomogeneity += constrained_line->inhomogeneity *
-               //                     weight;
-               //constrained_line->inhomogeneity = 0;
+               line->inhomogeneity += constrained_line->inhomogeneity *
+                                      weight;
                
-                                                // now that we're
-                                                // here, do not
-                                                // increase index by
-                                                // one but rather
-                                                // make another pass
-                                                // for the present
-                                                // entry because we
+                                                // now that we're here, do
+                                                // not increase index by
+                                                // one but rather make
+                                                // another pass for the
+                                                // present entry because we
                                                 // have replaced the
-                                                // present entry by
-                                                // another one, or
-                                                // because we have
-                                                // deleted it and
-                                                // shifted all
-                                                // following ones one
+                                                // present entry by another
+                                                // one, or because we have
+                                                // deleted it and shifted
+                                                // all following ones one
                                                 // forward
              }
            else
                                               // entry not further
-                                              // constrained. just
-                                              // move ahead by one
+                                              // constrained. just move
+                                              // ahead by one
              ++entry;
        }
 
@@ -342,41 +333,38 @@ void ConstraintMatrix::close ()
       if (chained_constraint_replaced == false)
        break;
 
-                                      // increase iteration
-                                      // count. note that we should
-                                      // not iterate more times than
-                                      // there are constraints, since
-                                      // this puts a natural upper
-                                      // bound on the length of
-                                      // constraint chains
+                                      // increase iteration count. note
+                                      // that we should not iterate more
+                                      // times than there are constraints,
+                                      // since this puts a natural upper
+                                      // bound on the length of constraint
+                                      // chains
       ++iteration;
       Assert (iteration <= lines.size(),
              ExcInternalError());
     }
 
-                                  // finally sort the entries and
-                                  // re-scale them if necessary. in
-                                  // this step, we also throw out
-                                  // duplicates as mentioned above
+                                  // finally sort the entries and re-scale
+                                  // them if necessary. in this step, we
+                                  // also throw out duplicates as mentioned
+                                  // above
   for (std::vector<ConstraintLine>::iterator line = lines.begin();
        line!=lines.end(); ++line)
     {
       std::sort (line->entries.begin(), line->entries.end());
 
-                                       // loop over the now sorted
-                                       // list and see whether any of
-                                       // the entries references the
-                                       // same dofs more than once
+                                       // loop over the now sorted list and
+                                       // see whether any of the entries
+                                       // references the same dofs more than
+                                       // once
       for (unsigned int i=1; i<line->entries.size(); ++i)
         if (line->entries[i].first == line->entries[i-1].first)
           {
                                              // ok, we've found a
-                                             // duplicate. go on to
-                                             // count how many
-                                             // duplicates there are
-                                             // so that we can
-                                             // allocate the right
-                                             // amount of memory
+                                             // duplicate. go on to count
+                                             // how many duplicates there
+                                             // are so that we can allocate
+                                             // the right amount of memory
             unsigned int duplicates = 1;
             for (unsigned int j=i+1; j<line->entries.size(); ++j)
               if (line->entries[j].first == line->entries[j-1].first)
@@ -461,8 +449,11 @@ void ConstraintMatrix::close ()
       for (unsigned int i=0; i<line->entries.size(); ++i)
        sum += line->entries[i].second;
       if ((sum != 1.0) && (std::fabs (sum-1.) < 1.e-13))
-       for (unsigned int i=0; i<line->entries.size(); ++i)
-         line->entries[i].second /= sum;
+       {
+         for (unsigned int i=0; i<line->entries.size(); ++i)
+           line->entries[i].second /= sum;
+         line->inhomogeneity /= sum;
+       }
     }
   
 #ifdef DEBUG
@@ -692,8 +683,9 @@ void ConstraintMatrix::merge (const ConstraintMatrix &other_constraints)
                   j!=tmp_other_lines[i]->entries.end(); ++j)
                tmp.push_back (std::make_pair(j->first, j->second*weight));
 
-             //line->inhomogeneity += tmp_other_lines[i]->inhomogeneity * weight;
            };
+         line->inhomogeneity += tmp_other_lines[i]->inhomogeneity *  
+           line->entries[i].second;
        };
                                       // finally exchange old and
                                       // newly resolved line

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.