]> https://gitweb.dealii.org/ - dealii.git/commitdiff
do some modifications to be able to use multigrid for Raviart Thomas
authorBaerbel Jannsen <baerbel.janssen@gmail.com>
Tue, 14 Sep 2010 14:53:00 +0000 (14:53 +0000)
committerBaerbel Jannsen <baerbel.janssen@gmail.com>
Tue, 14 Sep 2010 14:53:00 +0000 (14:53 +0000)
elements

git-svn-id: https://svn.dealii.org/trunk@21958 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_transfer_prebuilt.cc

index e5d0045dce570cafa827292e8fa027bc718956ca..91305357221b6313114067803bef7baae4505e48 100644 (file)
@@ -119,6 +119,12 @@ class MGConstrainedDoFs : public Subscriptor
     const std::vector<std::vector<bool> > &
       get_refinement_edge_boundary_indices () const;
 
+                                    /**
+                                     * Return if boundary_indices need to 
+                                      * be set or not.
+                                     */
+
+    const bool set_boundary_values () const;
   private:
 
                                     /**
@@ -263,6 +269,16 @@ MGConstrainedDoFs::get_refinement_edge_boundary_indices () const
   return refinement_edge_boundary_indices;
 }
 
+inline
+const bool 
+MGConstrainedDoFs::set_boundary_values () const
+{
+  const bool boundary_values_need_to_be_set 
+    = boundary_indices.size()!=0;
+  return boundary_values_need_to_be_set;
+}
+
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif
index dddff77ecb05b4e165921f2ee2e419c7364bbf84..1c29f66813f7bfd887ffcf89ab8a3f8081d00f49 100644 (file)
@@ -574,14 +574,14 @@ namespace MeshWorker
        
                                         /**
                                          * Initialize the matrices
-                                         * #interface_up and #interface_down
+                                         * #interface_in and #interface_out
                                          * used for local refinement
                                          * with continuous
                                          * Galerkin methods.
                                          */
 
-       void initialize_interfaces(MGLevelObject<MATRIX>& interface_up,
-                              MGLevelObject<MATRIX>& interface_down);
+       void initialize_interfaces(MGLevelObject<MATRIX>& interface_in,
+                              MGLevelObject<MATRIX>& interface_out);
                                         /**
                                          * Initialize the local data
                                          * in the
@@ -626,10 +626,19 @@ namespace MeshWorker
                                          * into a global matrix.
                                          */
     void assemble(MATRIX& G,
+                 const FullMatrix<double>& M,
+                 const std::vector<unsigned int>& i1,
+                 const std::vector<unsigned int>& i2);
+
+                                        /**
+                                         * Assemble a single matrix
+                                         * into a global matrix.
+                                         */
+    void assemble(MATRIX& G,
                  const FullMatrix<double>& M,
                  const std::vector<unsigned int>& i1,
                  const std::vector<unsigned int>& i2,
-                  const unsigned int level = -1);
+                  const unsigned int level);
        
                                         /**
                                          * Assemble a single matrix
@@ -638,14 +647,23 @@ namespace MeshWorker
     void assemble_transpose(MATRIX& G,
                            const FullMatrix<double>& M,
                            const std::vector<unsigned int>& i1,
-                           const std::vector<unsigned int>& i2,
-                            const unsigned int level = -1);
+                           const std::vector<unsigned int>& i2);
 
                                         /**
                                          * Assemble a single matrix
                                          * into a global matrix.
                                          */
-    void assemble_down(MATRIX& G,
+    void assemble_in(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_out(MATRIX& G,
                  const FullMatrix<double>& M,
                  const std::vector<unsigned int>& i1,
                  const std::vector<unsigned int>& i2,
@@ -680,7 +698,7 @@ namespace MeshWorker
                                          * refinement edge, coupling
                                          * coarse to fine.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_up;
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_in;
        
                                         /**
                                          * The matrix used for face
@@ -689,7 +707,7 @@ namespace MeshWorker
                                          * refinement edge, coupling
                                          * fine to coarse.
                                          */
-       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_down;
+       SmartPointer<MGLevelObject<MATRIX>,MGMatrixSimple<MATRIX> > interface_out;
       /**
        * A pointer to the object containing constraints.
        */
@@ -1463,10 +1481,10 @@ namespace MeshWorker
     template <class MATRIX>
     inline void
     MGMatrixSimple<MATRIX>::initialize_interfaces(
-      MGLevelObject<MATRIX>& up, MGLevelObject<MATRIX>& down)
+      MGLevelObject<MATRIX>& in, MGLevelObject<MATRIX>& out)
     {
-      interface_up = &up;
-      interface_down = &down;
+      interface_in = &in;
+      interface_out = &out;
     }
 
 
@@ -1480,6 +1498,23 @@ namespace MeshWorker
     }
 
 
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2)
+    {
+      AssertDimension(M.m(), i1.size());
+      AssertDimension(M.n(), i2.size());
+      
+        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));
+    }
+
 
     template <class MATRIX>
     inline void
@@ -1505,23 +1540,46 @@ 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])==false &&
-                  mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
-                if((mg_constrained_dofs->is_boundary_index(level, i1[j])==false &&
-                     mg_constrained_dofs->is_boundary_index(level, i2[k])==false)
-                    ||
-                    (mg_constrained_dofs->is_boundary_index(level, i1[j])==true &&
-                     mg_constrained_dofs->is_boundary_index(level, i2[k])==true
-                    &&
-                    i1[j] == i2[k])
-                    )
-                  G.add(i1[j], i2[k], M(j,k));
+              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())
+                {
+                  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));
+            }
       }
     }
     
+    
+    template <class MATRIX>
+    inline void
+    MGMatrixSimple<MATRIX>::assemble_transpose(
+      MATRIX& G,
+      const FullMatrix<double>& M,
+      const std::vector<unsigned int>& i1,
+      const std::vector<unsigned int>& i2)
+    {
+      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));
+    }
+
     template <class MATRIX>
     inline void
-    MGMatrixSimple<MATRIX>::assemble_down(
+    MGMatrixSimple<MATRIX>::assemble_in(
       MATRIX& G,
       const FullMatrix<double>& M,
       const std::vector<unsigned int>& i1,
@@ -1543,23 +1601,28 @@ 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])==true &&
-                  mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
-                if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false &&
-                     mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false)
-                    ||
-                    (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true &&
-                     mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true
-                     &&
-                    i1[j] == i2[k])
-                    )
-                  G.add(i1[j], i2[k], M(j,k));
+              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())
+                {
+                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]))
+                      ||
+                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                       mg_constrained_dofs->at_refinement_edge_boundary(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));
+              }
       }
     }
     
     template <class MATRIX>
     inline void
-    MGMatrixSimple<MATRIX>::assemble_transpose(
+    MGMatrixSimple<MATRIX>::assemble_out(
       MATRIX& G,
       const FullMatrix<double>& M,
       const std::vector<unsigned int>& i1,
@@ -1581,17 +1644,22 @@ 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])==true &&
-                  mg_constrained_dofs->at_refinement_edge(level, i2[k])==false)
-                if((mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==false &&
-                     mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==false) 
-                    ||
-                    (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j])==true &&
-                     mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])==true
-                    &&
-                    i1[j] == i2[k])
-                    )
-                  G.add(i1[j], i2[k], M(k,j));
+              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())
+                {
+                  if((!mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                        !mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k])) 
+                      ||
+                      (mg_constrained_dofs->at_refinement_edge_boundary(level, i1[j]) &&
+                       mg_constrained_dofs->at_refinement_edge_boundary(level, i2[k]) &&
+                       i1[j] == i2[k]))
+                    G.add(i1[j], i2[k], M(k,j));
+                }
+                else
+                    G.add(i1[j], i2[k], M(k,j));
+              }
       }
     }
     
@@ -1605,8 +1673,8 @@ namespace MeshWorker
       assemble((*matrix)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
       if(mg_constrained_dofs != 0)
       {
-        assemble_transpose((*interface_down)[level],info.matrix(0,false).matrix, info.indices, info.indices, level);
-        assemble_down((*interface_up)[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);
+        assemble_in((*interface_in)[level], info.matrix(0,false).matrix, info.indices, info.indices, level);
       }
     }
     
@@ -1629,6 +1697,13 @@ namespace MeshWorker
             assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices);
             assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices);
           }
+          else
+          {
+            assemble((*matrix)[level1], info1.matrix(0,false).matrix, info1.indices, info1.indices, level1);
+            assemble((*matrix)[level1], info1.matrix(0,true).matrix, info1.indices, info2.indices, level1);
+            assemble((*matrix)[level1], info2.matrix(0,false).matrix, info2.indices, info2.indices, level1);
+            assemble((*matrix)[level1], info2.matrix(0,true).matrix, info2.indices, info1.indices, level1);
+          }
        }
       else
        {
index 6a02a3d6db4b8cf9080901de8d1aa69c489962ca..d0f13cfae53b15ab8e3a6fa9ce52a9ee29e71163 100644 (file)
@@ -285,8 +285,8 @@ void MGTransferPrebuilt<VECTOR>::build_matrices (
                                // impose boundary conditions
                                // but only in the column of
                                // the prolongation matrix
-  if(mg_constrained_dofs != 0)
-  if (mg_constrained_dofs->get_boundary_indices().size() != 0)
+  if (mg_constrained_dofs != 0)
+  if (mg_constrained_dofs->set_boundary_values())
     {
       std::vector<unsigned int> constrain_indices;
       for (int level=n_levels-2; level>=0; --level)

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.