simplify logic

git-svn-id: https://svn.dealii.org/trunk@32166 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/include/deal.II/meshworker/simple.h b/deal.II/include/deal.II/meshworker/simple.h
index ec8241b64fb93b1ef9dd27a2f312b94162786ef1..224729da485d29a235ddda1cf043a6efe7b5c08c 100644 (file)
--- a/deal.II/include/deal.II/meshworker/simple.h
+++ b/deal.II/include/deal.II/meshworker/simple.h
@@ -824,7 +824,7 @@ namespace MeshWorker
     {
       AssertDimension(M.m(), i1.size());
       AssertDimension(M.n(), i2.size());
-      //Assert(mg_constrained_dofs == 0, ExcInternalError());
+      Assert(mg_constrained_dofs == 0, ExcInternalError());
 //TODO: Possibly remove this function all together
 
       for (unsigned int j=0; j<i1.size(); ++j)
@@ -857,28 +857,35 @@ namespace MeshWorker
         {
           for (unsigned int j=0; j<i1.size(); ++j)
             for (unsigned int k=0; k<i2.size(); ++k)
-              if (std::fabs(M(j,k)) >= threshold)
-                if (!mg_constrained_dofs->at_refinement_edge(level, i1[j]) &&
-                    !mg_constrained_dofs->at_refinement_edge(level, i2[k]))
-                  {
-                    if (mg_constrained_dofs->set_boundary_values())
-                      {
-                        // At the boundary, only enter the term on the
-                        // diagonal, but not the coupling terms
-                        if ((!mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                             !mg_constrained_dofs->is_boundary_index(level, i2[k]))
-                            ||
-                            (mg_constrained_dofs->is_boundary_index(level, i1[j]) &&
-                             mg_constrained_dofs->is_boundary_index(level, i2[k]) &&
-                             i1[j] == i2[k]))
-                          G.add(i1[j], i2[k], M(j,k));
-                      }
-                    else
-                      G.add(i1[j], i2[k], M(j,k));
-                  }
-        }
+             {
+               // Only enter the local values into the global matrix,
+               //  if the value is larger than the threshold
+               if (std::fabs(M(j,k)) < threshold)
+                 continue;
+               
+               // Do not enter, if either the row or the column
+               // corresponds to an index on the refinement edge. The
+               // level problems are solved with homogeneous
+               // Dirichlet boundary conditions, therefore we
+               // eliminate these rows and columns. The corresponding
+               // matrix entries are entered by assemble_in() and
+               // assemble_out().
+                if (mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
+                    mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+                 continue;
+               
+               // At the boundary, only enter the term on the
+               // diagonal, but not the coupling terms
+               if ((mg_constrained_dofs->is_boundary_index(level, i1[j]) ||
+                    mg_constrained_dofs->is_boundary_index(level, i2[k])) &&
+                   (i1[j] != i2[k]))
+                 continue;
+               
+               G.add(i1[j], i2[k], M(j,k));
+             }
+       }
     }
-
+    
 
     template <class MATRIX>
     inline void
@@ -1101,20 +1108,10 @@ namespace MeshWorker
 
             if (level1 == level2)
               {
-                if (mg_constrained_dofs == 0)
-                  {
-                    assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
-                  }
-                else
-                  {
-                    assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
-                    assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
-                  }
+               assemble((*matrix)[level1], info1.matrix(k,false).matrix, info1.indices_by_block[row], info1.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info1.matrix(k,true).matrix, info1.indices_by_block[row], info2.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info2.matrix(k,false).matrix, info2.indices_by_block[row], info2.indices_by_block[column], level1);
+               assemble((*matrix)[level1], info2.matrix(k,true).matrix, info2.indices_by_block[row], info1.indices_by_block[column], level1);
               }
             else
               {