]> https://gitweb.dealii.org/ - dealii.git/commitdiff
debugging MG, mglocal.testcase produces errors, sorry
authorGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 30 Apr 2000 15:51:40 +0000 (15:51 +0000)
committerGuido Kanschat <dr.guido.kanschat@gmail.com>
Sun, 30 Apr 2000 15:51:40 +0000 (15:51 +0000)
git-svn-id: https://svn.dealii.org/trunk@2775 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/multigrid/mg_base.h
deal.II/deal.II/include/multigrid/multigrid.h
deal.II/deal.II/include/multigrid/multigrid.templates.h
deal.II/deal.II/source/multigrid/mg_base.cc
tests/deal.II/Makefile.in
tests/deal.II/mglocal.cc

index 4d6b2e755e88149e9e0c3824c2af520c1c5273dc..4e0396528e5d6ab1b915afe5e640f2a44d24ea7c 100644 (file)
@@ -391,7 +391,15 @@ class MGBase : public Subscriptor
                              Vector<double>       &dst,
                              const Vector<double> &src,
                              const Vector<double> &rhs) = 0;
-  
+
+                                    /**
+                                     * Print a level vector using
+                                     * @ref{DoFHandler}.
+                                     */
+    virtual void print_vector (const unsigned int level,
+                              const Vector<double>& v,
+                              const char* name) const = 0;
+    
   private:
                                     /**
                                      * Auxiliary vector.
index 7cfa21e1aee799d10b4b2ca6886f9c2d116ad894..1098c4fb967cbcad7bd5fa98c04d24bfb2d6d164 100644 (file)
@@ -105,6 +105,14 @@ class Multigrid : public MGBase
                              const Vector<double> &src,
                              const Vector<double> &rhs);
 
+                                    /**
+                                     * Print a level vector using
+                                     * @ref{DoFHandler}.
+                                     */
+    virtual void print_vector (const unsigned int level,
+                              const Vector<double>& v,
+                              const char* name) const;
+    
                                     /**
                                      * Read-only access to level matrices.
                                      */
index 77d7be4f8b49ed8f4a35b409bae39798181681a2..c11731dc35be1023538dd341b7251af75b465bad 100644 (file)
@@ -15,6 +15,7 @@
 
 
 #include <dofs/dof_constraints.h>
+#include <numerics/data_out.h>
 #include <multigrid/multigrid.h>
 #include <algorithm>
 #include <fstream>
@@ -154,6 +155,9 @@ Multigrid<dim>::copy_from_mg(Vector<number> &dst) const
                                   // traverse all cells and copy the
                                   // data appropriately to the output
                                   // vector
+
+                                  // Is the level monotonuosly increasing?
+
   for (; level_cell != endc; ++level_cell, ++global_cell)
     {
       const unsigned int level = level_cell->level();
@@ -176,4 +180,51 @@ Multigrid<dim>::copy_from_mg(Vector<number> &dst) const
 }
 
 
+
+template <int dim>
+void
+Multigrid<dim>::print_vector (const unsigned int level,
+                             const Vector<double>& v,
+                             const char* name) const
+{
+  Vector<double> out_vector;
+  
+  const DoFHandler<dim>* dof = mg_dof_handler;
+  
+  const unsigned int dofs_per_cell = mg_dof_handler->get_fe().dofs_per_cell;
+
+  vector<unsigned int> global_dof_indices (dofs_per_cell);
+  vector<unsigned int> level_dof_indices (dofs_per_cell);
+
+  DoFHandler<dim>::cell_iterator
+    global_cell = dof->begin(level);
+  MGDoFHandler<dim>::cell_iterator
+    level_cell = mg_dof_handler->begin(level);
+  const MGDoFHandler<dim>::cell_iterator
+    endc = mg_dof_handler->end(level);
+
+                                  // traverse all cells and copy the
+                                  // data appropriately to the output
+                                  // vector
+  for (; level_cell != endc; ++level_cell, ++global_cell)
+    {
+      global_cell->get_dof_indices (global_dof_indices);
+      level_cell->get_mg_dof_indices(level_dof_indices);
+
+                                      // copy level-wise data to
+                                      // global vector
+      for (unsigned int i=0; i<dofs_per_cell; ++i)
+       out_vector(global_dof_indices[i])
+         = v(level_dof_indices[i]);
+    }
+
+  ofstream out_file(name);
+  DataOut<dim> out;
+  out.attach_dof_handler(*dof);
+  out.add_data_vector(out_vector, "v");
+  out.build_patches(5);
+  out.write_gnuplot(out_file);
+}
+
 #endif
+
index 9b9e927d3b0c33b5f182b8c01209301e7d47a22b..630e8c3ef1bd1b3c3d0a720f14a514497ca8ce4e 100644 (file)
 #include <cmath>
 
 
-//TODO: this function is only for debugging purposes and should be removed sometimes
-template<class VECTOR>
-static
-void print_vector(ostream& s, const VECTOR& v, const char* text)
-{
-  const unsigned int n = (unsigned int)(sqrt(v.size())+.3);
-  unsigned int k=0;
-
-  s << endl << "splot '-' title '" << text << "'" << endl;
-  
-                                  // write the vector entries in a
-                                  // kind of square
-  for (unsigned int i=0;i<n;++i)
-    {
-      for (unsigned int j=0;j<n;++j)
-       s << '\n' << v(k++);
-      s << endl;
-    }
-  s << "e\npause -1" << endl;
-}
-
-
-//////////////////////////////////////////////////////////////////////
-
-
 MGBase::~MGBase () 
-{};
+{}
 
 
 MGBase::MGBase(const MGTransferBase &transfer,
@@ -59,7 +34,7 @@ MGBase::MGBase(const MGTransferBase &transfer,
 {
   Assert(minlevel <= maxlevel,
         ExcSwitchedLevels(minlevel, maxlevel));
-};
+}
 
 
 void
@@ -67,11 +42,9 @@ MGBase::vcycle(const MGSmootherBase     &pre_smooth,
               const MGSmootherBase     &post_smooth,
               const MGCoarseGridSolver &coarse_grid_solver)
 {
-//  static int k=0;
-//  cout << "set title 'cycle " << ++k << "'\n";
-  
-  level_mgstep(maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
-};
+  level_mgstep (maxlevel, pre_smooth, post_smooth, coarse_grid_solver);
+  abort ();
+}
 
 
 void
@@ -80,6 +53,11 @@ MGBase::level_mgstep(const unsigned int        level,
                     const MGSmootherBase     &post_smooth,
                     const MGCoarseGridSolver &coarse_grid_solver)
 {
+  char *name = new char[100];
+
+  sprintf(name, "MG%d-defect",level);
+  print_vector(level, defect[level], name);
+  
   solution[level] = 0.;
   
   if (level == minlevel)
@@ -91,15 +69,17 @@ MGBase::level_mgstep(const unsigned int        level,
                           // smoothing of the residual by modifying s
   pre_smooth.smooth(level, solution[level], defect[level]);
                                   // t = d-As
+
+  sprintf(name, "MG%d-pre",level);
+  print_vector(level, solution[level], name);
+  
   t.reinit(solution[level].size());
   level_vmult(level, t, solution[level], defect[level]);
-//  print_vector(cout,t,"T");
   
                                   // make t rhs of lower level
 //TODO: this function adds the restricted t to defect[level-1].
 //TODO: why don't we have to clear it before?  
   transfer->restrict_and_add (level, defect[level-1], t);
-//  print_vector(cout,defect[level-1],"Dl-1");
   
                                   // do recursion
   level_mgstep(level-1, pre_smooth, post_smooth, coarse_grid_solver);
@@ -111,16 +91,22 @@ MGBase::level_mgstep(const unsigned int        level,
   t.reinit(solution[level].size());
 
                                   // do coarse grid correction
-//  print_vector(cout,solution[level-1],"Sl-1");
+
   transfer->prolongate(level, t, solution[level-1]);
-//  print_vector(cout,t,"T");
+
+  sprintf(name, "MG%d-cgc",level);
+  print_vector(level, t, name);
+
   solution[level] += t;
   
                                   // smoothing (modify solution again)
 //TODO: what happens here? smooth overwrites the solution[level],
 //TODO: so the previous two statements should have no effect. No?  
   post_smooth.smooth(level, solution[level], defect[level]);
-};
+
+  sprintf(name, "MG%d-post",level);
+  print_vector(level, solution[level], name);
+}
 
 
 //////////////////////////////////////////////////////////////////////
index d9a50cd784892c6a419dba6246e3a9015bc0452b..167c34ce069c2aa1e40c933b54e55d46d1a4c019 100644 (file)
@@ -65,7 +65,8 @@ endif
        @$(CXX) $(flags) -c $< -o $@
 
 
-all: grid_test.check dof_test.check fe_tables.check gradients.check \
+all: grid_generator.check grid_test.check dof_test.check \
+     fe_tables.check gradients.check \
      constraints.check mg.check mglocal.check wave-test-3.check second_derivatives.check
 
 exe: $(all:.check=.testcase) benchmark
@@ -171,6 +172,21 @@ grid_test.testcase: $(grid_test-o-files) $(libdeal3d) $(libraries)
        @$(CXX) $(flags) -o $@ $^
 
 
+############################################################
+
+grid_generator-cc-files = grid_generator.cc 
+
+ifeq ($(debug-mode),on)
+grid_generator-o-files = $(grid_generator-cc-files:.cc=.go)
+else
+grid_generator-o-files = $(grid_generator-cc-files:.cc=.o)
+endif
+
+grid_generator.testcase: $(grid_generator-o-files) $(libdeal3d) $(libraries)
+       @echo =====linking======= $<
+       @$(CXX) $(flags) -o $@ $^
+
+
 ############################################################
 
 second_derivatives-cc-files = second_derivatives.cc 
index 5c7ca307397574c9d223abfae4837362b99a140e..35312ffb1d068ba270e14291c2af2f899b095780 100644 (file)
@@ -28,7 +28,7 @@
 #include <dofs/dof_accessor.h>
 #include <multigrid/mg_dof_accessor.h>
 #include <grid/grid_generator.h>
-#include <numerics/data_io.h>
+#include <numerics/data_out.h>
 #include <fe/fe_lib.lagrange.h>
 #include <fe/fe_values.h>
 #include <multigrid/multigrid.h>
@@ -69,10 +69,10 @@ extern void write_gnuplot (const MGDoFHandler<2>& dofs,
 int main()
 {
   ofstream logfile("mglocal.output");
-  logfile.setf(ios::fixed);
-  logfile.precision (3);
+//  logfile.setf(ios::fixed);
+//  logfile.precision (3);
   deallog.attach(logfile);
-  deallog.depth_console(0);
+//  deallog.depth_console(0);
 
   Helmholtz equation;
   RHSFunction<2> rhs;
@@ -82,7 +82,7 @@ int main()
   FEQ2<2> fe2;
   FEQ3<2> fe3;
   FEQ4<2> fe4;
-  for (unsigned int degree=1;degree<=1;degree++)
+  for (unsigned int degree=1;degree<=4;degree++)
     {
       Triangulation<2> tr;
       MGDoFHandler<2> mgdof(tr);
@@ -104,7 +104,7 @@ int main()
       cell->set_refine_flag();
       tr.execute_coarsening_and_refinement();
 
-      tr.refine_global(2);
+      tr.refine_global(1);
       dof.distribute_dofs(*fe);
       const unsigned int size = dof.n_dofs();
       deallog << "DoFs " << size << endl;
@@ -133,8 +133,8 @@ int main()
          Vector<double> u;
          u.reinit(f);
          PrimitiveVectorMemory<> mem;
-         SolverControl control(20, 1.e-12, true);
-         SolverRichardson<> solver(control, mem);
+         SolverControl control(20, 1.e-8, true);
+         SolverCG<> solver(control, mem);
          
          MGLevelObject<SparsityPattern> mgstruct(0, tr.n_levels()-1);
          MGLevelObject<SparseMatrix<double> > mgA(0,tr.n_levels()-1);
@@ -173,6 +173,13 @@ Multigrid<2> multigrid(mgdof, hanging_nodes, mgstruct, mgA, transfer, tr.n_level
           
          solver.solve(A, u, f, mgprecondition);
          hanging_nodes.distribute(u);
+
+         DataOut<2> out;
+         ofstream ofile("out.gnuplot");
+         out.attach_dof_handler(dof);
+         out.add_data_vector(u,"u");
+         out.build_patches(4);
+         out.write_gnuplot(ofile);
        }
       deallog.pop();
     }

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.