geometric_multigrid.h
Go to the documentation of this file.
1 // Include guards
2 #ifndef OOMPH_GEOMETRIC_MULTIGRID_HEADER
3 #define OOMPH_GEOMETRIC_MULTIGRID_HEADER
4 
5 // Config header generated by autoconfig
6 #ifdef HAVE_CONFIG_H
7 #include <config.h>
8 #endif
9 
10 // Oomph-lib headers
11 #include "problem.h"
12 #include "matrices.h"
14 #include "preconditioner.h"
15 
16 // Namespace extension
17 namespace oomph
18 {
19 
20 //======================================================================
21 /// MGProblem class; subclass of Problem
22 //======================================================================
23  class MGProblem : public virtual Problem
24  {
25 
26  public:
27 
28  /// Constructor. Initialise pointers to coarser and finer levels
30  {}
31 
32  /// Destructor (empty)
33  virtual ~MGProblem()
34  {}
35 
36  /// \short This function needs to be implemented in the derived problem:
37  /// Returns a pointer to a new object of the same type as the derived
38  /// problem
39  virtual MGProblem* make_new_problem()=0;
40 
41  /// \short Function to get a pointer to the mesh we will be working
42  /// with. If there are flux elements present in the mesh this will
43  /// be overloaded to return a pointer to the bulk mesh otherwise
44  /// it can be overloaded to point to the global mesh but it must
45  /// be of type RefineableMeshBase
47 
48  }; // End of MGProblem class
49 
50 
51 //////////////////////////////////////////////////////////
52 //////////////////////////////////////////////////////////
53 //////////////////////////////////////////////////////////
54 
55 
56 //======================================================================
57 // MG solver class
58 //======================================================================
59  template<unsigned DIM>
61  {
62  public:
63 
64  /// \short typedef for a function that returns a pointer to an object
65  /// of the class Smoother to be used as the pre-smoother
66  typedef Smoother* (*PreSmootherFactoryFctPt)();
67 
68  /// \short typedef for a function that returns a pointer to an object
69  /// of the class Smoother to be used as the post-smoother
70  typedef Smoother* (*PostSmootherFactoryFctPt)();
71 
72  /// Access function to set the pre-smoother creation function.
74  {
75  // Assign the function pointer
76  Pre_smoother_factory_function_pt=pre_smoother_fn;
77  }
78 
79  /// Access function to set the post-smoother creation function.
81  PostSmootherFactoryFctPt post_smoother_fn)
82  {
83  // Assign the function pointer
84  Post_smoother_factory_function_pt=post_smoother_fn;
85  }
86 
87  /// \short Constructor: Set up default values for number of V-cycles
88  /// and pre- and post-smoothing steps.
89  MGSolver(MGProblem* mg_problem_pt) :
90  Nvcycle(200),
91  Mg_problem_pt(mg_problem_pt),
93  Suppress_all_output(false),
94  Stream_pt(0),
97  Npre_smooth(2),
98  Npost_smooth(2),
99  Doc_everything(false),
100  Has_been_setup(false),
101  Has_been_solved(false)
102  {
103  // Set the tolerance in the base class
104  this->Tolerance=1.0e-09;
105 
106  // Set the pointer to the finest level as the first
107  // entry in Mg_hierarchy
108  Mg_hierarchy.push_back(Mg_problem_pt);
109  } // End of MGSolver
110 
111  /// Delete any dynamically allocated data
113  {
114  // Run the function written to clean up the memory
115  clean_up_memory();
116  } // End of ~MGSolver
117 
118  /// Clean up anything that needs to be cleaned up
120  {
121  // We only need to destroy data if the solver has been set up and
122  // the data hasn't already been cleared
123  if (Has_been_setup)
124  {
125  // Loop over all of the levels in the hierarchy
126  for (unsigned i=0;i<Nlevel-1;i++)
127  {
128  // Delete the pre-smoother associated with this level
129  delete Pre_smoothers_storage_pt[i];
130 
131  // Make it a null pointer
133 
134  // Delete the post-smoother associated with this level
136 
137  // Make it a null pointer
139 
140  // Delete the system matrix associated with the i-th level
141  delete Mg_matrices_storage_pt[i];
142 
143  // Make it a null pointer
145  }
146 
147  // Loop over all but the coarsest of the levels in the hierarchy
148  for (unsigned i=0;i<Nlevel-1;i++)
149  {
150  // Delete the interpolation matrix associated with the i-th level
152 
153  // Make it a null pointer
155 
156  // Delete the restriction matrix associated with the i-th level
158 
159  // Make it a null pointer
161  }
162 
163  // If this solver has been set up then a hierarchy of problems
164  // will have been set up. If the user chose to document everything
165  // then the coarse-grid multigrid problems will have been kept alive
166  // which means we now have to loop over the coarse-grid levels and
167  // destroy them
168  if (Doc_everything)
169  {
170  // Loop over the levels
171  for (unsigned i=1;i<Nlevel;i++)
172  {
173  // Delete the i-th level problem
174  delete Mg_hierarchy[i];
175 
176  // Make the associated pointer a null pointer
177  Mg_hierarchy[i]=0;
178  }
179  } // if (Doc_everything)
180 
181  // Everything has been deleted now so we need to indicate that the
182  // solver is not set up
183  Has_been_setup=false;
184  } // if (Has_been_setup)
185  } // End of clean_up_memory
186 
187  /// \short Makes a vector which will be used in the self-test. Is currently
188  /// set to make the entries of the vector represent a plane wave propagating
189  /// at an angle of 45 degrees
190  void set_self_test_vector();
191 
192  /// \short Makes a vector, restricts it down the levels of the hierarchy
193  /// and documents it at each level. After this is done the vector is
194  /// interpolated up the levels of the hierarchy with the output
195  /// being documented at each level
196  void self_test();
197 
198  /// \short Make a self-test to make sure that the interpolation matrices
199  /// are doing the same thing to restrict the vectors down through the
200  /// heirachy.
201  void restriction_self_test();
202 
203  /// \short Make a self-test to make sure that the interpolation matrices
204  /// are doing the same thing to interpolate the vectors up.
206 
207  /// \short Given a level in the hierarchy, an input vector and a filename
208  /// this function will document the given vector according to the structure
209  /// of the mesh on the given level
210  void plot(const unsigned& hierarchy_level,
211  const DoubleVector& input_vector,
212  const std::string& filename);
213 
214  /// \short Disable all output from mg_solve apart from the number of
215  /// V-cycles used to solve the problem
217  {
218  // Set the value of Doc_time (inherited from LinearSolver) to false
219  Doc_time=false;
220 
221  // Enable the suppression of output from the V-cycle
223  } // End of disable_v_cycle_output
224 
225  /// \short Suppress anything that can be suppressed, i.e. any timings.
226  /// Things like mesh adaptation can not however be silenced using this
228  {
229  // Set the value of Doc_time (inherited from LinearSolver) to false
230  Doc_time=false;
231 
232  // Enable the suppression of output from the V-cycle
234 
235  // Enable the suppression of everything
236  Suppress_all_output=true;
237 
238  // Store the output stream pointer
240 
241  // Now set the oomph_info stream pointer to the null stream to
242  // disable all possible output
244  } // End of disable_output
245 
246  /// \short Enable the output of the V-cycle timings and other output
248  {
249  // Enable time documentation
250  Doc_time=true;
251 
252  // Enable output from the MG solver
254  } // End of enable_v_cycle_output
255 
256  /// \short Enable the output from anything that could have been suppressed
258  {
259  // Enable the documentation of everything (if this is set to TRUE then
260  // the function self_test() will be run which outputs a solution
261  // represented on each level of the hierarchy
262  Doc_everything=true;
263  } // End of enable_doc_everything
264 
265  /// \short Enable the output from anything that could have been suppressed
267  {
268  // Enable time documentation
269  Doc_time=true;
270 
271  // Enable output from everything during the full setup of the solver
272  Suppress_all_output=false;
273 
274  // Enable output from the MG solver
276  } // End of enable_output
277 
278  /// \short Suppress the output of both smoothers and SuperLU
280  {
281  // Loop over all levels of the hierarchy
282  for (unsigned i=0;i<Nlevel-1;i++)
283  {
284  // Disable time documentation on each level (for each pre-smoother)
285  Pre_smoothers_storage_pt[i]->disable_doc_time();
286 
287  // Disable time documentation on each level (for each post-smoother)
288  Post_smoothers_storage_pt[i]->disable_doc_time();
289  }
290 
291  // We only do a direct solve on the coarsest level so this is the only
292  // place we need to silence SuperLU
293  Mg_matrices_storage_pt[Nlevel-1]->linear_solver_pt()->disable_doc_time();
294  } // End of disable_smoother_and_superlu_doc_time
295 
296  /// Return the number of post-smoothing iterations (lvalue)
297  unsigned& npost_smooth()
298  {
299  // Return the number of post-smoothing iterations to be done on each
300  // level of the hierarchy
301  return Npost_smooth;
302  } // End of npost_smooth
303 
304  /// Return the number of pre-smoothing iterations (lvalue)
305  unsigned& npre_smooth()
306  {
307  // Return the number of pre-smoothing iterations to be done on each
308  // level of the hierarchy
309  return Npre_smooth;
310  } // End of npre_smooth
311 
312  /// \short Pre-smoother: Perform 'max_iter' smoothing steps on the
313  /// linear system Ax=b with current RHS vector, b, starting with
314  /// current solution vector, x. Return the residual vector r=b-Ax.
315  /// Uses the default smoother (set in the MGProblem constructor)
316  /// which can be overloaded for a specific problem.
317  void pre_smooth(const unsigned& level)
318  {
319  // Run pre-smoother 'max_iter' times
320  Pre_smoothers_storage_pt[level]->
321  smoother_solve(Rhs_mg_vectors_storage[level],
322  X_mg_vectors_storage[level]);
323 
324  // Calculate the residual r=b-Ax and assign it
325  Mg_matrices_storage_pt[level]->
326  residual(X_mg_vectors_storage[level],
327  Rhs_mg_vectors_storage[level],
329  } // End of pre_smooth
330 
331  /// \short Post-smoother: Perform max_iter smoothing steps on the
332  /// linear system Ax=b with current RHS vector, b, starting with
333  /// current solution vector, x. Uses the default smoother (set in
334  /// the MGProblem constructor) which can be overloaded for specific
335  /// problem.
336  void post_smooth(const unsigned& level)
337  {
338  // Run post-smoother 'max_iter' times
340  smoother_solve(Rhs_mg_vectors_storage[level],
341  X_mg_vectors_storage[level]);
342  } // End of post_smooth
343 
344  /// \short Return norm of residual r=b-Ax and the residual vector itself
345  /// on the level-th level
346  double residual_norm(const unsigned& level)
347  {
348  // And zero the entries of residual
349  Residual_mg_vectors_storage[level].initialise(0.0);
350 
351  // Get the residual
352  Mg_matrices_storage_pt[level]->residual(X_mg_vectors_storage[level],
353  Rhs_mg_vectors_storage[level],
355 
356  // Return the norm of the residual
357  return Residual_mg_vectors_storage[level].norm();
358  } // End of residual_norm
359 
360  /// \short Call the direct solver (SuperLU) to solve the problem exactly.
361  // The result is placed in X_mg
363  {
364  // Get solution by direct solve:
368  } // End of direct_solve
369 
370  /// \short Builds a CRDoubleMatrix that is used to interpolate the
371  /// residual between levels. The transpose can be used as the full
372  /// weighting restriction.
373  void interpolation_matrix_set(const unsigned& level,
374  double* value,
375  int* col_index,
376  int* row_st,
377  unsigned& ncol,
378  unsigned& nnz)
379  {
380  // Dynamically allocate the interpolation matrix
382 
383  // Build the matrix
385  build_without_copy(ncol,nnz,value,col_index,row_st);
386  } // End of interpolation_matrix_set
387 
388  /// \short Builds a CRDoubleMatrix that is used to interpolate the
389  /// residual between levels. The transpose can be used as the full
390  /// weighting restriction.
391  void interpolation_matrix_set(const unsigned& level,
392  Vector<double>& value,
393  Vector<int>& col_index,
394  Vector<int>& row_st,
395  unsigned& ncol,
396  unsigned& nrow)
397  {
398  // Dynamically allocate the interpolation matrix
400 
401  // Make the distribution pointer
402  LinearAlgebraDistribution* dist_pt=
404  communicator_pt(),nrow,false);
405 
406 #ifdef PARANOID
407 #ifdef OOMPH_HAS_MPI
408  // Set up the warning messages
409  std::string warning_message="Setup of interpolation matrix distribution ";
410  warning_message+="has not been tested with MPI.";
411 
412  // If we're not running the code in serial
413  if (dist_pt->communicator_pt()->nproc()>1)
414  {
415  // Throw a warning
416  OomphLibWarning(warning_message,
417  OOMPH_CURRENT_FUNCTION,
418  OOMPH_EXCEPTION_LOCATION);
419  }
420 #endif
421 #endif
422 
423  // Build the matrix itself
425  build(dist_pt,ncol,value,col_index,row_st);
426 
427  // Delete the newly created distribution pointer
428  delete dist_pt;
429 
430  // Make it a null pointer
431  dist_pt=0;
432  } // End of interpolation_matrix_set
433 
434  /// \short Builds a CRDoubleMatrix on each level that is used to
435  /// restrict the residual between levels. The transpose can be used
436  /// as the interpolation matrix
438  {
439  for (unsigned i=0;i<Nlevel-1;i++)
440  {
441  // Dynamically allocate the restriction matrix
443 
444  // Set the restriction matrix to be the transpose of the
445  // interpolation matrix
447  get_matrix_transpose(Restriction_matrices_storage_pt[i]);
448  }
449  } // End of set_restriction_matrices_as_interpolation_transposes
450 
451  /// \short Restrict residual (computed on level-th MG level) to the next
452  /// coarser mesh and stick it into the coarse mesh RHS vector.
453  void restrict_residual(const unsigned& level);
454 
455  /// \short Interpolate solution at current level onto next finer mesh
456  /// and correct the solution x at that level
457  void interpolate_and_correct(const unsigned& level);
458 
459  /// \short Given the son_type of an element and a local node number
460  /// j in that element with nnode_1d nodes per coordinate direction,
461  /// return the local coordinate s in its father element. Needed in
462  /// the setup of the interpolation matrices
463  void level_up_local_coord_of_node(const int& son_type,
464  Vector<double>& s);
465 
466  /// \short Setup the interpolation matrix on each level
468 
469  /// \short Setup the interpolation matrix on each level (used for
470  /// unstructured meshes)
472 
473  /// \short Setup the transfer matrices on each level
475 
476  /// \short Do a full setup (assumes everything will be setup around the
477  /// MGProblem pointer given in the constructor)
478  void full_setup();
479 
480  /// \short Virtual function in the base class that needs to be implemented
481  /// later but for now just leave it empty
482  void solve(Problem* const& problem_pt, DoubleVector& result)
483  {
484  // Dynamically cast problem_pt of type Problem to a MGProblem pointer
485  MGProblem* mg_problem_pt=dynamic_cast<MGProblem*>(problem_pt);
486 
487 #ifdef PARANOID
488  // PARANOID check - If the dynamic_cast produces a null pointer the
489  // input was not a MGProblem
490  if (0==mg_problem_pt)
491  {
492  throw OomphLibError("Input problem must be of type MGProblem.",
493  OOMPH_CURRENT_FUNCTION,
494  OOMPH_EXCEPTION_LOCATION);
495  }
496  // PARANOID check - If a node in the input problem has more than
497  // one value we cannot deal with it so arbitarily check the first one
498  if (problem_pt->mesh_pt()->node_pt(0)->nvalue()!=1)
499  {
500  // How many dofs are there in the first node
501  unsigned n_value=problem_pt->mesh_pt()->node_pt(0)->nvalue();
502 
503  // Make the error message
504  std::ostringstream error_message_stream;
505  error_message_stream << "Cannot currently deal with more than 1 dof"
506  << " per node. This problem has " << n_value
507  << " dofs per node."
508  << std::endl;
509 
510  // Throw the error message
511  throw OomphLibError(error_message_stream.str(),
512  OOMPH_CURRENT_FUNCTION,
513  OOMPH_EXCEPTION_LOCATION);
514  }
515 #endif
516 
517  // Assign the new MGProblem pointer to Mg_problem_pt
518  Mg_problem_pt=mg_problem_pt;
519 
520  // Set up all of the required MG structures
521  full_setup();
522 
523  // Run the MG method and assign the solution to result
524  mg_solve(result);
525 
526  // Only output if the V-cycle output isn't suppressed
528  {
529  // Notify user that the hierarchy of levels is complete
530  oomph_info << "\n================="
531  << "Multigrid Solve Complete"
532  << "=================\n"
533  << std::endl;
534  }
535 
536  // If the user did not request all output be suppressed
537  if (!Suppress_all_output)
538  {
539  // If the user requested all V-cycle output be suppressed
541  {
542  // Create an extra line spacing
543  oomph_info << std::endl;
544  }
545 
546  // Output number of V-cycles taken to solve
547  if (Has_been_solved)
548  {
549  oomph_info << "Total number of V-cycles required for solve: "
550  << V_cycle_counter << std::endl;
551  }
552  else
553  {
554  oomph_info << "Total number of V-cycles used: "
555  << V_cycle_counter << std::endl;
556  }
557  } // if (!Suppress_all_output)
558 
559  // Only enable and assign the stream pointer again if we originally
560  // suppressed everything otherwise it won't be set yet
562  {
563  // Now enable the stream pointer again
565  }
566  } // End of solve
567 
568  /// Number of iterations
569  unsigned iterations() const
570  {
571  // Return the number of V-cycles which have been done
572  return V_cycle_counter;
573  } // End of iterations
574 
575  protected:
576 
577  /// \short Do the actual solve -- this is called through the pure virtual
578  /// solve function in the LinearSolver base class. The function is stored
579  /// as protected to allow the MGPreconditioner derived class to use the
580  /// solver
581  void mg_solve(DoubleVector& result);
582 
583  /// \short Normalise the rows of the restriction matrices to avoid
584  /// amplifications when projecting to the coarser level
586 
587  /// \short Maximum number of V-cycles (this is set as a protected variable so
588  // that it can be changed in the MGPreconditioner class)
589  unsigned Nvcycle;
590 
591  /// \short Pointer to the MG problem (deep copy). This is protected to provide
592  /// access to the MG preconditioner
594 
595  /// \short Vector to store the RHS vectors (Rhs_mg). This is protected to
596  /// allow the multigrid preconditioner to assign the RHS vector during
597  /// preconditioner_solve()
599 
600  /// \short Indicates whether or not the V-cycle output should be
601  /// suppressed. Needs to be protected member data for the multigrid
602  /// preconditioner to know whether or not to output information
603  /// with each preconditioning step
605 
606  /// \short If this is set to true then all output from the solver is
607  /// suppressed. This is protected member data so that the multigrid
608  /// preconditioner knows whether or not to restore the stream pointer
610 
611  /// \short Pointer to the output stream -- defaults to std::cout.
612  /// This is protected member data to allow the preconditioner to
613  /// restore normal output if everything was chosen to be
614  /// suppressed by the user
615  std::ostream* Stream_pt;
616 
617  private:
618 
619  /// \short Function to create pre-smoothers
621 
622  /// \short Function to create post-smoothers
624 
625  /// \short Function to set up the hierachy of levels. Creates a vector
626  /// of pointers to each MG level
627  void setup_mg_hierarchy();
628 
629  /// \short Function to set up the hierachy of levels. Creates a vector
630  /// of pointers to each MG level
631  void setup_mg_structures();
632 
633  /// \short Function to set up all of the smoothers once the system matrices
634  /// have been set up
635  void setup_smoothers();
636 
637  /// The number of levels in the multigrid heirachy
638  unsigned Nlevel;
639 
640  /// Vector containing pointers to problems in hierarchy
642 
643  /// Vector to store the system matrices
645 
646  /// Vector to store the interpolation matrices
648 
649  /// Vector to store the restriction matrices
651 
652  /// Vector to store the solution vectors (X_mg)
654 
655  /// Vector to store the residual vectors
657 
658  /// \short Vector to store the result of interpolation on each level (only
659  /// required if the user wishes to document the output of interpolation
660  /// and restriction on each level)
662 
663  /// \short Vector to store the result of restriction on each level (only
664  /// required if the user wishes to document the output of interpolation
665  /// and restriction on each level)
667 
668  /// Vector to store the pre-smoothers
670 
671  /// Vector to store the post-smoothers
673 
674  /// Number of pre-smoothing steps
675  unsigned Npre_smooth;
676 
677  /// Number of post-smoothing steps
678  unsigned Npost_smooth;
679 
680  /// \short If this is set to true we document everything. In addition
681  /// to outputting the information of the setup timings and V-cycle
682  /// data we document the refinement and unrefinement patterns given
683  /// by the transfer operators which is done by keeping the coarser
684  /// MG problem pointers alive
686 
687  /// Boolean variable to indicate whether or not the solver has been setup
689 
690  /// Boolean variable to indicate whether or not the problem was
691  /// successfully solved
693 
694  /// Pointer to counter for V-cycles
695  unsigned V_cycle_counter;
696  };
697 
698 //====================================================================
699 /// An interface to allow scalar MG to be used as a Preconditioner
700 //====================================================================
701  template<unsigned DIM>
702  class MGPreconditioner : public MGSolver<DIM>, public Preconditioner
703  {
704  public:
705 
706  /// Constructor.
707  MGPreconditioner(MGProblem* mg_problem_pt) : MGSolver<DIM>(mg_problem_pt)
708  {
709  // Set the number of V-cycles to be 1 (as expected as a preconditioner)
710  this->Nvcycle=2;
711  } // End of MGPreconditioner (constructor)
712 
713  /// Destructor (empty)
715 
716  /// Broken copy constructor.
718  {
719  BrokenCopy::broken_copy("MGPreconditioner");
720  }
721 
722  /// Broken assignment operator.
724  {
725  BrokenCopy::broken_assign("MGPreconditioner");
726  }
727 
728  /// \short Function to set up a preconditioner for the linear system
729  void setup()
730  {
731  this->full_setup();
732 
733  // Only enable and assign the stream pointer again if we originally
734  // suppressed everything otherwise it won't be set yet
735  if (this->Suppress_all_output)
736  {
737  // Now enable the stream pointer again
739  }
740  } // End of setup
741 
742  /// \short Function applies MG to the vector r for a full solve
743  virtual void preconditioner_solve(const DoubleVector& rhs,DoubleVector &z)
744  {
745 #ifdef PARANOID
746  if (this->Mg_problem_pt->ndof()!=rhs.nrow())
747  {
748  throw OomphLibError("Matrix and RHS vector sizes incompatible.",
749  OOMPH_CURRENT_FUNCTION,
750  OOMPH_EXCEPTION_LOCATION);
751  }
752 #endif
753 
754  // Set the right-hand side vector on the finest level to r
755  this->Rhs_mg_vectors_storage[0]=rhs;
756 
757  // Run the MG method and assign the solution to z
758  this->mg_solve(z);
759 
760  // Only output if the V-cycle output isn't suppressed
761  if (!(this->Suppress_v_cycle_output))
762  {
763  // Notify user that the hierarchy of levels is complete
764  oomph_info << "\n==========Multigrid Preconditioner Solve Complete========="
765  << "\n" << std::endl;
766  }
767 
768  // Only enable and assign the stream pointer again if we originally
769  // suppressed everything otherwise it won't be set yet
770  if (this->Suppress_all_output)
771  {
772  // Now enable the stream pointer again
774  }
775  } // End of preconditioner_solve
776 
777  /// Clean up memory
779  };
780 
781 
782 //===================================================================
783 /// Runs a full setup of the MG solver
784 //===================================================================
785  template<unsigned DIM>
787  {
788  // Initialise the timer start variable
789  double t_fs_start=0.0;
790 
791  // If we're allowed to output
792  if (!Suppress_all_output)
793  {
794  // Start the timer
795  t_fs_start=TimingHelpers::timer();
796 
797  // Notify user that the hierarchy of levels is complete
798  oomph_info << "\n===============Starting Multigrid Full Setup=============="
799  << std::endl;
800 
801  // Notify user that the hierarchy of levels is complete
802  oomph_info << "\nStarting the full setup of the multigrid solver."
803  << std::endl;
804  }
805 
806 #ifdef PARANOID
807  // PARANOID check - Make sure the dimension of the solver matches the
808  // dimension of the elements used in the problems mesh
809  if (dynamic_cast<FiniteElement*>(Mg_problem_pt->
810  mesh_pt()->element_pt(0))->dim()!=DIM)
811  {
812  std::string err_strng="The dimension of the elements used in the mesh ";
813  err_strng+="does not match the dimension of the solver.";
814  throw OomphLibError(err_strng,
815  OOMPH_CURRENT_FUNCTION,
816  OOMPH_EXCEPTION_LOCATION);
817  }
818 
819  // PARANOID check - The elements of the bulk mesh must all be refineable
820  // elements otherwise we cannot deal with this
821  if (Mg_problem_pt->mg_bulk_mesh_pt()!=0)
822  {
823  // Find the number of elements in the bulk mesh
824  unsigned n_elements=Mg_problem_pt->mg_bulk_mesh_pt()->nelement();
825 
826  // Loop over the elements in the mesh and ensure that they are
827  // all refineable elements
828  for (unsigned el_counter=0;el_counter<n_elements;el_counter++)
829  {
830  // Upcast global mesh element to a refineable element
831  RefineableElement* el_pt=dynamic_cast<RefineableElement*>
832  (Mg_problem_pt->mg_bulk_mesh_pt()->element_pt(el_counter));
833 
834  // Check if the upcast worked or not; if el_pt is a null pointer the
835  // element is not refineable
836  if (el_pt==0)
837  {
838  throw OomphLibError
839  ("Element in global mesh could not be upcast to a refineable "
840  "element. We cannot deal with elements that are not refineable.",
841  OOMPH_CURRENT_FUNCTION,OOMPH_EXCEPTION_LOCATION);
842  }
843  }
844  }
845  else
846  {
847  throw OomphLibError
848  ("The provided bulk mesh pointer is set to be a null pointer. "
849  "The multigrid solver operates on the bulk mesh thus a pointer "
850  "to the correct mesh must be given.",
851  OOMPH_CURRENT_FUNCTION,
852  OOMPH_EXCEPTION_LOCATION);
853  }
854 #endif
855 
856  // If this is not the first Newton step then we will already have things
857  // in storage. If this is the case, delete them
858  clean_up_memory();
859 
860  // Resize the Mg_hierarchy vector
861  Mg_hierarchy.resize(1,0);
862 
863  // Set the pointer to the finest level as the first entry in Mg_hierarchy
864  Mg_hierarchy[0]=Mg_problem_pt;
865 
866  // Create the hierarchy of levels
867  setup_mg_hierarchy();
868 
869  // Set up the interpolation and restriction matrices
870  setup_transfer_matrices();
871 
872  // Set up the data structures on each level, i.e. the system matrix,
873  // LHS and RHS vectors and smoothers
874  setup_mg_structures();
875 
876  // Set up the smoothers on all of the levels
877  setup_smoothers();
878 
879  // If we do not want to document everything we want to delete all the
880  // coarse-grid problems
881  if (!Doc_everything)
882  {
883  // Loop over all of the coarser levels
884  for (unsigned i=1;i<Nlevel;i++)
885  {
886  // Delete the i-th coarse-grid MGProblem
887  delete Mg_hierarchy[i];
888 
889  // Set it to be a null pointer
890  Mg_hierarchy[i]=0;
891  }
892  }
893  // Otherwise, document everything!
894  else
895  {
896  // If the user wishes to document everything we run the self-test
897  self_test();
898  } // if (!Doc_everything)
899 
900  // Indicate that the full setup has been completed
901  Has_been_setup=true;
902 
903  // If we're allowed to output to the screen
904  if (!Suppress_all_output)
905  {
906  // Output the time taken to complete the full setup
907  double t_fs_end=TimingHelpers::timer();
908  double full_setup_time=t_fs_end-t_fs_start;
909 
910  // Output the CPU time
911  oomph_info << "\nCPU time for full setup [sec]: "
912  << full_setup_time << std::endl;
913 
914  // Notify user that the hierarchy of levels is complete
915  oomph_info << "\n===============Multigrid Full Setup Complete=============="
916  << std::endl;
917  } // if (!Suppress_all_output)
918  } // End of full_setup
919 
920  //===================================================================
921  /// \short Set up the MG hierarchy
922  //===================================================================
923  // Function to set up the hierachy of levels. Creates a vector of
924  // pointers to each MG level
925  template<unsigned DIM>
927  {
928  // Initialise the timer start variable
929  double t_m_start=0.0;
930 
931  // Notify the user if it is allowed
932  if (!Suppress_all_output)
933  {
934  // Notify user of progress
935  oomph_info << "\n===============Creating Multigrid Hierarchy==============="
936  << std::endl;
937 
938  // Start clock
939  t_m_start=TimingHelpers::timer();
940  }
941 
942  // Create a bool to indicate whether or not we could create an unrefined
943  // copy. This bool will be assigned the value FALSE when the current copy
944  // is the last level of the multigrid hierarchy
945  bool managed_to_create_unrefined_copy=true;
946 
947  // Now keep making copies and try to make an unrefined copy of
948  // the mesh
949  unsigned level=0;
950 
951  // Set up all of the levels by making a completely unrefined copy
952  // of the problem using the function make_new_problem
953  while (managed_to_create_unrefined_copy)
954  {
955  // Make a new object of the same type as the derived problem
956  MGProblem* new_problem_pt=Mg_problem_pt->make_new_problem();
957 
958  // Do anything that needs to be done before we can refine the mesh
959  new_problem_pt->actions_before_adapt();
960 
961  // To create the next level in the hierarchy we need to create a mesh
962  // which matches the refinement of the current problem and then unrefine
963  // the mesh. This can alternatively be done using the function
964  // refine_base_mesh_as_in_reference_mesh_minus_one which takes a
965  // reference mesh to do precisely the above with
966  managed_to_create_unrefined_copy=
967  new_problem_pt->mg_bulk_mesh_pt()->
968  refine_base_mesh_as_in_reference_mesh_minus_one(
969  Mg_hierarchy[level]->mg_bulk_mesh_pt());
970 
971  // If we were able to unrefine the problem on the current level
972  // then add the unrefined problem to a vector of the levels
973  if (managed_to_create_unrefined_copy)
974  {
975  // Another level has been created so increment the level counter
976  level++;
977 
978  // If the documentation of everything has not been suppressed
979  // then tell the user we managed to create another level
980  if (!Suppress_all_output)
981  {
982  // Notify user that unrefinement was successful
983  oomph_info << "\nSuccess! Level " << level
984  << " has been created." << std::endl;
985  }
986 
987  // Do anything that needs to be done after refinement
988  new_problem_pt->actions_after_adapt();
989 
990  // Do the equation numbering for the new problem
991  oomph_info << "\n - Number of equations: "
992  << new_problem_pt->assign_eqn_numbers() << std::endl;
993 
994  // Add the new problem pointer onto the vector of MG levels
995  // and increment the value of level by 1
996  Mg_hierarchy.push_back(new_problem_pt);
997  }
998  // If we weren't able to create an unrefined copy
999  else
1000  {
1001  // Delete the new problem
1002  delete new_problem_pt;
1003 
1004  // Make it a null pointer
1005  new_problem_pt=0;
1006 
1007  // Assign the number of levels to Nlevel
1008  Nlevel=Mg_hierarchy.size();
1009 
1010  // If we're allowed to document then tell the user we've reached
1011  // the coarsest level of the hierarchy
1012  if (!Suppress_all_output)
1013  {
1014  // Notify the user
1015  oomph_info << "\n Reached the coarsest level! "
1016  << "Number of levels: " << Nlevel << std::endl;
1017  }
1018  } // if (managed_to_unrefine)
1019  } // while (managed_to_unrefine)
1020 
1021  // Given that we know the number of levels in the hierarchy we can resize the
1022  // vectors which will store all the information required for our solver:
1023  // Resize the vector storing all of the system matrices
1024  Mg_matrices_storage_pt.resize(Nlevel,0);
1025 
1026  // Resize the vector storing all of the solution vectors (X_mg)
1027  X_mg_vectors_storage.resize(Nlevel);
1028 
1029  // Resize the vector storing all of the RHS vectors (Rhs_mg)
1030  Rhs_mg_vectors_storage.resize(Nlevel);
1031 
1032  // Resize the vector storing all of the residual vectors
1033  Residual_mg_vectors_storage.resize(Nlevel);
1034 
1035  // Allocate space for the pre-smoother storage vector (remember, we do
1036  // not need a smoother on the coarsest level we use a direct solve there)
1037  Pre_smoothers_storage_pt.resize(Nlevel-1,0);
1038 
1039  // Allocate space for the post-smoother storage vector (remember, we do
1040  // not need a smoother on the coarsest level we use a direct solve there)
1041  Post_smoothers_storage_pt.resize(Nlevel-1,0);
1042 
1043  // Resize the vector storing all of the interpolation matrices
1044  Interpolation_matrices_storage_pt.resize(Nlevel-1,0);
1045 
1046  // Resize the vector storing all of the restriction matrices
1047  Restriction_matrices_storage_pt.resize(Nlevel-1,0);
1048 
1049  if (!Suppress_all_output)
1050  {
1051  // Stop clock
1052  double t_m_end=TimingHelpers::timer();
1053  double total_setup_time=double(t_m_end-t_m_start);
1054  oomph_info << "\nCPU time for creation of hierarchy of MG problems [sec]: "
1055  << total_setup_time << std::endl;
1056 
1057  // Notify user that the hierarchy of levels is complete
1058  oomph_info << "\n===============Hierarchy Creation Complete================"
1059  << "\n" << std::endl;
1060  }
1061  } // End of setup_mg_hierarchy
1062 
1063  //===================================================================
1064  /// \short Set up the transfer matrices. Both the pure injection and
1065  /// full weighting method have been implemented here but it is highly
1066  /// recommended that full weighting is used in general. In both
1067  /// methods the transpose of the transfer matrix is used to transfer
1068  /// a vector back
1069  //===================================================================
1070  template<unsigned DIM>
1072  {
1073  // Initialise the timer start variable
1074  double t_r_start=0.0;
1075 
1076  // Notify the user (if we're allowed)
1077  if (!Suppress_all_output)
1078  {
1079  // Notify user of progress
1080  oomph_info << "Creating the transfer matrices ";
1081 
1082  // Start the clock!
1083  t_r_start=TimingHelpers::timer();
1084  }
1085 
1086  // If we're allowed to output information
1087  if (!Suppress_all_output)
1088  {
1089  // Say what method we're using
1090  oomph_info << "using full weighting (recommended).\n" << std::endl;
1091  }
1092 
1093  // Using full weighting so use setup_interpolation_matrices.
1094  // Note: There are two methods to choose from here, the ideal choice is
1095  // setup_interpolation_matrices() but that requires a refineable mesh base
1096  if (dynamic_cast<TreeBasedRefineableMeshBase*>
1097  (Mg_problem_pt->mg_bulk_mesh_pt()))
1098  {
1099  setup_interpolation_matrices();
1100  }
1101  // If the mesh is unstructured we have to use the locate_zeta function
1102  // to set up the interpolation matrices
1103  else
1104  {
1105  setup_interpolation_matrices_unstructured();
1106  }
1107 
1108  // Loop over all levels that will be assigned a restriction matrix
1109  set_restriction_matrices_as_interpolation_transposes();
1110 
1111  // If we're allowed
1112  if (!Suppress_all_output)
1113  {
1114  // Stop the clock
1115  double t_r_end=TimingHelpers::timer();
1116  double total_G_setup_time=double(t_r_end-t_r_start);
1117  oomph_info << "CPU time for transfer matrices setup [sec]: "
1118  << total_G_setup_time << std::endl;
1119 
1120  // Notify user that the hierarchy of levels is complete
1121  oomph_info << "\n============Transfer Matrices Setup Complete=============="
1122  << "\n" << std::endl;
1123  }
1124  } // End of setup_transfer_matrices function
1125 
1126  //===================================================================
1127  /// \short Set up the MG hierarchy structures
1128  //===================================================================
1129  // Function to set up the hierachy of levels. Creates a vector of
1130  // pointers to each MG level
1131  template<unsigned DIM>
1133  {
1134  // Initialise the timer start variable
1135  double t_m_start=0.0;
1136 
1137  // Start the clock (if we're allowed to time things)
1138  if (!Suppress_all_output)
1139  {
1140  // Start the clock
1141  t_m_start=TimingHelpers::timer();
1142  }
1143 
1144  // Allocate space for the system matrix on each level
1145  for (unsigned i=0;i<Nlevel;i++)
1146  {
1147  // Dynamically allocate a new CRDoubleMatrix
1148  Mg_matrices_storage_pt[i]=new CRDoubleMatrix;
1149  }
1150 
1151  // Loop over each level and extract the system matrix, solution vector
1152  // right-hand side vector and residual vector (to store the value of r=b-Ax)
1153  for (unsigned i=0;i<Nlevel;i++)
1154  {
1155  // If we're allowed to output
1156  if (!Suppress_all_output)
1157  {
1158  // Output the level we're working on
1159  oomph_info << "Setting up MG structures on level: " << i
1160  << "\n" << std::endl;
1161  }
1162 
1163  // Resize the solution and RHS vector
1164  unsigned n_dof=Mg_hierarchy[i]->ndof();
1165  LinearAlgebraDistribution* dist_pt=
1166  new LinearAlgebraDistribution(Mg_hierarchy[i]->communicator_pt(),
1167  n_dof,false);
1168 
1169 #ifdef PARANOID
1170 #ifdef OOMPH_HAS_MPI
1171  // Set up the warning messages
1172  std::string warning_message="Setup of distribution has not been ";
1173  warning_message+="tested with MPI.";
1174 
1175  // If we're not running the code in serial
1176  if (dist_pt->communicator_pt()->nproc()>1)
1177  {
1178  // Throw a warning
1179  OomphLibWarning(warning_message,
1180  OOMPH_CURRENT_FUNCTION,
1181  OOMPH_EXCEPTION_LOCATION);
1182  }
1183 #endif
1184 #endif
1185 
1186  // Build the approximate solution
1187  X_mg_vectors_storage[i].clear();
1188  X_mg_vectors_storage[i].build(dist_pt);
1189 
1190  // Build the point source function
1191  Rhs_mg_vectors_storage[i].clear();
1192  Rhs_mg_vectors_storage[i].build(dist_pt);
1193 
1194  // Build the residual vector (r=b-Ax)
1195  Residual_mg_vectors_storage[i].clear();
1196  Residual_mg_vectors_storage[i].build(dist_pt);
1197 
1198  // Delete the distribution pointer
1199  delete dist_pt;
1200 
1201  // Make it a null pointer
1202  dist_pt=0;
1203 
1204  // Build the matrix distribution
1205  Mg_matrices_storage_pt[i]->clear();
1206  Mg_matrices_storage_pt[i]->distribution_pt()->
1207  build(Mg_hierarchy[i]->communicator_pt(),n_dof,false);
1208 
1209  // Compute system matrix on the current level. On the finest level of the
1210  // hierarchy the system matrix and RHS vector is given by the Jacobian and
1211  // vector of residuals which define the original problem which the Galerkin
1212  // approximation to the system matrix is used on the subsequent levels so
1213  // that the correct contributions are taken from each dof on the level above
1214  // (that is to say, it should match the contribution taken from the solution
1215  // vector and RHS vector on the level above)
1216  if (i==0)
1217  {
1218  // Initialise the timer start variable
1219  double t_jac_start=0.0;
1220 
1221  // If we're allowed to output things
1222  if (!Suppress_all_output)
1223  {
1224  // Start timer for Jacobian setup
1225  t_jac_start=TimingHelpers::timer();
1226  }
1227 
1228  // The system matrix on the finest level is the Jacobian and the RHS
1229  // vector is given by the residual vector which accompanies the Jacobian
1230  Mg_hierarchy[0]->get_jacobian(Rhs_mg_vectors_storage[0],
1231  *Mg_matrices_storage_pt[0]);
1232 
1233  if (!Suppress_all_output)
1234  {
1235  // Document the time taken
1236  double t_jac_end=TimingHelpers::timer();
1237  double jacobian_setup_time=t_jac_end-t_jac_start;
1238  oomph_info << " - Time for setup of Jacobian [sec]: "
1239  << jacobian_setup_time << "\n" << std::endl;
1240  }
1241  }
1242  else
1243  {
1244  // Initialise the timer start variable
1245  double t_gal_start=0.0;
1246 
1247  // If we're allowed
1248  if (!Suppress_all_output)
1249  {
1250  // Start timer for Galerkin matrix calculation
1251  t_gal_start=TimingHelpers::timer();
1252  }
1253 
1254  // The system matrix on the coarser levels must be formed using the
1255  // Galerkin approximation which we do by calculating the product
1256  // A^2h = I^2h_h * A^h * I^h_2h, i.e. the coarser version of the
1257  // finer grid system matrix is formed by multiplying by the (fine grid)
1258  // restriction matrix from the left and the (fine grid) interpolation
1259  // matrix from the left
1260 
1261  // First we need to calculate A^h * I^h_2h which we store as A^2h
1262  Mg_matrices_storage_pt[i-1]->
1263  multiply(*Interpolation_matrices_storage_pt[i-1],
1264  *Mg_matrices_storage_pt[i]);
1265 
1266  // Now calculate I^2h_h * (A^h * I^h_2h) where the quantity in brackets
1267  // was just calculated. This updates A^2h to give us the true
1268  // Galerkin approximation to the finer grid matrix
1269  Restriction_matrices_storage_pt[i-1]->multiply(*Mg_matrices_storage_pt[i],
1270  *Mg_matrices_storage_pt[i]);
1271 
1272  // If the user did not choose to suppress everything
1273  if (!Suppress_all_output)
1274  {
1275  // End timer for Galerkin matrix calculation
1276  double t_gal_end=TimingHelpers::timer();
1277 
1278  // Calculate setup time
1279  double galerkin_matrix_calculation_time=t_gal_end-t_gal_start;
1280 
1281  // Document the time taken
1282  oomph_info << " - Time for system matrix formation using the Galerkin "
1283  << "approximation [sec]: " << galerkin_matrix_calculation_time
1284  << "\n" << std::endl;
1285  }
1286  } // if (i==0) else
1287  } // for (unsigned i=0;i<Nlevel;i++)
1288 
1289  // If we're allowed to output
1290  if (!Suppress_all_output)
1291  {
1292  // Stop clock
1293  double t_m_end=TimingHelpers::timer();
1294  double total_setup_time=double(t_m_end-t_m_start);
1295  oomph_info << "Total CPU time for setup of MG structures [sec]: "
1296  << total_setup_time << std::endl;
1297 
1298  // Notify user that the hierarchy of levels is complete
1299  oomph_info << "\n============"
1300  << "Multigrid Structures Setup Complete"
1301  << "===========\n" << std::endl;
1302  }
1303  } // End of setup_mg_structures
1304 
1305 
1306 //=========================================================================
1307 /// \short Set up the smoothers on all levels
1308 //=========================================================================
1309  template<unsigned DIM>
1311  {
1312  // Initialise the timer start variable
1313  double t_m_start=0.0;
1314 
1315  // Start the clock (if we're allowed to time things)
1316  if (!Suppress_all_output)
1317  {
1318  // Notify user
1319  oomph_info << "Starting the setup of all smoothers.\n" << std::endl;
1320 
1321  // Start the clock
1322  t_m_start=TimingHelpers::timer();
1323  }
1324 
1325  // Loop over the levels and assign the pre- and post-smoother points
1326  for (unsigned i=0;i<Nlevel-1;i++)
1327  {
1328  // If the pre-smoother factory function pointer hasn't been assigned
1329  // then we simply create a new instance of the DampedJacobi smoother
1330  // which is the default pre-smoother
1331  if (0==Pre_smoother_factory_function_pt)
1332  {
1333  Pre_smoothers_storage_pt[i]=new DampedJacobi<CRDoubleMatrix>;
1334  }
1335  // Otherwise we use the pre-smoother factory function pointer to
1336  // generate a new pre-smoother
1337  else
1338  {
1339  // Get a pointer to an object of the same type as the pre-smoother
1340  Pre_smoothers_storage_pt[i]=(*Pre_smoother_factory_function_pt)();
1341  }
1342 
1343  // If the post-smoother factory function pointer hasn't been assigned
1344  // then we simply create a new instance of the DampedJacobi smoother
1345  // which is the default post-smoother
1346  if (0==Post_smoother_factory_function_pt)
1347  {
1348  Post_smoothers_storage_pt[i]=new DampedJacobi<CRDoubleMatrix>;
1349  }
1350  // Otherwise we use the post-smoother factory function pointer to
1351  // generate a new post-smoother
1352  else
1353  {
1354  // Get a pointer to an object of the same type as the post-smoother
1355  Post_smoothers_storage_pt[i]=(*Post_smoother_factory_function_pt)();
1356  }
1357  }
1358 
1359  // Set the tolerance for the pre- and post-smoothers. The norm of the
1360  // solution which is compared against the tolerance is calculated
1361  // differently to the multigrid solver. To ensure that the smoother
1362  // continues to solve for the prescribed number of iterations we
1363  // lower the tolerance
1364  for (unsigned i=0;i<Nlevel-1;i++)
1365  {
1366  // Set the tolerance of the i-th level pre-smoother
1367  Pre_smoothers_storage_pt[i]->tolerance()=1.0e-16;
1368 
1369  // Set the tolerance of the i-th level post-smoother
1370  Post_smoothers_storage_pt[i]->tolerance()=1.0e-16;
1371  }
1372 
1373  // Set the number of pre- and post-smoothing iterations in each
1374  // pre- and post-smoother
1375  for (unsigned i=0;i<Nlevel-1;i++)
1376  {
1377  // Set the number of pre-smoothing iterations as the value of Npre_smooth
1378  Pre_smoothers_storage_pt[i]->max_iter()=Npre_smooth;
1379 
1380  // Set the number of pre-smoothing iterations as the value of Npost_smooth
1381  Post_smoothers_storage_pt[i]->max_iter()=Npost_smooth;
1382  }
1383 
1384  // Complete the setup of all of the smoothers
1385  for (unsigned i=0;i<Nlevel-1;i++)
1386  {
1387  // Pass a pointer to the system matrix on the i-th level to the i-th
1388  // level pre-smoother
1389  Pre_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
1390 
1391  // Pass a pointer to the system matrix on the i-th level to the i-th
1392  // level post-smoother
1393  Post_smoothers_storage_pt[i]->smoother_setup(Mg_matrices_storage_pt[i]);
1394  }
1395 
1396  // Set up the distributions of each smoother
1397  for (unsigned i=0;i<Nlevel-1;i++)
1398  {
1399  // Get the number of dofs on the i-th level of the hierarchy.
1400  // This will be the same as the size of the solution vector
1401  // associated with the i-th level
1402  unsigned n_dof=X_mg_vectors_storage[i].nrow();
1403 
1404  // Create a LinearAlgebraDistribution which will be passed to the
1405  // linear solver
1406  LinearAlgebraDistribution dist(Mg_hierarchy[i]->communicator_pt(),
1407  n_dof,false);
1408 
1409 #ifdef PARANOID
1410 #ifdef OOMPH_HAS_MPI
1411  // Set up the warning messages
1412  std::string warning_message="Setup of pre- and post-smoother distribution ";
1413  warning_message+="has not been tested with MPI.";
1414 
1415  // If we're not running the code in serial
1416  if (dist.communicator_pt()->nproc()>1)
1417  {
1418  // Throw a warning
1419  OomphLibWarning(warning_message,
1420  OOMPH_CURRENT_FUNCTION,
1421  OOMPH_EXCEPTION_LOCATION);
1422  }
1423 #endif
1424 #endif
1425 
1426  // Build the distribution of the pre-smoother
1427  Pre_smoothers_storage_pt[i]->build_distribution(dist);
1428 
1429  // Build the distribution of the post-smoother
1430  Post_smoothers_storage_pt[i]->build_distribution(dist);
1431  }
1432 
1433  // Disable the smoother output on this level
1434  if (!Doc_time)
1435  {
1436  disable_smoother_and_superlu_doc_time();
1437  }
1438 
1439  // If we're allowed to output
1440  if (!Suppress_all_output)
1441  {
1442  // Stop clock
1443  double t_m_end=TimingHelpers::timer();
1444  double total_setup_time=double(t_m_end-t_m_start);
1445  oomph_info << "CPU time for setup of smoothers on all levels [sec]: "
1446  << total_setup_time << std::endl;
1447 
1448  // Notify user that the extraction is complete
1449  oomph_info << "\n==================Smoother Setup Complete================="
1450  << std::endl;
1451  }
1452  } // End of setup_smoothers
1453 
1454 
1455 //===================================================================
1456 /// \short Setup the interpolation matrices
1457 //===================================================================
1458  template<unsigned DIM>
1460  {
1461  // Variable to hold the number of sons in an element
1462  unsigned n_sons;
1463 
1464  // Number of son elements
1465  if (DIM==2)
1466  {
1467  n_sons=4;
1468  }
1469  else if (DIM==3)
1470  {
1471  n_sons=8;
1472  }
1473  else
1474  {
1475  std::ostringstream error_message_stream;
1476  error_message_stream << "DIM should be 2 or 3 not "
1477  << DIM << std::endl;
1478  throw OomphLibError(error_message_stream.str(),
1479  OOMPH_CURRENT_FUNCTION,
1480  OOMPH_EXCEPTION_LOCATION);
1481  }
1482 
1483  // Vector of local coordinates in the element
1484  Vector<double> s(DIM,0.0);
1485 
1486  // Loop over each level (apart from the coarsest level; an interpolation
1487  // matrix and thus a restriction matrix is not needed here), with 0 being
1488  // the finest level and Nlevel-1 being the coarsest level
1489  for (unsigned level=0;level<Nlevel-1;level++)
1490  {
1491  // Assign values to a couple of variables to help out
1492  unsigned coarse_level=level+1;
1493  unsigned fine_level=level;
1494 
1495  // Make a pointer to the mesh on the finer level and dynamic_cast
1496  // it as an object of the refineable mesh class
1497  RefineableMeshBase* ref_fine_mesh_pt=
1498  dynamic_cast<RefineableMeshBase*>
1499  (Mg_hierarchy[fine_level]->mg_bulk_mesh_pt());
1500 
1501  // Make a pointer to the mesh on the coarse level and dynamic_cast
1502  // it as an object of the refineable mesh class
1503  RefineableMeshBase* ref_coarse_mesh_pt=
1504  dynamic_cast<RefineableMeshBase*>
1505  (Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
1506 
1507  // Access information about the number of elements in the fine mesh
1508  // from the pointer to the fine mesh (to loop over the rows of the
1509  // interpolation matrix)
1510  unsigned fine_n_element=ref_fine_mesh_pt->nelement();
1511 
1512  // The numbers of rows and columns in the interpolation matrix. The
1513  // number of unknowns has been divided by 2 since there are 2 dofs at
1514  // each node in the mesh corresponding to the real and imaginary part
1515  // of the solution
1516  unsigned n_rows=Mg_hierarchy[fine_level]->ndof();
1517  unsigned n_cols=Mg_hierarchy[coarse_level]->ndof();
1518 
1519  // Mapping relating the pointers to related elements in the coarse and
1520  // fine meshes: coarse_mesh_element_pt[fine_mesh_element_pt]
1521  // If the element in the fine mesh has been unrefined between these two
1522  // levels, this map returns the father element in the coarsened mesh.
1523  // If this element in the fine mesh has not been unrefined,
1524  // the map returns the pointer to the same-sized equivalent
1525  // element in the coarsened mesh.
1526  std::map<RefineableQElement<DIM>*,RefineableQElement<DIM>*>
1527  coarse_mesh_reference_element_pt;
1528 
1529  // Counter of elements in coarse and fine meshes: Start with element
1530  // zero in both meshes.
1531  unsigned e_coarse=0;
1532  unsigned e_fine=0;
1533 
1534  // While loop over fine elements (while because we're
1535  // incrementing the counter internally if the element was
1536  // unrefined...)
1537  while(e_fine<fine_n_element)
1538  {
1539  // Pointer to element in fine mesh
1540  RefineableQElement<DIM>* el_fine_pt=
1541  dynamic_cast<RefineableQElement<DIM>*>
1542  (ref_fine_mesh_pt->finite_element_pt(e_fine));
1543 
1544  // Pointer to element in coarse mesh
1545  RefineableQElement<DIM>* el_coarse_pt=
1546  dynamic_cast<RefineableQElement<DIM>*>
1547  (ref_coarse_mesh_pt->finite_element_pt(e_coarse));
1548 
1549  // If the levels are different then the element in the fine
1550  // mesh has been unrefined between these two levels
1551  if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level())
1552  {
1553  // The element in the fine mesh has been unrefined between these two
1554  // levels. Hence it and its three brothers (ASSUMED to be stored
1555  // consecutively in the Mesh's vector of pointers to its constituent
1556  // elements -- we'll check this!) share the same father element in
1557  // the coarse mesh, currently pointed to by el_coarse_pt.
1558  for(unsigned i=0;i<n_sons;i++)
1559  {
1560  // Set mapping to father element in coarse mesh
1561  coarse_mesh_reference_element_pt[
1562  dynamic_cast<RefineableQElement<DIM>*>
1563  (ref_fine_mesh_pt->finite_element_pt(e_fine))]=el_coarse_pt;
1564 
1565  // Increment counter for elements in fine mesh
1566  e_fine++;
1567  }
1568  }
1569  // The element in the fine mesh has not been unrefined between
1570  // these two levels, so the reference element is the same-sized
1571  // equivalent element in the coarse mesh
1572  else
1573  {
1574  // Set the mapping between the two elements since they are
1575  // the same element
1576  coarse_mesh_reference_element_pt[el_fine_pt]=el_coarse_pt;
1577 
1578  // Increment counter for elements in fine mesh
1579  e_fine++;
1580  }
1581  // Increment counter for elements in coarse mesh
1582  e_coarse++;
1583  } // End of while loop for setting up element-coincidences
1584 
1585  // To allow update of a row only once we use stl vectors for bools
1586  std::vector<bool> contribution_made(n_rows,false);
1587 
1588  // Make storage vectors to form the interpolation matrix using a
1589  // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
1590  // row_start contains entries which tells us at which entry of the
1591  // vector column_start we start on the i-th row of the actual matrix.
1592  // The entries of column_start indicate the column position of each
1593  // non-zero entry in every row. This runs through the first row
1594  // from left to right then the second row (again, left to right)
1595  // and so on until the end. The entries in value are the entries in
1596  // the interpolation matrix. The vector column_start has the same length
1597  // as the vector value.
1598  Vector<double> value;
1599  Vector<int> column_index;
1600  Vector<int> row_start(n_rows+1);
1601 
1602  // The value of index will tell us which row of the interpolation matrix
1603  // we're working on in the following for loop
1604  unsigned index=0;
1605 
1606  // New loop to go over each element in the fine mesh
1607  for (unsigned k=0;k<fine_n_element;k++)
1608  {
1609  // Set a pointer to the element in the fine mesh
1610  RefineableQElement<DIM>* el_fine_pt=
1611  dynamic_cast<RefineableQElement<DIM>*>
1612  (ref_fine_mesh_pt->finite_element_pt(k));
1613 
1614  // Get the reference element (either the father element or the
1615  // same-sized element) in the coarse mesh
1616  RefineableQElement<DIM>* el_coarse_pt=
1617  coarse_mesh_reference_element_pt[el_fine_pt];
1618 
1619  // Find out what type of son it is (set to OMEGA if no unrefinement
1620  // took place)
1621  int son_type=0;
1622 
1623  // Check if the elements are different on both levels (i.e. to check
1624  // if any unrefinement took place)
1625  if (el_fine_pt->tree_pt()->level()!=el_coarse_pt->tree_pt()->level())
1626  {
1627  // If there was refinement we need to find the son type
1628  son_type=el_fine_pt->tree_pt()->son_type();
1629  }
1630  else
1631  {
1632  // If there was no refinement then the son_type is given by the
1633  // value of Tree::OMEGA
1634  son_type=Tree::OMEGA;
1635  }
1636 
1637  // Find the number of nodes in the fine mesh element
1638  unsigned nnod_fine=el_fine_pt->nnode();
1639 
1640  // Loop through all the nodes in an element in the fine mesh
1641  for(unsigned i=0;i<nnod_fine;i++)
1642  {
1643  // Row number in interpolation matrix: Global equation number
1644  // of the d.o.f. stored at this node in the fine element
1645  int ii=el_fine_pt->node_pt(i)->eqn_number(0);
1646 
1647  // Check whether or not the node is a proper d.o.f.
1648  if (ii>=0)
1649  {
1650  // Only assign values to the given row of the interpolation
1651  // matrix if they haven't already been assigned
1652  if(contribution_made[ii]==false)
1653  {
1654  // The value of index was initialised when we allocated space
1655  // for the three vectors to store the CRDoubleMatrix. We use
1656  // index to go through the entries of the row_start vector.
1657  // The next row starts at the value.size() (draw out the entries
1658  // in value if this doesn't make sense noting that the storage
1659  // for the vector 'value' is dynamically allocated)
1660  row_start[index]=value.size();
1661 
1662  // Calculate the local coordinates of the given node
1663  el_fine_pt->local_coordinate_of_node(i,s);
1664 
1665  // Find the local coordinates s, of the present node, in the
1666  // reference element in the coarse mesh, given the element's
1667  // son_type (e.g. SW,SE... )
1668  level_up_local_coord_of_node(son_type,s);
1669 
1670  // Allocate space for shape functions in the coarse mesh
1671  Shape psi(el_coarse_pt->nnode());
1672 
1673  // Set the shape function in the reference element
1674  el_coarse_pt->shape(s,psi);
1675 
1676  // Auxiliary storage
1677  std::map<unsigned,double> contribution;
1678  Vector<unsigned> keys;
1679 
1680  // Find the number of nodes in an element in the coarse mesh
1681  unsigned nnod_coarse=el_coarse_pt->nnode();
1682 
1683  // Loop through all the nodes in the reference element
1684  for (unsigned j=0;j<nnod_coarse;j++)
1685  {
1686  // Column number in interpolation matrix: Global equation
1687  // number of the d.o.f. stored at this node in the coarse
1688  // element
1689  int jj=el_coarse_pt->node_pt(j)->eqn_number(0);
1690 
1691  // If the value stored at this node is pinned or hanging
1692  if (jj<0)
1693  {
1694  // Hanging node: In this case we need to accumulate the
1695  // contributions from the master nodes
1696  if (el_coarse_pt->node_pt(j)->is_hanging())
1697  {
1698  // Find the number of master nodes of the hanging
1699  // the node in the reference element
1700  HangInfo* hang_info_pt=el_coarse_pt->node_pt(j)->hanging_pt();
1701  unsigned nmaster=hang_info_pt->nmaster();
1702 
1703  // Loop over the master nodes
1704  for (unsigned i_master=0;i_master<nmaster;i_master++)
1705  {
1706  // Set up a pointer to the master node
1707  Node* master_node_pt=hang_info_pt->master_node_pt(i_master);
1708 
1709  // The column number in the interpolation matrix: the
1710  // global equation number of the d.o.f. stored at this master
1711  // node for the coarse element
1712  int master_jj=master_node_pt->eqn_number(0);
1713 
1714  // Is the master node a proper d.o.f.?
1715  if (master_jj>=0)
1716  {
1717  // If the weight of the master node is non-zero
1718  if (psi(j)*hang_info_pt->master_weight(i_master)!=0.0)
1719  {
1720  contribution[master_jj]+=
1721  psi(j)*hang_info_pt->master_weight(i_master);
1722  }
1723  } // if (master_jj>=0)
1724  } // for (unsigned i_master=0;i_master<nmaster;i_master++)
1725  } // if (el_coarse_pt->node_pt(j)->is_hanging())
1726  }
1727  // In the case that the node is not pinned or hanging
1728  else
1729  {
1730  // If we can get a nonzero contribution from the shape function
1731  if (psi(j)!=0.0)
1732  {
1733  contribution[jj]+=psi(j);
1734  }
1735  } // if (jj<0) else
1736  } // for (unsigned j=0;j<nnod_coarse;j++)
1737 
1738  // Put the contributions into the value vector
1739  for (std::map<unsigned,double>::iterator it=contribution.begin();
1740  it!=contribution.end(); ++it)
1741  {
1742  if (it->second!=0)
1743  {
1744  column_index.push_back(it->first);
1745  value.push_back(it->second);
1746  }
1747  } // for (std::map<unsigned,double>::iterator it=...)
1748 
1749  // Increment the value of index by 1 to indicate that we have
1750  // finished the row, or equivalently, covered another global
1751  // node in the fine mesh
1752  index++;
1753 
1754  // Change the entry in contribution_made to true now to indicate
1755  // that the row has been filled
1756  contribution_made[ii]=true;
1757  } // if(contribution_made[ii]==false)
1758  } // if (ii>=0)
1759  } // for(unsigned i=0;i<nnod_element;i++)
1760  } // for (unsigned k=0;k<fine_n_element;k++)
1761 
1762  // Set the last entry in the row_start vector
1763  row_start[n_rows]=value.size();
1764 
1765  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
1766  // using the vectors value, row_start, column_index and the value
1767  // of fine_n_unknowns and coarse_n_unknowns
1768  interpolation_matrix_set(level,
1769  value,
1770  column_index,
1771  row_start,
1772  n_cols,
1773  n_rows);
1774  } // for (unsigned level=0;level<Nlevel-1;level++)
1775  } // End of setup_interpolation_matrices
1776 
1777 //===================================================================
1778 /// \short Setup the interpolation matrices
1779 //===================================================================
1780  template<unsigned DIM>
1782  {
1783  // Vector of local coordinates in the element
1784  Vector<double> s(DIM,0.0);
1785 
1786  // Loop over each level (apart from the coarsest level; an interpolation
1787  // matrix and thus a restriction matrix is not needed here), with 0 being
1788  // the finest level and Nlevel-1 being the coarsest level
1789  for (unsigned level=0;level<Nlevel-1;level++)
1790  {
1791  // Assign values to a couple of variables to help out
1792  unsigned coarse_level=level+1;
1793  unsigned fine_level=level;
1794 
1795  // Make a pointer to the mesh on the finer level and dynamic_cast
1796  // it as an object of the refineable mesh class
1797  Mesh* ref_fine_mesh_pt=Mg_hierarchy[fine_level]->mg_bulk_mesh_pt();
1798 
1799  // To use the locate zeta functionality the coarse mesh must be of the
1800  // type MeshAsGeomObject
1801  MeshAsGeomObject* coarse_mesh_from_obj_pt=
1802  new MeshAsGeomObject(Mg_hierarchy[coarse_level]->mg_bulk_mesh_pt());
1803 
1804  // Access information about the number of degrees of freedom
1805  // from the pointers to the problem on each level
1806  unsigned coarse_n_unknowns=Mg_hierarchy[coarse_level]->ndof();
1807  unsigned fine_n_unknowns=Mg_hierarchy[fine_level]->ndof();
1808 
1809  // Make storage vectors to form the interpolation matrix using a
1810  // condensed row matrix (CRDoubleMatrix). The i-th value of the vector
1811  // row_start contains entries which tells us at which entry of the
1812  // vector column_start we start on the i-th row of the actual matrix.
1813  // The entries of column_start indicate the column position of each
1814  // non-zero entry in every row. This runs through the first row
1815  // from left to right then the second row (again, left to right)
1816  // and so on until the end. The entries in value are the entries in
1817  // the interpolation matrix. The vector column_start has the same length
1818  // as the vector value.
1819  Vector<double> value;
1820  Vector<int> column_index;
1821  Vector<int> row_start(fine_n_unknowns+1);
1822 
1823  // Vector to contain the (Eulerian) spatial location of the fine node
1824  Vector<double> fine_node_position(DIM);
1825 
1826  // Find the number of nodes in the mesh
1827  unsigned n_node_fine_mesh=ref_fine_mesh_pt->nnode();
1828 
1829  // Loop over the unknowns in the mesh
1830  for (unsigned i_fine_node=0;i_fine_node<n_node_fine_mesh;i_fine_node++)
1831  {
1832  // Set a pointer to the i_fine_unknown-th node in the fine mesh
1833  Node* fine_node_pt=ref_fine_mesh_pt->node_pt(i_fine_node);
1834 
1835  // Get the global equation number
1836  int i_fine=fine_node_pt->eqn_number(0);
1837 
1838  // If the node is a proper d.o.f.
1839  if (i_fine>=0)
1840  {
1841  // Row number in interpolation matrix: Global equation number
1842  // of the d.o.f. stored at this node in the fine element
1843  row_start[i_fine]=value.size();
1844 
1845  // Get the (Eulerian) spatial location of the fine node
1846  fine_node_pt->position(fine_node_position);
1847 
1848  // Create a null pointer to the GeomObject class
1849  GeomObject* el_pt=0;
1850 
1851  // Get the reference element (either the father element or the
1852  // same-sized element) in the coarse mesh using locate_zeta
1853  coarse_mesh_from_obj_pt->locate_zeta(fine_node_position,el_pt,s);
1854 
1855  // Upcast GeomElement as a FiniteElement
1856  FiniteElement* el_coarse_pt=dynamic_cast<FiniteElement*>(el_pt);
1857 
1858  // Find the number of nodes in the element
1859  unsigned n_node=el_coarse_pt->nnode();
1860 
1861  // Allocate space for shape functions in the coarse mesh
1862  Shape psi(n_node);
1863 
1864  // Calculate the geometric shape functions at local coordinate s
1865  el_coarse_pt->shape(s,psi);
1866 
1867  // Auxiliary storage
1868  std::map<unsigned,double> contribution;
1869  Vector<unsigned> keys;
1870 
1871  // Loop through all the nodes in the (coarse mesh) element containing the
1872  // node pointed to by fine_node_pt (fine mesh)
1873  for(unsigned j_node=0;j_node<n_node;j_node++)
1874  {
1875  // Get the j_coarse_unknown-th node in the coarse element
1876  Node* coarse_node_pt=el_coarse_pt->node_pt(j_node);
1877 
1878  // Column number in interpolation matrix: Global equation number of
1879  // the d.o.f. stored at this node in the coarse element
1880  int j_coarse=coarse_node_pt->eqn_number(0);
1881 
1882  // If the value stored at this node is pinned or hanging
1883  if (j_coarse<0)
1884  {
1885  // Hanging node: In this case we need to accumulate the
1886  // contributions from the master nodes
1887  if (el_coarse_pt->node_pt(j_node)->is_hanging())
1888  {
1889  // Find the number of master nodes of the hanging
1890  // the node in the reference element
1891  HangInfo* hang_info_pt=coarse_node_pt->hanging_pt();
1892  unsigned nmaster=hang_info_pt->nmaster();
1893 
1894  // Loop over the master nodes
1895  for (unsigned i_master=0;i_master<nmaster;i_master++)
1896  {
1897  // Set up a pointer to the master node
1898  Node* master_node_pt=hang_info_pt->master_node_pt(i_master);
1899 
1900  // The column number in the interpolation matrix: the
1901  // global equation number of the d.o.f. stored at this master
1902  // node for the coarse element
1903  int master_jj=master_node_pt->eqn_number(0);
1904 
1905  // Is the master node a proper d.o.f.?
1906  if (master_jj>=0)
1907  {
1908  // If the weight of the master node is non-zero
1909  if (psi(j_node)*hang_info_pt->master_weight(i_master)!=0.0)
1910  {
1911  contribution[master_jj]+=
1912  psi(j_node)*hang_info_pt->master_weight(i_master);
1913  }
1914  } // End of if statement (check it's not a boundary node)
1915  } // End of the loop over the master nodes
1916  } // End of the if statement for only hanging nodes
1917  } // End of the if statement for pinned or hanging nodes
1918  // In the case that the node is not pinned or hanging
1919  else
1920  {
1921  // If we can get a nonzero contribution from the shape function
1922  // at the j_node-th node in the element
1923  if (psi(j_node)!=0.0)
1924  {
1925  contribution[j_coarse]+=psi(j_node);
1926  }
1927  } // End of the if-else statement (check if the node was pinned/hanging)
1928  } // Finished loop over the nodes j in the reference element (coarse)
1929 
1930  // Put the contributions into the value vector
1931  for (std::map<unsigned,double>::iterator it=contribution.begin();
1932  it!=contribution.end(); ++it)
1933  {
1934  if (it->second!=0)
1935  {
1936  value.push_back(it->second);
1937  column_index.push_back(it->first);
1938  }
1939  } // End of putting contributions into the value vector
1940  } // End check (whether or not the fine node was a d.o.f.)
1941  } // End of the for-loop over nodes in the fine mesh
1942 
1943  // Set the last entry of row_start
1944  row_start[fine_n_unknowns]=value.size();
1945 
1946  // Set the interpolation matrix to be that formed as the CRDoubleMatrix
1947  // using the vectors value, row_start, column_index and the value
1948  // of fine_n_unknowns and coarse_n_unknowns
1949  interpolation_matrix_set(level,
1950  value,
1951  column_index,
1952  row_start,
1953  coarse_n_unknowns,
1954  fine_n_unknowns);
1955  } // End of loop over each level
1956  } // End of setup_interpolation_matrices_unstructured
1957 
1958  //=========================================================================
1959  /// \short Given the son type of the element and the local coordinate s of
1960  /// a given node in the son element, return the local coordinate s in its
1961  /// father element. 3D case.
1962  //=========================================================================
1963  template<>
1965  Vector<double>& s)
1966  {
1967  // If the element is unrefined between the levels the local coordinate
1968  // of the node in one element is the same as that in the other element
1969  // therefore we only need to perform calculations if the levels are
1970  // different (i.e. son_type is not OMEGA)
1971  if (son_type!=Tree::OMEGA)
1972  {
1973  // Scale the local coordinate from the range [-1,1]x[-1,1]x[-1,1]
1974  // to the range [0,1]x[0,1]x[0,1] to match the width of the local
1975  // coordinate range of the fine element from the perspective of
1976  // the father element. This then simply requires a shift of the
1977  // coordinates to match which type of son element we're dealing with
1978  s[0]=(s[0]+1.0)/2.0;
1979  s[1]=(s[1]+1.0)/2.0;
1980  s[2]=(s[2]+1.0)/2.0;
1981 
1982  // Cases: The son_type determines how the local coordinates should be
1983  // shifted to give the local coordinates in the coarse mesh element
1984  switch(son_type)
1985  {
1986  case OcTreeNames::LDF:
1987  s[0]-=1;
1988  s[1]-=1;
1989  break;
1990 
1991  case OcTreeNames::LDB:
1992  s[0]-=1;
1993  s[1]-=1;
1994  s[2]-=1;
1995  break;
1996 
1997  case OcTreeNames::LUF:
1998  s[0]-=1;
1999  break;
2000 
2001  case OcTreeNames::LUB:
2002  s[0]-=1;
2003  s[2]-=1;
2004  break;
2005 
2006  case OcTreeNames::RDF:
2007  s[1]-=1;
2008  break;
2009 
2010  case OcTreeNames::RDB:
2011  s[1]-=1;
2012  s[2]-=1;
2013  break;
2014 
2015  case OcTreeNames::RUF:
2016  break;
2017 
2018  case OcTreeNames::RUB:
2019  s[2]-=1;
2020  break;
2021  }
2022  } // if son_type!=Tree::OMEGA
2023  } // End of level_up_local_coord_of_node
2024 
2025  //=========================================================================
2026  /// \short Given the son type of the element and the local coordinate s of
2027  /// a given node in the son element, return the local coordinate s in its
2028  /// father element. 2D case.
2029  //=========================================================================
2030  template<>
2032  Vector<double>& s)
2033  {
2034  // If the element is unrefined between the levels the local coordinate
2035  // of the node in one element is the same as that in the other element
2036  // therefore we only need to perform calculations if the levels are
2037  // different (i.e. son_type is not OMEGA)
2038  if (son_type!=Tree::OMEGA)
2039  {
2040  // Scale the local coordinate from the range [-1,1]x[-1,1] to the range
2041  // [0,1]x[0,1] to match the width of the local coordinate range of the
2042  // fine element from the perspective of the father element. This
2043  // then simply requires a shift of the coordinates to match which type
2044  // of son element we're dealing with
2045  s[0]=(s[0]+1.0)/2.0;
2046  s[1]=(s[1]+1.0)/2.0;
2047 
2048  // Cases: The son_type determines how the local coordinates should be
2049  // shifted to give the local coordinates in the coarse mesh element
2050  switch(son_type)
2051  {
2052  // If we're dealing with the bottom-left element we need to shift
2053  // the range [0,1]x[0,1] to [-1,0]x[-1,0]
2054  case QuadTreeNames::SW:
2055  s[0]-=1;
2056  s[1]-=1;
2057  break;
2058 
2059  // If we're dealing with the bottom-right element we need to shift
2060  // the range [0,1]x[0,1] to [0,1]x[-1,0]
2061  case QuadTreeNames::SE:
2062  s[1]-=1;
2063  break;
2064 
2065  // If we're dealing with the top-right element we need to shift the
2066  // range [0,1]x[0,1] to [0,1]x[0,1], i.e. nothing needs to be done
2067  case QuadTreeNames::NE:
2068  break;
2069 
2070  // If we're dealing with the top-left element we need to shift
2071  // the range [0,1]x[0,1] to [-1,0]x[0,1]
2072  case QuadTreeNames::NW:
2073  s[0]-=1;
2074  break;
2075  }
2076  } // if son_type!=Tree::OMEGA
2077  } // End of level_up_local_coord_of_node
2078 
2079  //===================================================================
2080  /// \short Restrict residual (computed on current MG level) to
2081  /// next coarser mesh and stick it into the coarse mesh RHS vector
2082  /// using the restriction matrix (if restrict_flag=1) or the transpose
2083  /// of the interpolation matrix (if restrict_flag=2)
2084  //===================================================================
2085  template<unsigned DIM>
2086  void MGSolver<DIM>::restrict_residual(const unsigned& level)
2087  {
2088 #ifdef PARANOID
2089  // Check to make sure we can actually restrict the vector
2090  if (!(level<Nlevel-1))
2091  {
2092  throw OomphLibError("Input exceeds the possible parameter choice.",
2093  OOMPH_CURRENT_FUNCTION,
2094  OOMPH_EXCEPTION_LOCATION);
2095  }
2096 #endif
2097 
2098  // Multiply the residual vector by the restriction matrix on the level-th
2099  // level (to restrict the vector down to the next coarser level)
2100  Restriction_matrices_storage_pt[level]->
2101  multiply(Residual_mg_vectors_storage[level],
2102  Rhs_mg_vectors_storage[level+1]);
2103  } // End of restrict_residual
2104 
2105  //===================================================================
2106  /// \short Interpolate solution at current level onto
2107  /// next finer mesh and correct the solution x at that level
2108  //===================================================================
2109  template<unsigned DIM>
2110  void MGSolver<DIM>::interpolate_and_correct(const unsigned& level)
2111  {
2112 #ifdef PARANOID
2113  // Check to make sure we can actually restrict the vector
2114  if (!(level>0))
2115  {
2116  throw OomphLibError("Input level exceeds the possible parameter choice.",
2117  OOMPH_CURRENT_FUNCTION,
2118  OOMPH_EXCEPTION_LOCATION);
2119  }
2120 #endif
2121 
2122  // Build distribution of a temporary vector
2123  DoubleVector temp_soln(X_mg_vectors_storage[level-1].distribution_pt());
2124 
2125  // Interpolate the solution vector
2126  Interpolation_matrices_storage_pt[level-1]->
2127  multiply(X_mg_vectors_storage[level],temp_soln);
2128 
2129  // Update
2130  X_mg_vectors_storage[level-1]+=temp_soln;
2131  } // End of interpolate_and_correct
2132 
2133 //===================================================================
2134 /// \short Modify the restriction matrices
2135 //===================================================================
2136  template<unsigned DIM>
2138  {
2139  // Create a null pointer
2140  CRDoubleMatrix* restriction_matrix_pt=0;
2141 
2142  // Loop over the levels
2143  for (unsigned level=0;level<Nlevel-1;level++)
2144  {
2145  // Store a pointer to the (level)-th restriction matrix
2146  restriction_matrix_pt=Restriction_matrices_storage_pt[level];
2147 
2148  // Get access to the row start data
2149  const int* row_start_pt=restriction_matrix_pt->row_start();
2150 
2151  // Get access to the matrix entries
2152  double* value_pt=restriction_matrix_pt->value();
2153 
2154  // Initialise an auxiliary variable to store the index of the start
2155  // of the i-th row
2156  unsigned start_index=0;
2157 
2158  // Initialise an auxiliary variable to store the index of the start
2159  // of the (i+1)-th row
2160  unsigned end_index=0;
2161 
2162  // Store the number of rows in the matrix
2163  unsigned n_row=restriction_matrix_pt->nrow();
2164 
2165  // Loop over the rows of the matrix
2166  for (unsigned i=0;i<n_row;i++)
2167  {
2168  // Index associated with the start of the i-th row
2169  start_index=row_start_pt[i];
2170 
2171  // Index associated with the start of the (i+1)-th row
2172  end_index=row_start_pt[i+1];
2173 
2174  // Variable to store the sum of the absolute values of the off-diagonal
2175  // entries in the i-th row
2176  double row_sum=0.0;
2177 
2178  // Loop over the entries in the row
2179  for (unsigned j=start_index;j<end_index;j++)
2180  {
2181  // Add the absolute value of this entry to the off-diagonal sum
2182  row_sum+=value_pt[j];
2183  }
2184 
2185  // Loop over the entries in the row
2186  for (unsigned j=start_index;j<end_index;j++)
2187  {
2188  // Normalise the row entry
2189  value_pt[j]/=row_sum;
2190  }
2191  } // for (unsigned i=0;i<n_row;i++)
2192  } // for (unsigned level=0;level<Nlevel-1;level++)
2193  } // End of modify_restriction_matrices
2194 
2195  //===================================================================
2196  /// \short Makes a vector, restricts it down the levels of the hierarchy
2197  /// and documents it at each level. After this is done the vector is
2198  /// interpolated up the levels of the hierarchy with the output
2199  /// being documented at each level
2200  //===================================================================
2201  template<unsigned DIM>
2203  {
2204  // Start the timer
2205  double t_st_start=TimingHelpers::timer();
2206 
2207  // Notify user that the hierarchy of levels is complete
2208  oomph_info << "\nStarting the multigrid solver self-test."
2209  << std::endl;
2210 
2211  // Resize the vector storing all of the restricted vectors
2212  Restriction_self_test_vectors_storage.resize(Nlevel);
2213 
2214  // Resize the vector storing all of the interpolated vectors
2215  Interpolation_self_test_vectors_storage.resize(Nlevel);
2216 
2217  // Find the number of dofs on the finest level
2218  unsigned n_dof=X_mg_vectors_storage[0].nrow();
2219 
2220  // Need to set up the distribution of the finest level vector
2221  LinearAlgebraDistribution* dist_pt=
2222  new LinearAlgebraDistribution(Mg_problem_pt->communicator_pt(),n_dof,false);
2223 
2224 #ifdef PARANOID
2225 #ifdef OOMPH_HAS_MPI
2226  // Set up the warning messages
2227  std::string warning_message="Setup of distribution has not been ";
2228  warning_message+="tested with MPI.";
2229 
2230  // If we're not running the code in serial
2231  if (dist_pt->communicator_pt()->nproc()>1)
2232  {
2233  // Throw a warning
2234  OomphLibWarning(warning_message,
2235  OOMPH_CURRENT_FUNCTION,
2236  OOMPH_EXCEPTION_LOCATION);
2237  }
2238 #endif
2239 #endif
2240 
2241  // Build the vector using the distribution of the restriction matrix
2242  Restriction_self_test_vectors_storage[0].build(dist_pt);
2243 
2244  // Delete the distribution pointer
2245  delete dist_pt;
2246 
2247  // Make it a null pointer
2248  dist_pt=0;
2249 
2250  // Now assign the values to the first vector. This will be restricted down
2251  // the levels of the hierarchy then back up
2252  set_self_test_vector();
2253 
2254  // Call the restriction self-test
2255  restriction_self_test();
2256 
2257  // Set the coarest level vector in Restriction_self_test_vectors_storage
2258  // to be the last entry in Interpolation_self_test_vectors_storage. This
2259  // will be interpolated up to the finest level in interpolation_self_test()
2260  Interpolation_self_test_vectors_storage[Nlevel-1]=
2261  Restriction_self_test_vectors_storage[Nlevel-1];
2262 
2263  // Call the interpolation self-test
2264  interpolation_self_test();
2265 
2266  // Destroy all of the created restriction self-test vectors
2267  Restriction_self_test_vectors_storage.resize(0);
2268 
2269  // Destroy all of the created interpolation self-test vectors
2270  Interpolation_self_test_vectors_storage.resize(0);
2271 
2272  // Stop the clock
2273  double t_st_end=TimingHelpers::timer();
2274  double total_st_time=double(t_st_end-t_st_start);
2275  oomph_info << "\nCPU time for self-test [sec]: " << total_st_time << std::endl;
2276 
2277  // Notify user that the hierarchy of levels is complete
2278  oomph_info << "\n====================Self-Test Complete===================="
2279  << std::endl;
2280  } // End of self_test
2281 
2282  //===================================================================
2283  /// \short Sets the initial vector to be used in the restriction and
2284  /// interpolation self-tests
2285  //===================================================================
2286  template<unsigned DIM>
2288  {
2289  // Set up a pointer to the refineable mesh
2290  TreeBasedRefineableMeshBase* bulk_mesh_pt=
2291  Mg_problem_pt->mg_bulk_mesh_pt();
2292 
2293  // Find the number of elements in the bulk mesh
2294  unsigned n_el=bulk_mesh_pt->nelement();
2295 
2296  // Get the dimension of the problem (assumed to be the dimension of any
2297  // node in the mesh)
2298  unsigned n_dim=bulk_mesh_pt->node_pt(0)->ndim();
2299 
2300  // Find the number of nodes in an element
2301  unsigned nnod=bulk_mesh_pt->finite_element_pt(0)->nnode();
2302 
2303  // Loop over all elements
2304  for (unsigned e=0;e<n_el;e++)
2305  {
2306  // Get pointer to element
2307  FiniteElement* el_pt=bulk_mesh_pt->finite_element_pt(e);
2308 
2309  // Loop over nodes
2310  for (unsigned j=0;j<nnod;j++)
2311  {
2312  // Pointer to node
2313  Node* nod_pt=el_pt->node_pt(j);
2314 
2315  // Sanity check
2316  if (nod_pt->nvalue()!=1)
2317  {
2318  // Set up an output stream
2319  std::ostringstream error_stream;
2320 
2321  // Output the error message
2322  error_stream << "Sorry, not sure what to do here! I can't deal with "
2323  << nod_pt->nvalue() << " values!" << std::endl;
2324 
2325  // Throw an error to indicate that it was not possible to plot the data
2326  throw OomphLibError(error_stream.str(),
2327  OOMPH_CURRENT_FUNCTION,
2328  OOMPH_EXCEPTION_LOCATION);
2329  }
2330 
2331  // Free or pinned?
2332  int eqn_num=nod_pt->eqn_number(0);
2333 
2334  // If it's a free node
2335  if (eqn_num>=0)
2336  {
2337  // Initialise the variable coordinate_sum
2338  double coordinate_sum=0.0;
2339 
2340  // Loop over the coordinates of the node
2341  for (unsigned i=0;i<n_dim;i++)
2342  {
2343  // Increment coordinate_sum by the value of x(i)
2344  coordinate_sum+=nod_pt->x(i);
2345  }
2346 
2347  // Store the value of pi in cache
2348  double pi=MathematicalConstants::Pi;
2349 
2350  // Make the vector represent a constant function
2351  Restriction_self_test_vectors_storage[0][eqn_num]=
2352  sin(pi*(nod_pt->x(0)))*sin(pi*(nod_pt->x(1)));
2353  }
2354  } // for (unsigned j=0;j<nnod;j++)
2355  } // for (unsigned e=0;e<n_el;e++)
2356 
2357  // // Get the size of the vector
2358  // unsigned n_row=Restriction_self_test_vectors_storage[0].nrow();
2359 
2360  // // Loop over the entries of the vector
2361  // for (unsigned i=0;i<n_row;i++)
2362  // {
2363  // // Make the vector represent a constant function
2364  // Restriction_self_test_vectors_storage[0][i]=1.0;
2365  // }
2366  } // End of set_self_test_vector
2367 
2368  //===================================================================
2369  /// \short Function which implements a self-test to restrict the
2370  /// residual vector down all of the coarser grids and output
2371  /// the restricted vectors to file
2372  //===================================================================
2373  template<unsigned DIM>
2375  {
2376  // Set the start of the output file name. The full names will
2377  // set after we calculate the restricted vectors
2378  std::string outputfile="RESLT/restriction_self_test";
2379 
2380  // Loop over the levels of the hierarchy
2381  for (unsigned level=0;level<Nlevel-1;level++)
2382  {
2383  // Restrict the vector down to the next level
2384  Restriction_matrices_storage_pt[level]->
2385  multiply(Restriction_self_test_vectors_storage[level],
2386  Restriction_self_test_vectors_storage[level+1]);
2387  } // End of the for loop over the hierarchy levels
2388 
2389  // Loop over the levels of hierarchy to plot the restricted vectors
2390  for (unsigned level=0;level<Nlevel;level++)
2391  {
2392  // Set output file name
2393  std::string filename=outputfile;
2394  std::stringstream string;
2395  string << level;
2396  filename+=string.str();
2397  filename+=".dat";
2398 
2399  // Plot the RHS vector on the current level
2400  plot(level,Restriction_self_test_vectors_storage[level],filename);
2401 
2402  } // End of for-loop to output the RHS vector on each level
2403  } // End of restriction_self_test
2404 
2405  //=======================================================================
2406  /// \short Function which implements a self-test to interpolate a
2407  /// vector up all of levels of the MG hierarchy and outputs the
2408  /// restricted vectors to file
2409  //=======================================================================
2410  template<unsigned DIM>
2412  {
2413  // Set the start of the output file name. The full names will
2414  // set after we calculate the interpolated vectors
2415  std::string outputfile="RESLT/interpolation_self_test";
2416 
2417  // Loop over the levels of the hierarchy in reverse order
2418  for (unsigned level=Nlevel-1;level>0;level--)
2419  {
2420  // Interpolate the vector up a level
2421  Interpolation_matrices_storage_pt[level-1]->
2422  multiply(Interpolation_self_test_vectors_storage[level],
2423  Interpolation_self_test_vectors_storage[level-1]);
2424  } // End of the for loop over the hierarchy levels
2425 
2426  for (unsigned level=0;level<Nlevel;level++)
2427  {
2428  // Set output file name
2429  std::string filename=outputfile;
2430  std::stringstream string;
2431  string << level;
2432  filename+=string.str();
2433  filename+=".dat";
2434 
2435  // Plot the RHS vector on the current level
2436  plot(level,Interpolation_self_test_vectors_storage[level],filename);
2437 
2438  } // End of for-loop to output the RHS vector on each level
2439  } // End of interpolation_self_test
2440 
2441  //===================================================================
2442  /// \short Plots the input vector (assuming we're dealing with scalar
2443  /// nodal data, otherwise I don't know how to implement this...)
2444  //===================================================================
2445  template<unsigned DIM>
2446  void MGSolver<DIM>::plot(const unsigned& hierarchy_level,
2447  const DoubleVector& input_vector,
2448  const std::string& filename)
2449  {
2450  // Setup an output file
2451  std::ofstream some_file;
2452  some_file.open(filename.c_str());
2453 
2454  // Set up a pointer to the refineable mesh
2455  TreeBasedRefineableMeshBase* bulk_mesh_pt=
2456  Mg_hierarchy[hierarchy_level]->mg_bulk_mesh_pt();
2457 
2458  // Find the number of elements in the bulk mesh
2459  unsigned n_el=bulk_mesh_pt->nelement();
2460 
2461  // Get the dimension of the problem (assumed to be the dimension of any
2462  // node in the mesh)
2463  unsigned n_dim=bulk_mesh_pt->node_pt(0)->ndim();
2464 
2465  // Find the number of nodes in an element
2466  unsigned nnod=bulk_mesh_pt->finite_element_pt(0)->nnode();
2467 
2468  // Find the number of nodes in an element in one direction
2469  unsigned nnod_1d=bulk_mesh_pt->finite_element_pt(0)->nnode_1d();
2470 
2471  // Loop over all elements
2472  for (unsigned e=0;e<n_el;e++)
2473  {
2474  // Get pointer to element
2475  FiniteElement* el_pt=bulk_mesh_pt->finite_element_pt(e);
2476 
2477  // Input string for tecplot zone header (when plotting nnode_1d
2478  // points in each coordinate direction)
2479  some_file << el_pt->tecplot_zone_string(nnod_1d) << std::endl;
2480 
2481  // Loop over nodes
2482  for (unsigned j=0;j<nnod;j++)
2483  {
2484  // Pointer to node
2485  Node* nod_pt=el_pt->node_pt(j);
2486 
2487  // Sanity check
2488  if (nod_pt->nvalue()!=1)
2489  {
2490  std::ostringstream error_stream;
2491 
2492  error_stream << "Sorry, not sure what to do here!"
2493  << nod_pt->nvalue() << std::endl;
2494 
2495  // Output the value of dimension to screen
2496  error_stream << "The dimension is: " << n_dim << std::endl;
2497 
2498  // Output the values of x_i for all i
2499  for (unsigned i=0;i<n_dim;i++)
2500  {
2501  error_stream << nod_pt->x(i) << " ";
2502  }
2503 
2504  // End line
2505  error_stream << std::endl;
2506 
2507  // Throw an error to indicate that it was
2508  // not possible to plot the data
2509  throw OomphLibError(error_stream.str(),
2510  OOMPH_CURRENT_FUNCTION,
2511  OOMPH_EXCEPTION_LOCATION);
2512  }
2513 
2514  // Output the coordinates of the node
2515  for (unsigned i=0;i<n_dim;i++)
2516  {
2517  some_file << nod_pt->x(i) << " ";
2518  }
2519 
2520  // Free or pinned?
2521  int e=nod_pt->eqn_number(0);
2522  if (e>=0)
2523  {
2524  // Output the e-th entry if the node is free
2525  some_file << input_vector[e] << std::endl;
2526  }
2527  else
2528  {
2529  // Boundary node
2530  if (nod_pt->is_on_boundary())
2531  {
2532  // On the finest level the boundary data is pinned so set to 0
2533  if (hierarchy_level==0)
2534  {
2535  some_file << 0.0 << std::endl;
2536  }
2537  // On any other level we output this data since it corresponds
2538  // to the correction on the boundary nodes
2539  else
2540  {
2541  some_file << input_vector[e] << std::endl;
2542  }
2543  }
2544  // Hanging node
2545  else if (nod_pt->is_hanging())
2546  {
2547  // Allocate space for a double to store the weighted sum
2548  // of the surrounding master nodes of the hanging node
2549  double hang_sum=0.0;
2550 
2551  // Find the number of master nodes of the hanging
2552  // the node in the reference element
2553  HangInfo* hang_info_pt=nod_pt->hanging_pt();
2554  unsigned nmaster=hang_info_pt->nmaster();
2555 
2556  // Loop over the master nodes
2557  for (unsigned i_master=0;i_master<nmaster;i_master++)
2558  {
2559  // Set up a pointer to the master node
2560  Node* master_node_pt=hang_info_pt->master_node_pt(i_master);
2561 
2562  // Get the global equation number of this node
2563  int master_jj=master_node_pt->eqn_number(0);
2564 
2565  // Is the master node a proper d.o.f.?
2566  if (master_jj>=0)
2567  {
2568  // If the weight of the master node is non-zero
2569  if (hang_info_pt->master_weight(i_master)!=0.0)
2570  {
2571  hang_sum+=hang_info_pt->
2572  master_weight(i_master)*input_vector[master_jj];
2573  }
2574  } // if (master_jj>=0)
2575  } // for (unsigned i_master=0;i_master<nmaster;i_master++)
2576 
2577  // Output the weighted sum of the surrounding master nodes
2578  some_file << hang_sum << std::endl;
2579  }
2580  } // if (e>=0)
2581  } // for (unsigned j=0;j<nnod;j++)
2582  } // for (unsigned e=0;e<n_el;e++)
2583 
2584  // Close output file
2585  some_file.close();
2586  } // End of plot
2587 
2588  //===================================================================
2589  /// \short Linear solver
2590  //===================================================================
2591  template<unsigned DIM>
2593  {
2594  // If we're allowed to time things
2595  double t_mg_start=0.0;
2596  if (!Suppress_v_cycle_output)
2597  {
2598  // Start the clock!
2599  t_mg_start=TimingHelpers::timer();
2600  }
2601 
2602  // Current level
2603  unsigned finest_level=0;
2604 
2605  // Initialise the V-cycle counter
2606  V_cycle_counter=0;
2607 
2608  // Calculate the norm of the residual then output
2609  double normalised_residual_norm=residual_norm(finest_level);
2610  if (!Suppress_v_cycle_output)
2611  {
2612  oomph_info << "\nResidual on finest level for V-cycle: "
2613  << normalised_residual_norm << std::endl;
2614  }
2615 
2616  // Outer loop over V-cycles
2617  //-------------------------
2618  // While the tolerance is not met and the maximum number of
2619  // V-cycles has not been completed
2620  while((normalised_residual_norm>(this->Tolerance)) &&
2621  (V_cycle_counter!=Nvcycle))
2622  {
2623  if (!Suppress_v_cycle_output)
2624  {
2625  // Output the V-cycle counter
2626  oomph_info << "\nStarting V-cycle: " << V_cycle_counter << std::endl;
2627  }
2628 
2629  // Loop downwards over all levels that have coarser levels beneath them
2630  //---------------------------------------------------------------------
2631  for (unsigned i=0;i<Nlevel-1;i++)
2632  {
2633  // Initialise x_mg and residual_mg to 0.0 except for the finest level
2634  // since x_mg contains the current approximation to the solution and
2635  // residual_mg contains the RHS vector on the finest level
2636  if(i!=0)
2637  {
2638  // Loop over the entries in x_mg and residual_mg
2639  X_mg_vectors_storage[i].initialise(0.0);
2640  Residual_mg_vectors_storage[i].initialise(0.0);
2641  }
2642 
2643  // Perform a few pre-smoothing steps and return vector that contains
2644  // the residuals of the linear system at this level.
2645  pre_smooth(i);
2646 
2647  // Restrict the residual to the next coarser mesh and
2648  // assign it to the RHS vector at that level
2649  restrict_residual(i);
2650  } // Moving down the V-cycle
2651 
2652  // Reached the lowest level: Do a direct solve, using the RHS vector
2653  //------------------------------------------------------------------
2654  // obtained by restriction from above.
2655  //------------------------------------
2656  direct_solve();
2657 
2658  // Loop upwards over all levels that have finer levels above them
2659  //---------------------------------------------------------------
2660  for (unsigned i=Nlevel-1;i>0;i--)
2661  {
2662  // Interpolate solution at current level onto
2663  // next finer mesh and correct the solution x at that level
2664  interpolate_and_correct(i);
2665 
2666  // Perform a few post-smoothing steps (ignore
2667  // vector that contains the residuals of the linear system
2668  // at this level)
2669  post_smooth(i-1);
2670  }
2671 
2672  // Calculate the new residual norm then output (if allowed)
2673  normalised_residual_norm=residual_norm(finest_level);
2674 
2675  // If required, we will document the convergence history to screen or file
2676  // (if stream open)
2677  if (Doc_convergence_history)
2678  {
2679  if (!Output_file_stream.is_open())
2680  {
2681  oomph_info << V_cycle_counter << " "
2682  << normalised_residual_norm << std::endl;
2683  }
2684  else
2685  {
2686  Output_file_stream << V_cycle_counter << " "
2687  << normalised_residual_norm << std::endl;
2688  }
2689  } // if (Doc_convergence_history)
2690 
2691  // Set counter for number of cycles (for doc)
2692  V_cycle_counter++;
2693 
2694  // Print the residual on the finest level
2695  if (!Suppress_v_cycle_output)
2696  {
2697  oomph_info << "Residual on finest level of V-cycle: "
2698  << normalised_residual_norm << std::endl;
2699  }
2700  } // End of the V-cycles
2701 
2702  // Copy the solution into the result vector
2703  result=X_mg_vectors_storage[finest_level];
2704 
2705  // Need an extra line space if V-cycle output is suppressed
2706  if (!Suppress_v_cycle_output)
2707  {
2708  oomph_info << std::endl;
2709  }
2710 
2711  // If all output is to be suppressed
2712  if (!Suppress_all_output)
2713  {
2714  // Output number of V-cycles taken to solve
2715  if (normalised_residual_norm<(this->Tolerance))
2716  {
2717  // Indicate that the problem has been solved
2718  Has_been_solved=true;
2719  }
2720  } // if (!Suppress_all_output)
2721 
2722  // If the V-cycle output isn't suppressed
2723  if (!Suppress_v_cycle_output)
2724  {
2725  // Stop the clock
2726  double t_mg_end=TimingHelpers::timer();
2727  double total_G_setup_time=double(t_mg_end-t_mg_start);
2728  oomph_info << "CPU time for MG solver [sec]: "
2729  << total_G_setup_time << std::endl;
2730  }
2731  } // End of mg_solve
2732 } // End of namespace oomph
2733 
2734 #endif
MGPreconditioner(const MGPreconditioner &)
Broken copy constructor.
MGProblem * Mg_problem_pt
Pointer to the MG problem (deep copy). This is protected to provide access to the MG preconditioner...
void setup_smoothers()
Function to set up all of the smoothers once the system matrices have been set up.
void setup_mg_hierarchy()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level...
bool Doc_everything
If this is set to true we document everything. In addition to outputting the information of the setup...
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
virtual void actions_before_adapt()
Actions that are to be performed before a mesh adaptation. These might include removing any additiona...
Definition: problem.h:1004
double & tolerance()
Access to convergence tolerance.
virtual TreeBasedRefineableMeshBase * mg_bulk_mesh_pt()=0
Function to get a pointer to the mesh we will be working with. If there are flux elements present in ...
Nullstream oomph_nullstream
Single (global) instantiation of the Nullstream.
bool Suppress_v_cycle_output
Indicates whether or not the V-cycle output should be suppressed. Needs to be protected member data f...
unsigned long nnode() const
Return number of nodes in the mesh.
Definition: mesh.h:590
Vector< MGProblem * > Mg_hierarchy
Vector containing pointers to problems in hierarchy.
void disable_smoother_and_superlu_doc_time()
Suppress the output of both smoothers and SuperLU.
double Tolerance
Convergence tolerance.
std::ostream * Stream_pt
Pointer to the output stream – defaults to std::cout. This is protected member data to allow the prec...
virtual ~MGProblem()
Destructor (empty)
void restrict_residual(const unsigned &level)
Restrict residual (computed on level-th MG level) to the next coarser mesh and stick it into the coar...
cstr elem_len * i
Definition: cfortran.h:607
FiniteElement * finite_element_pt(const unsigned &e) const
Upcast (downcast?) to FiniteElement (needed to access FiniteElement member functions).
Definition: mesh.h:477
Vector< Smoother * > Pre_smoothers_storage_pt
Vector to store the pre-smoothers.
void pre_smooth(const unsigned &level)
Pre-smoother: Perform 'max_iter' smoothing steps on the linear system Ax=b with current RHS vector...
const double Pi
50 digits from maple
unsigned nrow() const
access function to the number of global rows.
virtual void actions_after_adapt()
Actions that are to be performed after a mesh adaptation.
Definition: problem.h:1007
Vector< Smoother * > Post_smoothers_storage_pt
Vector to store the post-smoothers.
void direct_solve()
Call the direct solver (SuperLU) to solve the problem exactly.
The Problem class.
Definition: problem.h:152
MGSolver(MGProblem *mg_problem_pt)
Constructor: Set up default values for number of V-cycles and pre- and post-smoothing steps...
bool Has_been_setup
Boolean variable to indicate whether or not the solver has been setup.
void setup_interpolation_matrices_unstructured()
Setup the interpolation matrix on each level (used for unstructured meshes)
A general Finite Element class.
Definition: elements.h:1271
Smoother *(* PostSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class Smoother to be used as the po...
virtual bool is_on_boundary() const
Test whether the Node lies on a boundary. The "bulk" Node cannot lie on a boundary, so return false. This will be overloaded by BoundaryNodes.
Definition: nodes.h:1293
void setup_mg_structures()
Function to set up the hierachy of levels. Creates a vector of pointers to each MG level...
PostSmootherFactoryFctPt Post_smoother_factory_function_pt
Function to create post-smoothers.
void multiply(const CRDoubleMatrix *matrix, const DoubleVector &x, DoubleVector &soln)
Function to perform a matrix-vector multiplication on a oomph-lib matrix and vector using Trilinos fu...
void set_post_smoother_factory_function(PostSmootherFactoryFctPt post_smoother_fn)
Access function to set the post-smoother creation function.
void locate_zeta(const Vector< double > &zeta, GeomObject *&sub_geom_object_pt, Vector< double > &s, const bool &use_coordinate_as_initial_guess=false)
Find the sub geometric object and local coordinate therein that corresponds to the intrinsic coordina...
void clean_up_memory()
Clean up memory.
Nodes are derived from Data, but, in addition, have a definite (Eulerian) position in a space of a gi...
Definition: nodes.h:852
void plot(const unsigned &hierarchy_level, const DoubleVector &input_vector, const std::string &filename)
Given a level in the hierarchy, an input vector and a filename this function will document the given ...
Mesh *& mesh_pt()
Return a pointer to the global mesh.
Definition: problem.h:1264
HangInfo *const & hanging_pt() const
Return pointer to hanging node data (this refers to the geometric hanging node status) (const version...
Definition: nodes.h:1148
OomphInfo oomph_info
void setup()
Function to set up a preconditioner for the linear system.
long & eqn_number(const unsigned &i)
Return the equation number of the i-th stored variable.
Definition: nodes.h:365
unsigned Npre_smooth
Number of pre-smoothing steps.
unsigned Nvcycle
Maximum number of V-cycles (this is set as a protected variable so.
e
Definition: cfortran.h:575
double * value()
Access to C-style value array.
Definition: matrices.h:1062
unsigned nvalue() const
Return number of values stored in data object (incl pinned ones).
Definition: nodes.h:448
void setup_interpolation_matrices()
Setup the interpolation matrix on each level.
virtual unsigned nnode_1d() const
Return the number of nodes along one edge of the element Default is to return zero — must be overload...
Definition: elements.h:2139
PreSmootherFactoryFctPt Pre_smoother_factory_function_pt
Function to create pre-smoothers.
Vector< DoubleVector > X_mg_vectors_storage
Vector to store the solution vectors (X_mg)
unsigned nmaster() const
Return the number of master nodes.
Definition: nodes.h:733
virtual MGProblem * make_new_problem()=0
This function needs to be implemented in the derived problem: Returns a pointer to a new object of th...
void set_pre_smoother_factory_function(PreSmootherFactoryFctPt pre_smoother_fn)
Access function to set the pre-smoother creation function.
unsigned long nelement() const
Return number of elements in the mesh.
Definition: mesh.h:587
unsigned iterations() const
Number of iterations.
double & x(const unsigned &i)
Return the i-th nodal coordinate.
Definition: nodes.h:995
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void enable_doc_everything()
Enable the output from anything that could have been suppressed.
virtual void preconditioner_solve(const DoubleVector &rhs, DoubleVector &z)
Function applies MG to the vector r for a full solve.
~MGSolver()
Delete any dynamically allocated data.
bool is_hanging() const
Test whether the node is geometrically hanging.
Definition: nodes.h:1207
void full_setup()
Do a full setup (assumes everything will be setup around the MGProblem pointer given in the construct...
unsigned self_test()
Self-test: Return 0 for OK.
std::ostream *& stream_pt()
Access function for the stream pointer.
bool Doc_time
Boolean flag that indicates whether the time taken.
Definition: linear_solver.h:84
Base class for tree-based refineable meshes.
void mg_solve(DoubleVector &result)
Do the actual solve – this is called through the pure virtual solve function in the LinearSolver base...
double residual_norm(const unsigned &level)
Return norm of residual r=b-Ax and the residual vector itself on the level-th level.
static char t char * s
Definition: cfortran.h:572
Smoother *(* PreSmootherFactoryFctPt)()
typedef for a function that returns a pointer to an object of the class Smoother to be used as the pr...
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:992
void set_restriction_matrices_as_interpolation_transposes()
Builds a CRDoubleMatrix on each level that is used to restrict the residual between levels...
Vector< CRDoubleMatrix * > Mg_matrices_storage_pt
Vector to store the system matrices.
unsigned long ndof() const
Return the number of dofs.
Definition: problem.h:1574
void level_up_local_coord_of_node(const int &son_type, Vector< double > &s)
Given the son_type of an element and a local node number j in that element with nnode_1d nodes per co...
Vector< DoubleVector > Interpolation_self_test_vectors_storage
Vector to store the result of interpolation on each level (only required if the user wishes to docume...
void interpolate_and_correct(const unsigned &level)
Interpolate solution at current level onto next finer mesh and correct the solution x at that level...
void set_self_test_vector()
Makes a vector which will be used in the self-test. Is currently set to make the entries of the vecto...
void modify_restriction_matrices()
Normalise the rows of the restriction matrices to avoid amplifications when projecting to the coarser...
void enable_v_cycle_output()
Enable the output of the V-cycle timings and other output.
double const & master_weight(const unsigned &i) const
Return weight for dofs on i-th master node.
Definition: nodes.h:753
Vector< DoubleVector > Residual_mg_vectors_storage
Vector to store the residual vectors.
void interpolation_matrix_set(const unsigned &level, double *value, int *col_index, int *row_st, unsigned &ncol, unsigned &nnz)
Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be...
void post_smooth(const unsigned &level)
Post-smoother: Perform max_iter smoothing steps on the linear system Ax=b with current RHS vector...
Node *& node_pt(const unsigned &n)
Return a pointer to the local node n.
Definition: elements.h:2097
double timer()
returns the time in seconds after some point in past
bool Suppress_all_output
If this is set to true then all output from the solver is suppressed. This is protected member data s...
void clean_up_memory()
Clean up anything that needs to be cleaned up.
Class that contains data for hanging nodes.
Definition: nodes.h:684
An interface to allow scalar MG to be used as a Preconditioner.
Node *& node_pt(const unsigned long &n)
Return pointer to global node n.
Definition: mesh.h:456
unsigned Nlevel
The number of levels in the multigrid heirachy.
void restriction_self_test()
Make a self-test to make sure that the interpolation matrices are doing the same thing to restrict th...
void interpolation_matrix_set(const unsigned &level, Vector< double > &value, Vector< int > &col_index, Vector< int > &row_st, unsigned &ncol, unsigned &nrow)
Builds a CRDoubleMatrix that is used to interpolate the residual between levels. The transpose can be...
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
static const int OMEGA
Default value for an unassigned neighbour.
Definition: tree.h:241
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
Vector< DoubleVector > Restriction_self_test_vectors_storage
Vector to store the result of restriction on each level (only required if the user wishes to document...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
MGPreconditioner(MGProblem *mg_problem_pt)
Constructor.
Vector< DoubleVector > Rhs_mg_vectors_storage
Vector to store the RHS vectors (Rhs_mg). This is protected to allow the multigrid preconditioner to ...
unsigned nnode() const
Return the number of nodes.
Definition: elements.h:2134
void solve(Problem *const &problem_pt, DoubleVector &result)
Virtual function in the base class that needs to be implemented later but for now just leave it empty...
unsigned & npost_smooth()
Return the number of post-smoothing iterations (lvalue)
~MGPreconditioner()
Destructor (empty)
virtual std::string tecplot_zone_string(const unsigned &nplot) const
Return string for tecplot zone header (when plotting nplot points in each "coordinate direction") ...
Definition: elements.h:2949
unsigned long assign_eqn_numbers(const bool &assign_local_eqn_numbers=true)
Assign all equation numbers for problem: Deals with global data (= data that isn't attached to any el...
Definition: problem.cc:1965
Node *const & master_node_pt(const unsigned &i) const
Return a pointer to the i-th master node.
Definition: nodes.h:736
unsigned ndim() const
Return (Eulerian) spatial dimension of the node.
Definition: nodes.h:992
Vector< CRDoubleMatrix * > Interpolation_matrices_storage_pt
Vector to store the interpolation matrices.
void self_test()
Makes a vector, restricts it down the levels of the hierarchy and documents it at each level...
OomphCommunicator * communicator_pt() const
const access to the communicator pointer
void operator=(const MGPreconditioner &)
Broken assignment operator.
A vector in the mathematical sense, initially developed for linear algebra type applications. If MPI then this vector can be distributed - its distribution is described by the LinearAlgebraDistribution object at Distribution_pt. Data is stored in a C-style pointer vector (double*)
Definition: double_vector.h:61
void clear()
clears the distribution
unsigned V_cycle_counter
Pointer to counter for V-cycles.
MGProblem class; subclass of Problem.
int * row_start()
Access to C-style row_start array.
Definition: matrices.h:1050
Vector< CRDoubleMatrix * > Restriction_matrices_storage_pt
Vector to store the restriction matrices.
void interpolation_self_test()
Make a self-test to make sure that the interpolation matrices are doing the same thing to interpolate...
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
void enable_output()
Enable the output from anything that could have been suppressed.
void disable_output()
Suppress anything that can be suppressed, i.e. any timings. Things like mesh adaptation can not howev...
virtual void shape(const Vector< double > &s, Shape &psi) const =0
Calculate the geometric shape functions at local coordinate s. This function must be overloaded for e...
A general mesh class.
Definition: mesh.h:74
unsigned & npre_smooth()
Return the number of pre-smoothing iterations (lvalue)
MGProblem()
Constructor. Initialise pointers to coarser and finer levels.
void setup_transfer_matrices()
Setup the transfer matrices on each level.
Base class for all linear iterative solvers. This merely defines standard interfaces for linear itera...
void position(Vector< double > &pos) const
Compute Vector of nodal positions either directly or via hanging node representation.
Definition: nodes.cc:2399
unsigned Npost_smooth
Number of post-smoothing steps.
void disable_v_cycle_output()
Disable all output from mg_solve apart from the number of V-cycles used to solve the problem...