]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
some more modifications to use dg and cg in multigrid at the same time
authorjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 27 Sep 2010 10:13:50 +0000 (10:13 +0000)
committerjanssen <janssen@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 27 Sep 2010 10:13:50 +0000 (10:13 +0000)
git-svn-id: https://svn.dealii.org/trunk@22169 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_constrained_dofs.h
deal.II/deal.II/include/numerics/mesh_worker_assembler.h
deal.II/deal.II/source/multigrid/mg_tools.cc

index 42ec2f510f83bccdeb82d3a5d13ef13f11f3a1a2..b9d82d14b80d741923a3773906bb4e2d99adf6ff 100644 (file)
@@ -125,6 +125,7 @@ class MGConstrainedDoFs : public Subscriptor
                                      */
 
     bool set_boundary_values () const;
+    bool continuity_across_edges () const;
   private:
 
                                     /**
@@ -278,6 +279,20 @@ MGConstrainedDoFs::set_boundary_values () const
   return boundary_values_need_to_be_set;
 }
 
+inline
+bool 
+MGConstrainedDoFs::continuity_across_edges () const
+{
+  bool is_continuous = false;
+  for(unsigned int l=0; l<refinement_edge_indices.size(); ++l)
+    for(unsigned int i=0; i<refinement_edge_indices[l].size(); ++i)
+      if(refinement_edge_indices[l][i])
+      {
+        is_continuous = true;
+        return is_continuous;
+      }
+  return is_continuous;
+}
 
 DEAL_II_NAMESPACE_CLOSE
 
index 1c29f66813f7bfd887ffcf89ab8a3f8081d00f49..82dc1a1752fe0c37bc74ffcfad2d3009087c5ade 100644 (file)
@@ -644,15 +644,28 @@ namespace MeshWorker
                                          * Assemble a single matrix
                                          * into a global matrix.
                                          */
-    void assemble_transpose(MATRIX& G,
+
+    void assemble_up(MATRIX& G,
                            const FullMatrix<double>& M,
                            const std::vector<unsigned int>& i1,
-                           const std::vector<unsigned int>& i2);
+                           const std::vector<unsigned int>& i2,
+                            const unsigned int level = -1);
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+
+    void assemble_down(MATRIX& G,
+                           const FullMatrix<double>& M,
+                           const std::vector<unsigned int>& i1,
+                           const std::vector<unsigned int>& i2,
+                            const unsigned int level = -1);
 
                                         /**
                                          * Assemble a single matrix
                                          * into a global matrix.
                                          */
+
     void assemble_in(MATRIX& G,
                  const FullMatrix<double>& M,
                  const std::vector<unsigned int>& i1,
@@ -663,6 +676,7 @@ namespace MeshWorker
                                          * Assemble a single matrix
                                          * into a global matrix.
                                          */
+
     void assemble_out(MATRIX& G,
                  const FullMatrix<double>& M,
                  const std::vector<unsigned int>& i1,
@@ -1002,14 +1016,14 @@ namespace MeshWorker
                                          * the fine and the coarse
                                          * level at refinement edges.
                                          */
-       std::vector<MatrixPtr> interface_down;
+       std::vector<MatrixPtr> interface_out;
        
                                         /**
                                          * The interface matrix between
                                          * the coarse and the fine
                                          * level at refinement edges.
                                          */
-       std::vector<MatrixPtr> interface_up;    
+       std::vector<MatrixPtr> interface_in;    
       
       /**
        * A pointer to the object containing the block structure.
@@ -1562,19 +1576,72 @@ namespace MeshWorker
     
     template <class MATRIX>
     inline void
-    MGMatrixSimple<MATRIX>::assemble_transpose(
+    MGMatrixSimple<MATRIX>::assemble_up(
       MATRIX& G,
       const FullMatrix<double>& M,
       const std::vector<unsigned int>& i1,
-      const std::vector<unsigned int>& i2)
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
     {
       AssertDimension(M.n(), i1.size());
       AssertDimension(M.m(), i2.size());
       
-      for (unsigned int j=0; j<i1.size(); ++j)
-       for (unsigned int k=0; k<i2.size(); ++k)
-         if (std::fabs(M(k,j)) >= threshold)
-           G.add(i1[j], i2[k], M(k,j));
+      if(mg_constrained_dofs == 0)
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(k,j)) >= threshold)
+              G.add(i1[j], i2[k], M(k,j));
+      }
+      else
+      {
+        for (unsigned int j=0; j<i1.size(); ++j)
+          for (unsigned int k=0; k<i2.size(); ++k)
+            if (std::fabs(M(k,j)) >= threshold)
+            {
+              if(!mg_constrained_dofs->at_refinement_edge(level, i1[j]) ||
+                  mg_constrained_dofs->at_refinement_edge(level, i2[k]))
+              {
+                if(!mg_constrained_dofs->continuity_across_edges())
+                  G.add(i1[j], i2[k], M(k,j));
+              }
+            }
+      }
+    }
+
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_down(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2,
+      const unsigned int level)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+      if(mg_constrained_dofs == 0)
+      {
+        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)
+              G.add(i1[j], i2[k], M(j,k));
+      }
+      else
+      {
+        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->continuity_across_edges())
+                  G.add(i1[j], i2[k], M(j,k));
+              }
+            }
+      }
     }
 
     template <class MATRIX>
@@ -1615,7 +1682,7 @@ namespace MeshWorker
                     G.add(i1[j], i2[k], M(j,k));
                 }
                 else
-                    G.add(i1[j], i2[k], M(j,k));
+                  G.add(i1[j], i2[k], M(j,k));
               }
       }
     }
@@ -1671,10 +1738,12 @@ namespace MeshWorker
     {
       const unsigned int level = info.cell->level();
       assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
+
       if(mg_constrained_dofs != 0)
+      //if(interface_in != 0 && interface_out != 0 && level>0)
       {
-        assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
         assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
+        assemble_out((*interface_out)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
       }
     }
     
@@ -1712,8 +1781,13 @@ namespace MeshWorker
                                           // which is done by
                                           // the coarser cell
           assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices);
-          assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices);
-          assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
+          //assemble_transpose((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices);
+          //assemble((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
+          if(level1>0)
+          {
+          assemble_up((*flux_up)[level1],info1.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          assemble_down((*flux_down)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          }
        }
     }
     
index 5a7ea92d23743914ffc3612fe1d70cc1c7f42171..a0706b1a71caffd084d3793e672dc6ea4d69e652 100644 (file)
@@ -1495,6 +1495,9 @@ extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
       std::fill (cell_dofs.begin(), cell_dofs.end(), false);
       std::fill (boundary_cell_dofs.begin(), boundary_cell_dofs.end(), false);
 
+      if(fe.conforms(FiniteElementData<dim>::H1))
+      {
+        deallog << "conforms H1" << std::endl;
       for (unsigned int face_nr=0; face_nr<GeometryInfo<dim>::faces_per_cell; ++face_nr)
        {
          const typename DoFHandler<dim,spacedim>::face_iterator face = cell->face(face_nr);
@@ -1509,7 +1512,7 @@ extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
              if (neighbor->level() < cell->level())
                {
                  for (unsigned int j=0; j<dofs_per_face; ++j)
-                   cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
+                      cell_dofs[fe.face_to_cell_index(j,face_nr)] = true;
 
                  has_coarser_neighbor = true;
                }
@@ -1536,6 +1539,7 @@ extract_inner_interface_dofs (const MGDoFHandler<dim,spacedim> &mg_dof_handler,
            boundary_interface_dofs[level][local_dof_indices[i]] = true;
        }
     }
+    }
 }
 
 #endif

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.