block_preconditioner.h
Go to the documentation of this file.
1 //LIC// ====================================================================
2 //LIC// This file forms part of oomph-lib, the object-oriented,
3 //LIC// multi-physics finite-element library, available
4 //LIC// at http://www.oomph-lib.org.
5 //LIC//
6 //LIC// Version 1.0; svn revision $LastChangedRevision: 1225 $
7 //LIC//
8 //LIC// $LastChangedDate: 2016-10-10 13:02:38 +0100 (Mon, 10 Oct 2016) $
9 //LIC//
10 //LIC// Copyright (C) 2006-2016 Matthias Heil and Andrew Hazel
11 //LIC//
12 //LIC// This library is free software; you can redistribute it and/or
13 //LIC// modify it under the terms of the GNU Lesser General Public
14 //LIC// License as published by the Free Software Foundation; either
15 //LIC// version 2.1 of the License, or (at your option) any later version.
16 //LIC//
17 //LIC// This library is distributed in the hope that it will be useful,
18 //LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
19 //LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 //LIC// Lesser General Public License for more details.
21 //LIC//
22 //LIC// You should have received a copy of the GNU Lesser General Public
23 //LIC// License along with this library; if not, write to the Free Software
24 //LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
25 //LIC// 02110-1301 USA.
26 //LIC//
27 //LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
28 //LIC//
29 //LIC//====================================================================
30 //Include guards
31 #ifndef OOMPH_BLOCK_PRECONDITION_HEADER
32 #define OOMPH_BLOCK_PRECONDITION_HEADER
33 
34 
35 // Config header generated by autoconfig
36 #ifdef HAVE_CONFIG_H
37 #include <oomph-lib-config.h>
38 #endif
39 
40 // c++ include
41 #include <math.h>
42 #include <typeinfo>
43 
44 // oomph-lib includes
45 #include "matrices.h"
46 #include "mesh.h"
47 #include "vector_matrix.h"
48 
49 // #include "problem.h"
50 #include "preconditioner.h"
51 #include "SuperLU_preconditioner.h"
52 #include "matrix_vector_product.h"
53 
54 namespace oomph
55 {
56 
57 //=============================================================================
58 /// \short Data structure to store information about a certain "block" or
59 /// sub-matrix from the overall matrix in the block preconditioning framework.
60 ///
61 /// Example of use: Let's assume we want to form a concatenated matrix
62 /// from the blocks of a Jacobian matrix that contains the following blocks:
63 ///
64 /// [J_00, J_01, J_02
65 /// J_10, J_11, J_12
66 /// J_20, J_21, J_22]
67 ///
68 /// so that the new matrix has the entries
69 ///
70 /// [ J_01, J_00
71 /// J_21, 0 ]
72 ///
73 /// where "0" indicates zero matrix of the required size.
74 ///
75 /// To do this we create a 2x2 (the block size of the new concatenated
76 /// matrix) VectorMatrix of BlockSelectors and then declare for each
77 /// entry the block indices in the original matrix and if the entry
78 /// is to be included (and copied from the corresponding entry in
79 /// the Jacobian (final boolean argument true) or if the block is
80 /// to be omitted and replaced by an appropriately sized zero matrix.
81 /// For the example above this would be done as follows:
82 ///
83 /// VectorMatrix<BlockSelector> required_block(2,2);
84 /// required_block[0][0].select_block(0,1,true);
85 /// required_block[0][1].select_block(0,0,true);
86 /// required_block[1][0].select_block(2,1,true);
87 /// required_block[1][1].select_block(2,0,false);
88 ///
89 /// and the concatenated matrix would then be built as
90 ///
91 /// CRDoubleMatrix concatenated_block1
92 /// = get_concatenated_block(required_block);
93 ///
94 /// Note that it is necessary to identify the row and column indices of any
95 /// omitted blocks (here block J_20 in the original matrix) to enable
96 /// the correct setup of the sparse matrix storage.
97 ///
98 /// The initial assignment of the boolean may be over-written with the
99 /// do_not_want_block() member function; this can again be reversed
100 /// with the want_block() counterpart. So if we call
101 ///
102 /// required_block[0][0].do_not_want_block();
103 ///
104 /// and the build a new conctatenated matrix with
105 ///
106 /// CRDoubleMatrix concatenated_block2
107 /// = get_concatenated_block(required_block);
108 ///
109 /// the resulting matrix would the anti-diagonal matrix
110 ///
111 /// [ 0 , J_00
112 /// J_21 , 0 ]
113 ///
114 /// Finally it is possible to specify a replacement block
115 /// by specifying a pointer to an appropriately sized matrix
116 /// that is to be used instead of the block in the Jacobian
117 /// matrix, so if replacement_block_pt points to a matrix, R, say,
118 /// of the same size as J_01, then
119 ///
120 /// selected_block[0][0].select_block(0,1,true,replacement_block_pt);
121 ///
122 /// then the resulting concatenated matrix would contain
123 ///
124 /// [ R , J_00
125 /// J_21 , 0 ]
126 ///
127 //=============================================================================
129 {
130  public:
131 
132  /// \short Default constructor,
133  /// initialise block index i, j to 0 and bool to false.
135  {
136  // Needs to be set to zero because if the build function leaves the
137  // Replacement_block_pt alone if replacement_block_pt = 0 (the default argument).
139  this->build(0,0,false);
140  }
141 
142  /// \short Constructor, takes the row and column indices
143  /// and a boolean indicating if the block is required or not. The optional
144  /// parameter replacement_block_pt is set to null. If the block is not required a block
145  /// of the correct dimensions full of 0s is used.
146  BlockSelector(const unsigned& row_index,
147  const unsigned& column_index,
148  const bool& wanted,
150  {
151 #ifdef PARANOID
152  if((wanted == false) && (replacement_block_pt !=0))
153  {
154  std::ostringstream err_msg;
155  err_msg << "Trying to construct a BlockSelector object with:\n"
156  << "replacement_block_pt != 0 and wanted == false"
157  << "If you require the block, please set wanted == true.\n";
158  throw OomphLibError(err_msg.str(),
159  OOMPH_CURRENT_FUNCTION,
160  OOMPH_EXCEPTION_LOCATION);
161  }
162 #endif
163 
164  // Needs to be set to zero because if the build function leaves the
165  // Replacement_block_pt alone if replacement_block_pt = 0 (the default argument).
166  // Thus if it is not set here, it would not be initialised to null.
168 
169  this->build(row_index,column_index,wanted,replacement_block_pt);
170  }
171 
172  /// \short Default destructor.
173  virtual ~BlockSelector()
174  {
175 #ifdef PARANOID
176  if(Replacement_block_pt != 0)
177  {
178  std::ostringstream warning_msg;
179  warning_msg << "Warning: BlockSelector destructor is called but...\n"
180  << "replacement_block_pt() is not null.\n"
181  << "Please remember to null this via the function\n"
182  << "BlockSelector::null_replacement_block_pt()\n"
183  << "Row_index: " << Row_index << "\n"
184  << "Column_index: " << Column_index << std::endl;
185 
186  OomphLibWarning(warning_msg.str(),
187  OOMPH_CURRENT_FUNCTION,
188  OOMPH_EXCEPTION_LOCATION);
189  }
190 #endif
191  }
192 
193  /// \short Select a block.
194  void select_block(const unsigned& row_index,
195  const unsigned& column_index,
196  const bool& wanted,
198  {
199 #ifdef PARANOID
200  if((wanted == false) && (replacement_block_pt !=0))
201  {
202  std::ostringstream err_msg;
203  err_msg << "Trying to select block with:\n"
204  << "replacement_block_pt != 0 and wanted == false"
205  << "If you require the block, please set wanted == true.\n"
206  << "row_index: " << row_index << "\n"
207  << "column_index: " << column_index << "\n";
208  throw OomphLibError(err_msg.str(),
209  OOMPH_CURRENT_FUNCTION,
210  OOMPH_EXCEPTION_LOCATION);
211  }
212 #endif
213 
214  this->build(row_index,column_index,wanted,replacement_block_pt);
215  }
216 
217 
218  /// \short Indicate that we require the block (set Wanted to true).
219  void want_block()
220  {
221  Wanted = true;
222  }
223 
224  /// \short Indicate that we do not want the block (set Wanted to false).
226  {
227 #ifdef PARANOID
228  if(Replacement_block_pt !=0)
229  {
230  std::ostringstream warning_msg;
231  warning_msg << "Trying to set Wanted = false, but replacement_block_pt is not null.\n"
232  << "Please call null_replacement_block_pt()\n"
233  << "(remember to free memory if necessary)\n"
234  << "Row_index: " << Row_index << "\n"
235  << "Column_index: " << Column_index << "\n";
236  OomphLibWarning(warning_msg.str(),
237  OOMPH_CURRENT_FUNCTION,
238  OOMPH_EXCEPTION_LOCATION);
239  }
240 #endif
241 
243 
244  Wanted = false;
245  }
246 
247  /// \short Set Replacement_block_pt to null.
249  {
251  }
252 
253  /// \short set Replacement_block_pt.
255  {
256 #ifdef PARANOID
257  if(Wanted == false)
258  {
259  std::ostringstream err_msg;
260  err_msg << "Trying to set replacement_block_pt, but Wanted == false.\n"
261  << "Please call want_block()\n"
262  << "Row_index: " << Row_index << "\n"
263  << "Column_index: " << Column_index << "\n";
264  throw OomphLibError(err_msg.str(),
265  OOMPH_CURRENT_FUNCTION,
266  OOMPH_EXCEPTION_LOCATION);
267  }
268 #endif
269 
271  }
272 
273  /// \short Returns Replacement_block_pt
275  {
276  return Replacement_block_pt;
277  }
278 
279  /// \short Set the row index.
280  void set_row_index(const unsigned& row_index)
281  {
283  }
284 
285  /// \short returns the row index.
286  const unsigned& row_index() const
287  {
288  return Row_index;
289  }
290 
291  /// \short Set the column index.
292  void set_column_index(const unsigned& column_index)
293  {
295  }
296 
297  /// \short returns the column index.
298  const unsigned& column_index() const
299  {
300  return Column_index;
301  }
302 
303  /// \short returns whether the block is wanted or not.
304  const bool& wanted() const
305  {
306  return Wanted;
307  }
308 
309 
310  /// \short Output function, outputs the Row_index, Column_index, Wanted and
311  /// the address of the Replacement_block_pt.
312  friend std::ostream& operator<<(std::ostream& o_stream,
313  const BlockSelector& block_selector)
314  {
315  o_stream << "Row_index = " << block_selector.row_index() << "\n"
316  << "Column_index = " << block_selector.column_index() << "\n"
317  << "Wanted = " << block_selector.wanted() << "\n"
318  << "Replacement_block_pt = " << block_selector.replacement_block_pt();
319 
320  return o_stream;
321  }
322 
323  private:
324 
325  /// Build function, sets the Row_index, Column_index and Wanted variables.
326  /// the Replacement_block_pt is only set if it is not null. Otherwise it is left alone.
327  void build(const unsigned& row_index,
328  const unsigned& column_index,
329  const bool& wanted,
331  {
334  Wanted = wanted;
335 
336  // Only set the replacement_block_pt if it is wanted. Otherwise we leave it alone.
337  // All constructors should set Replacement_block_pt to 0.
338  if(replacement_block_pt != 0)
339  {
340 #ifdef PARANOID
341  if(Wanted == false)
342  {
343  std::ostringstream err_msg;
344  err_msg << "Trying to set replacement_block_pt, but Wanted == false.\n"
345  << "Please either not set the replacement_block_pt or call the function\n"
346  << "do_not_want_block()\n"
347  << "Row_index: " << Row_index << "\n"
348  << "Column_index: " << Column_index << "\n";
349  throw OomphLibError(err_msg.str(),
350  OOMPH_CURRENT_FUNCTION,
351  OOMPH_EXCEPTION_LOCATION);
352  }
353 #endif
354 
356  }
357  }
358 
359  /// Row index of the block.
360  unsigned Row_index;
361 
362  /// Column index of the block.
363  unsigned Column_index;
364 
365  /// Bool to indicate if we require this block.
366  bool Wanted;
367 
368  /// Pointer to the block.
370 
371 };
372 
373 
374  //============================================================================
375  /// Block Preconditioner base class. The block structure of the
376  /// overall problem is determined from the \c Mesh's constituent
377  /// elements. Each constituent element must be block-preconditionable - i.e
378  /// must implement the \c GeneralisedElements functions \c ndof_types() and
379  /// get_dof_numbers_for_unknowns(...). A \c Problem can have several
380  /// \c Meshes, but each \c Mesh must contain elements with the same DOF types.
381  /// The association between global degrees of freedom and their unique local
382  /// dof numbers is therefore based on information provided by the elements.
383  /// We refer to the local dof numbers provided by the elements as the
384  /// elemental dof numbers.
385  ///
386  /// By default each type of DOF is assumed to be unique type of block,
387  /// but DOF types can be grouped together in a single block when
388  /// block_setup(...) is called.
389  ///
390  /// This class can function in one of two ways. Either it acts as a
391  /// stand-alone block preconditioner which computes and stores
392  /// the association between global degrees of freedom and their unique global
393  /// block numbers itself. Alternatively, the block preconditioner can act as
394  /// a subsidiary block preconditioner within a (larger) master block
395  /// preconditioner (pointed to by Master_block_preconditioner_pt).
396  /// The master block preconditioner
397  /// must have an equal or greater number of block types. Examples
398  /// are the FSI preconditioner which is the 3x3 "master block preconditioner"
399  /// for the Navier-Stokes preconditioners which deals with the
400  /// 2x2 fluid-blocks within the overall structure. In this case, \b only
401  /// the master block preconditioner computes and stores the master
402  /// lookup schemes. All block preconditioners compute and store their own
403  /// optimised lookup schemes.
404  ///
405  /// In cases where a \c Problem contains elements of different element types
406  /// (e.g. fluid and solid elements in a fluid-structure interaction problem),
407  /// access to the elements of the same type must be provided via pointers to
408  /// (possibly auxiliary) \c Meshes that only contain elements of a single
409  /// type. The block preconditioner will then create global block
410  /// numbers by concatenating the block types. Consider, e.g. a fluid-structure
411  /// interaction problem in which the first \c Mesh contains (fluid)
412  /// elements whose degrees of freedom have been subdivided into
413  /// types "0" (the velocity, say) and "1" (the pressure say), while
414  /// the second \c Mesh contains (solid) elements whose degrees of freedom
415  /// are the nodal displacements, classified as its type "0".
416  /// The combined block preconditioner then has three "block types":
417  /// "0": Fluid velocity, "1": Fluid pressure, "2": Solid nodal positions.
418  /// NOTE: currently this preconditioner uses the same communicator as the
419  /// underlying problem. We may need to change this in the future.
420  //============================================================================
421  template<typename MATRIX>
423  {
424  public:
425 
426  /// \short Constructor
428  : Ndof_types_in_mesh(0)
429  {
430  // Initially set the master block preconditioner pointer to zero
431  // indicating that this is stand-alone preconditioner (i.e. the upper most
432  // level preconditioner) that will set up its own block lookup schemes etc.
434 
435  // The distribution of the concatenation of the internal block
436  // distributions.
437  // I.e. LinearAlgebraDistributionHelpers::concatenate
438  // (distributions of internal blocks).
439  //
440  // The concatenation of the internal block distributions is stored in two
441  // places depending on if this is the upper-most master block preconditioner
442  // or not.
443  //
444  // If this is a master block preconditioner (Master_block_preconditioner_pt
445  // is null), then it is stored in the variable
446  // Internal_preconditioner_matrix_distribution_pt (below).
447  // For subsidiary block preconditioners, this remains null.
448  //
449  // Because BlockPreconditioners are DistributedLinearAlgebraObjects, they
450  // have a distribution. For the upper-most master block preconditioner,
451  // this is the distribution of the underlying Jacobian.
452  //
453  // For all subsidiary block preconditioners, this remains null. The
454  // concatenation of the distribution of the internal blocks are stored
455  // as the distribution of the BlockPreconditioner.
456  //
457  // This seems inconsistent and I cannot figure out why this is done.
459 
460  // The concatenation of the external block distributions.
462 
463  // Initialise number of rows in this block preconditioner.
464  // This information is maintained if used as subsidiary or
465  // stand-alone block preconditioner (in the latter case it
466  // obviously stores the number of rows within the subsidiary
467  // preconditioner.
468  Nrow=0;
469 
470  // Initialise number of different block types in this preconditioner.
471  // This information is maintained if used as subsidiary or
472  // stand-alone block preconditioner (in the latter case it
473  // obviously stores the number of rows within the subsidiary
474  // preconditioner.
476 
477  // Initialise number of different dof types in this preconditioner.
478  // This information is maintained if used as subsidiary or
479  // stand-alone block preconditioner (in the latter case it
480  // obviously stores the number of rows within the subsidiary
481  // preconditioner.
483 
484  // There are no blocks to start off with.
485  Block_distribution_pt.resize(0);
486 
487  // The distributions of the underlying internal blocks.
489 
490  // The distribution of the dof-level blocks, these are used during the
491  // concatenation process to create the distribution of the blocks.
492  Dof_block_distribution_pt.resize(0);
493 
494  // Clear both the Block_to_dof_map_coarse and Block_to_dof_map_fine
495  // vectors.
496  Block_to_dof_map_coarse.resize(0);
497  Block_to_dof_map_fine.resize(0);
498 
499  // Default the debug flag to false.
500  Recursive_debug_flag = false;
501 
502  // Default the debug flag to false.
503  Debug_flag = false;
504  } // EOFunc constructor
505 
506 
507 
508  /// Destructor
510  {
512  } // EOFunc destructor
513 
514  /// Broken copy constructor
516  {
517  BrokenCopy::broken_copy("BlockPreconditioner");
518  }
519 
520  /// Broken assignment operator
522  {
523  BrokenCopy::broken_assign("BlockPreconditioner");
524  }
525 
526  /// \short Access function to matrix_pt. If this is the master then cast
527  /// the matrix pointer to MATRIX*, error check and return. Otherwise ask
528  /// the master for its matrix pointer.
529  MATRIX* matrix_pt() const
530  {
532  {
533  return master_block_preconditioner_pt()->matrix_pt();
534  }
535  else
536  {
537  MATRIX* m_pt = dynamic_cast<MATRIX*>(Preconditioner::matrix_pt());
538 #ifdef PARANOID
539  if(m_pt == 0)
540  {
541  std::ostringstream error_msg;
542  error_msg << "Matrix is not correct type.";
543  throw OomphLibError(error_msg.str(),
544  OOMPH_CURRENT_FUNCTION,
545  OOMPH_EXCEPTION_LOCATION);
546  }
547 #endif
548  return m_pt;
549  }
550  } // EOFunc matrix_pt()
551 
552 
553  /// \short Toggles on the recursive debug flag. The change goes
554  /// up the block preconditioning hierarchy.
556  {
557  Recursive_debug_flag = true;
560  ->turn_on_recursive_debug_flag();
561  }
562 
563  /// \short Toggles off the recursive debug flag. The change goes
564  /// up the block preconditioning hierarchy.
566  {
567  Recursive_debug_flag = false;
570  ->turn_off_recursive_debug_flag();
571  }
572 
573  /// \short Toggles on the debug flag.
575  {
576  Debug_flag = true;
577  }
578 
579  /// \short Toggles off the debug flag.
581  {
582  Debug_flag = false;
583  }
584 
585  /// \short Function to turn this preconditioner into a
586  /// subsidiary preconditioner that operates within a bigger
587  /// "master block preconditioner (e.g. a Navier-Stokes 2x2 block
588  /// preconditioner dealing with the fluid sub-blocks within a
589  /// 3x3 FSI preconditioner. Once this is done the master block
590  /// preconditioner deals with the block setup etc.
591  /// The vector doftype_in_master_preconditioner_coarse must specify the
592  /// dof number in the master preconditioner that corresponds to a dof number
593  /// in this preconditioner.
594  /// \b 1. The length of the vector is used to determine the number of
595  /// blocks in this preconditioner therefore it must be correctly sized.
596  /// \b 2. block_setup(...) should be called in the master preconditioner
597  /// before this method is called.
598  /// \b 3. block_setup(...) should be called in the corresponding subsidiary
599  /// preconditioner after this method is called.
600  ///
601  /// This calls the other turn_into_subsidiary_block_preconditioner
602  /// function with the identity mapping for doftype_coarsen_map_coarse
603  /// vector.
605  (BlockPreconditioner<MATRIX>* master_block_prec_pt,
606  const Vector<unsigned>& doftype_in_master_preconditioner_coarse);
607 
608  /// \short Function to turn this preconditioner into a
609  /// subsidiary preconditioner that operates within a bigger
610  /// "master block preconditioner (e.g. a Navier-Stokes 2x2 block
611  /// preconditioner dealing with the fluid sub-blocks within a
612  /// 3x3 FSI preconditioner. Once this is done the master block
613  /// preconditioner deals with the block setup etc.
614  /// The vector doftype_in_master_preconditioner_coarse must specify the
615  /// dof number in the master preconditioner that corresponds to a dof number
616  /// in this preconditioner.
617  /// \b 1. The length of the vector is used to determine the number of
618  /// blocks in this preconditioner therefore it must be correctly sized.
619  /// \b 2. block_setup(...) should be called in the master preconditioner
620  /// before this method is called.
621  /// \b 3. block_setup(...) should be called in the corresponding subsidiary
622  /// preconditioner after this method is called.
623  ///
624  /// The doftype_coarsen_map_coarse is a mapping of the dof numbers in the
625  /// master preconditioner to the dof numbers REQUIRED by THIS preconditioner.
626  /// This allows for coarsening of the dof types if the master preconditioner
627  /// has a more fine grain dof type splitting.
628  ///
629  /// For example, the Lagrangian preconditioner (in 3D with one constraint)
630  /// has doftypes:
631  /// 0 1 2 3 4 5 6 7
632  /// ub vb wb uc vc wc p Lc
633  ///
634  /// We wish to use an existing Navier-Stokes preconditioner, for example,
635  /// LSC, to solve the sub-system associated with the dof numbers
636  /// 0, 1, 2, 3, 4, 5, 6. But the existing LSC preconditioner only works
637  /// with four dof types (3 velocity dof types and 1 pressure dof types).
638  /// We need to coarsen the number of dof types in the master preconditioner.
639  ///
640  /// If the LSC preconditioner requires the dof ordering: u, v, w, p. Then
641  /// the doftype_coarsen_map_coarse will be:
642  /// [0 3] -> u velocity dof type
643  /// [1 4] -> v velocity dof type
644  /// [2 5] -> w velocity dof type
645  /// [6] -> pressure dof type.
647  (BlockPreconditioner<MATRIX>* master_block_prec_pt,
648  const Vector<unsigned>& doftype_in_master_preconditioner_coarse,
649  const Vector<Vector<unsigned> > & doftype_coarsen_map_coarse);
650 
651 
652 
653  /// \short Determine the size of the matrix blocks and setup the
654  /// lookup schemes relating the global degrees of freedom with
655  /// their "blocks" and their indices (row/column numbers) in those
656  /// blocks.
657  /// The distributions of the preconditioner and the internal blocks are
658  /// automatically specified (and assumed to be uniform) at this
659  /// stage.
660  /// This method should be used if the identity dof-to-block mapping is okay,
661  /// i.e.
662  /// dof number 0 corresponds to block number 0
663  /// dof number 1 corresponds to block number 1
664  /// dof number 2 corresponds to block number 2
665  /// etc...
666  virtual void block_setup();
667 
668  /// \short Determine the size of the matrix blocks and setup the
669  /// lookup schemes relating the global degrees of freedom with
670  /// their "blocks" and their indices (row/column numbers) in those
671  /// blocks.
672  /// The distributions of the preconditioner and the blocks are
673  /// automatically specified (and assumed to be uniform) at this
674  /// stage.
675  /// This method should be used if anything other than the identity
676  /// dof-to-block mapping is required. The argument vector dof_to_block_map
677  /// should be of length ndof. The indices represents the dof types whilst the
678  /// value represents the block types. In general we want:
679  ///
680  /// dof_to_block_map[dof_number] = block_number.
681  void block_setup(const Vector<unsigned>& dof_to_block_map);
682 
683  /// \short Put block (i,j) into output_matrix. This block accounts for any
684  /// coarsening of dof types and any replaced dof-level blocks above this
685  /// preconditioner.
686  void get_block(const unsigned& i, const unsigned& j,
687  MATRIX& output_matrix,
688  const bool& ignore_replacement_block = false) const
689  {
690 #ifdef PARANOID
691  // Check the range of i and j, they should not be greater than nblock_types
692  unsigned n_block_types = this->nblock_types();
693  if((i > n_block_types) || (j > n_block_types))
694  {
695  std::ostringstream err_msg;
696  err_msg << "Requested block("<< i <<","<< j <<")"<<"\n"
697  << "but nblock_types() is " << n_block_types << ".\n"
698  << "Maybe your argument to block_setup(...) is incorrect?"
699  << std::endl;
700  throw OomphLibError(err_msg.str(),
701  OOMPH_CURRENT_FUNCTION,
702  OOMPH_EXCEPTION_LOCATION);
703  }
704 #endif
705 
706  // The logic is this:
707  //
708  // Block_to_dof_map_coarse tells us which dof types belongs in each block.
709  // This is relative to this preconditioner, and describes the external
710  // block and dof type mappings (what the preconditioner writer
711  // expects/sees).
712  //
713  // As such, the dof types in Block_to_dof_map_coarse describes the
714  // dof-level blocks needed to be concatenated to produce block(i,j). These
715  // dofs may have been coarsened.
716  //
717  // Now, the dof blocks to concatenate comes from:
718  // If the dof block exists in Replacement_dof_block_pt, then we make a
719  // deep copy of it. Otherwise, if this is the upper-most master block
720  // preconditioner then we get it from the original matrix via function
721  // internal_get_block(...) otherwise, if this is a subsidiary
722  // block preconditioner, we go one level up the hierarchy and repeat the
723  // process.
724  //
725  //
726  // A small note about indirections which caused me some headache.
727  // Thought I would mention it here in case the below code needs to be
728  // re-visited.
729  //
730  // A small subtlety with indirections:
731  //
732  // The underlying ordering of the dof-level blocks is STILL AND ALWAYS the
733  // `natural' ordering determined by first the elements then the order of
734  // the meshes.
735  //
736  // But during the process of block_setup(...), the external (perceived)
737  // block ordering may have changed. So some indirection has to take place,
738  // this mapping is maintained in Block_to_dof_map_coarse.
739  //
740  // For example, take the Lagrangian preconditioner, which expects the
741  // natural dof type ordering:
742  // 0 1 2 3 4 5
743  // u v p uc vc L
744  //
745  // If the mapping given to block_setup(...) is:
746  // dof_to_block_map = [0, 1, 4, 2, 3, 5]
747  // saying that: dof type 0 goes to block 0
748  // dof type 1 goes to block 1
749  // dof type 2 goes to block 4
750  // dof type 3 goes to block 2
751  // dof type 4 goes to block 3
752  // dof type 5 goes to block 5
753  //
754  // The function block_setup(...) will populate the vector
755  // Block_to_dof_map_coarse with [[0][1][3][4][2][5]],
756  // which says that get_block(0,0) will get the u block
757  // get_block(1,1) will get the v block
758  // get_block(2,2) will get the uc block
759  // get_block(3,3) will get the vc block
760  // get_block(4,4) will get the p block
761  // get_block(5,5) will get the L block
762  //
763  // i.e. the ordering of the block types is a permutation of the dof types,
764  // so that we now have the following block ordering:
765  // 0 1 2 3 4 5 <- block ordering
766  // u v uc vc p L
767  //
768  // Now, the writer expects to work with the block ordering. I.e. when we
769  // replace a dof-level block, say the pressure block, we call
770  // set_replacement_dof_block(4,4,Matrix);
771  // Similarly, when we want the pressure block, we call
772  // get_block(4,4);
773  //
774  // Now, the below code uses the indirection maintained in
775  // Block_to_dof_map_coarse. I.e. When we call get_block(4,4), we use the
776  // mapping Block_to_dof_map_coarse[4] = 2, we get the block (2,2)
777  // (from Replacement_dof_block_pt or internal_get_block), since the
778  // underlying dof_to_block mapping is still the identity, i.e. it has not
779  // changed from:
780  // 0 1 2 3 4 5
781  // u v p uc vc L
782  //
783  // So, block (4,4) (pressure block) maps to the block (2,2).
784 
785  // How many external dof types are in this block?
786  const unsigned ndofblock_in_row = Block_to_dof_map_coarse[i].size();
787  const unsigned ndofblock_in_col = Block_to_dof_map_coarse[j].size();
788 
789  // If there is only one dof block in this block then there is no need to
790  // concatenate.
791  if(ndofblock_in_row == 1 && ndofblock_in_col == 1)
792  {
793  // Get the indirection for the dof we want.
794  const unsigned wanted_dof_row = Block_to_dof_map_coarse[i][0];
795  const unsigned wanted_dof_col = Block_to_dof_map_coarse[j][0];
796 
797  // If the block has NOT been replaced or if we want to ignore the
798  // replacement, then we call get_dof_level_block(...) which will get the
799  // dof-level blocks up the preconditioning hierarchy, dragging down
800  // any replacement dof blocks of parent preconditioners if required.
801  if((Replacement_dof_block_pt.get(wanted_dof_row,wanted_dof_col) == 0) ||
802  ignore_replacement_block)
803  {
804  get_dof_level_block(wanted_dof_row, wanted_dof_col,
805  output_matrix, ignore_replacement_block);
806  }
807  else
808  // Replacement_dof_block_pt.get(wanted_dof_row,wanted_dof_col) is not
809  // null, this means that the block has been replaced. We simply make
810  // a deep copy of it.
811  {
813  Replacement_dof_block_pt.get(wanted_dof_row, wanted_dof_col),
814  output_matrix);
815  }
816  }
817  else
818  // This block contains more than one dof-level block. So we need to
819  // concatenate the (external) dof-level blocks.
820  {
821  // The CRDoubleMatrixHelpers::concatenate_without_communication(...)
822  // takes a DenseMatrix of pointers to CRDoubleMatrices to concatenate.
823  DenseMatrix<CRDoubleMatrix*> tmp_blocks_pt(ndofblock_in_row,
824  ndofblock_in_col,0);
825 
826  // Vector of Vectors of unsigns to indicate whether we have created
827  // CRDoubleMatrices with new or not... so we can delete it later.
828  // 0 - no new CRDoubleMatrix is created.
829  // 1 - a new CRDoubleMatrix is created.
830  // If we ever use C++11, remove this and use smart pointers.
831  Vector<Vector<unsigned> >new_block(ndofblock_in_row,
832  Vector<unsigned>(ndofblock_in_col,0));
833 
834  // Loop through the number of dof block rows and then the number of dof
835  // block columns, either get the pointer from Replacement_dof_block_pt
836  // or from get_dof_level_block(...).
837  for (unsigned row_dofblock = 0; row_dofblock < ndofblock_in_row;
838  row_dofblock++)
839  {
840  // Indirection for the row, as discuss in the large chunk of text
841  // previously.
842  const unsigned wanted_dof_row
843  = Block_to_dof_map_coarse[i][row_dofblock];
844 
845  for (unsigned col_dofblock = 0; col_dofblock < ndofblock_in_col;
846  col_dofblock++)
847  {
848  // Indirection for the column as discussed in the large chunk of text
849  // above.
850  const unsigned wanted_dof_col
851  = Block_to_dof_map_coarse[j][col_dofblock];
852 
853  // Get the pointer from Replacement_dof_block_pt.
854  tmp_blocks_pt(row_dofblock,col_dofblock)
855  = Replacement_dof_block_pt.get(wanted_dof_row,wanted_dof_col);
856 
857  // If the pointer from Replacement_dof_block_pt is null, or if
858  // we have to ignore replacement blocks, build a new CRDoubleMatrix
859  // via get_dof_level_block.
860  if((tmp_blocks_pt(row_dofblock,col_dofblock) == 0) ||
861  ignore_replacement_block)
862  {
863  // We have to create a new CRDoubleMatrix, so put in 1 to indicate
864  // that we have to delete it later.
865  new_block[row_dofblock][col_dofblock] = 1;
866 
867  // Create the new CRDoubleMatrix. Note that we do not use the
868  // indirection, since the indirection is only used one way.
869  tmp_blocks_pt(row_dofblock,col_dofblock) = new CRDoubleMatrix;
870 
871  // Get the dof-level block, as discussed above.
872  get_dof_level_block(wanted_dof_row,
873  wanted_dof_col,
874  *tmp_blocks_pt(row_dofblock,col_dofblock),
875  ignore_replacement_block);
876  }
877  }
878  }
879 
880  // Concatenation of matrices require the distribution of the individual
881  // sub-matrices (for both row and column). This is because concatenation
882  // is implemented without communication in such a way that the rows
883  // and column values are both permuted, the permutation is defined by
884  // the individual distributions of the sub blocks.
885  // Without a vector of distributions describing the distributions of
886  // of the rows, we do not know how to permute the rows. For the columns,
887  // although CRDoubleMatrices does not have a column distribution, the
888  // concatenated matrix must have it's columns permuted, so we mirror
889  // the diagonal and get the corresponding row distribution.
890  //
891  // Confused? - Example: Say we have a matrix with dof blocking
892  //
893  // | a | b | c | d | e |
894  // --|---|---|---|---|---|
895  // a | | | | | |
896  // --|---|---|---|---|---|
897  // b | | | | | |
898  // --|---|---|---|---|---|
899  // c | | | | | |
900  // --|---|---|---|---|---|
901  // d | | | | | |
902  // --|---|---|---|---|---|
903  // e | | | | | |
904  // --|---|---|---|---|---|
905  //
906  // We wish to concatenate the blocks
907  //
908  // | d | e |
909  // --|---|---|
910  // a | | |
911  // --|---|---|
912  // b | | |
913  // --|---|---|
914  //
915  // Then clearly the row distributions required are the distributions for
916  // the dof blocks a and b. But block(a,d) does not have a column
917  // distribution since it is a CRDoubleMatrix! - We use the distribution
918  // mirrored by the diagonal, so the column distributions required to
919  // concatenate these blocks is the same as the distributions of the rows
920  // for dof block d and e.
921 
922  // First we do the row distribution.
923 
924  // Storage for the row distribution pointers.
925  Vector<LinearAlgebraDistribution*> tmp_row_dist_pt(ndofblock_in_row,0);
926 
927  // Loop through the number of dof blocks in the row. For the upper-most
928  // master block preconditioner, the external dof distributions is the
929  // same as the internal BLOCK distributions. Recall what we said above
930  // about the underlying blocks maintaining it's natural ordering.
931  //
932  // If this is a subsidiary block preconditioner, then the distributions
933  // for the dof blocks are stored in Dof_block_distribution_pt. The reason
934  // why this is different for subsidiary block preconditioners is because
935  // subsidiary block preconditioners would have it's dof types coarsened.
936  // Then a single dof distribution in a subsidiary block preconditioner
937  // could be a concatenation of many dof distributions of the master
938  // dof distributions.
939  for (unsigned row_dof_i = 0; row_dof_i < ndofblock_in_row; row_dof_i++)
940  {
941  const unsigned mapped_dof_i = Block_to_dof_map_coarse[i][row_dof_i];
943  {
944  tmp_row_dist_pt[row_dof_i]
945  = Internal_block_distribution_pt[mapped_dof_i];
946  }
947  else
948  {
949  tmp_row_dist_pt[row_dof_i]
950  = Dof_block_distribution_pt[mapped_dof_i];
951  }
952  }
953 
954  // Storage for the column distribution pointers.
955  Vector<LinearAlgebraDistribution*> tmp_col_dist_pt(ndofblock_in_col,0);
956 
957  // We do the same thing as before.
958  for (unsigned col_dof_i = 0; col_dof_i < ndofblock_in_col; col_dof_i++)
959  {
960  const unsigned mapped_dof_j = Block_to_dof_map_coarse[j][col_dof_i];
962  {
963  tmp_col_dist_pt[col_dof_i]
964  = Internal_block_distribution_pt[mapped_dof_j];
965  }
966  else
967  {
968  tmp_col_dist_pt[col_dof_i]
969  = Dof_block_distribution_pt[mapped_dof_j];
970  }
971  }
972 
973  // Perform the concatenation.
975  tmp_col_dist_pt,
976  tmp_blocks_pt,
977  output_matrix);
978 
979  // Delete any new CRDoubleMatrices we have created.
980  for (unsigned row_i = 0; row_i < ndofblock_in_row; row_i++)
981  {
982  for (unsigned col_i = 0; col_i < ndofblock_in_col; col_i++)
983  {
984  if(new_block[row_i][col_i])
985  {
986  delete tmp_blocks_pt(row_i,col_i);
987  }
988  }
989  }
990  }// else need to concatenate
991  } // EOFunc get_block(...)
992 
993 
994  /// \short Return block (i,j). If the optional argument
995  /// ignore_replacement_block is true, then any blocks in
996  /// Replacement_dof_block_pt will be ignored throughout the preconditioning
997  /// hierarchy.
998  MATRIX get_block(const unsigned& i, const unsigned& j,
999  const bool& ignore_replacement_block = false) const
1000  {
1001  MATRIX output_matrix;
1002  get_block(i, j, output_matrix, ignore_replacement_block);
1003  return output_matrix;
1004  } // EOFunc get_block(...)
1005 
1006  /// \short Set the matrix_pt in the upper-most master preconditioner.
1007  void set_master_matrix_pt(MATRIX* in_matrix_pt)
1008  {
1010  {
1012  ->set_master_matrix_pt(in_matrix_pt);
1013  }
1014  else
1015  {
1016  this->set_matrix_pt(in_matrix_pt);
1017  }
1018  }
1019 
1020  /// \short Get a block from a different matrix using the blocking scheme
1021  /// that has already been set up.
1022  void get_block_other_matrix(const unsigned& i, const unsigned& j,
1023  MATRIX* in_matrix_pt,
1024  MATRIX& output_matrix)
1025  {
1026  MATRIX* backup_matrix_pt = matrix_pt();
1027  set_master_matrix_pt(in_matrix_pt);
1028  get_block(i, j, output_matrix, true);
1029  set_master_matrix_pt(backup_matrix_pt);
1030  } // EOFunc get_block_other_matrix(...)
1031 
1032  /// \short Get all the block matrices required by the block
1033  /// preconditioner. Takes a pointer to a matrix of bools that indicate
1034  /// if a specified sub-block is required for the preconditioning
1035  /// operation. Computes the required block matrices, and stores pointers
1036  /// to them in the matrix block_matrix_pt. If an entry in block_matrix_pt
1037  /// is equal to NULL on return, that sub-block has not been requested and
1038  /// is therefore not available.
1039  ///
1040  /// WARNING: the matrix pointers are created using new so you must delete
1041  /// them all manually!
1042  ///
1043  /// WARNING 2: the matrix pointers in block_matrix_pt MUST be null
1044  /// because Richard in all his wisdom decided to call delete on any
1045  /// non-null pointers. Presumably to avoid fixing his memory leaks
1046  /// properly...
1047  void get_blocks(DenseMatrix<bool>& required_blocks,
1048  DenseMatrix<MATRIX*>& block_matrix_pt) const;
1049 
1050  /// \short Gets dof-level block (i,j).
1051  /// If Replacement_dof_block_pt(i,j) is not null, then the replacement
1052  /// block is returned via a deep copy.
1053  ///
1054  /// Otherwise if this is the uppermost block preconditioner then it calls
1055  /// internal_get_block(i,j), else if it is a subsidiary
1056  /// block preconditioner, it will call it's master block preconditioners'
1057  /// get_dof_level_block function.
1058  void get_dof_level_block(const unsigned& i, const unsigned& j,
1059  MATRIX& output_block,
1060  const bool& ignore_replacement_block = false) const;
1061 
1062 
1063  /// \short Returns a concatenation of the block matrices specified by the
1064  /// argument selected_block. The VectorMatrix selected_block must be
1065  /// correctly sized as it is used to determine the number of sub block
1066  /// matrices to concatenate.
1067  ///
1068  /// For each entry in the VectorMatrix, the following variables must
1069  /// correctly set:
1070  /// BlockSelector::Row_index - Refers to the row index of the block.
1071  /// BlockSelector::Column_index - Refers to the column index of the block.
1072  /// BlockSelector::Wanted - Do we want the block?
1073  /// BlockSelector::Replacement_block_pt - If not null, this block will be used instead of
1074  /// get_block(Row_index,Column_index).
1075  ///
1076  /// For example, assume that we have a matrix of the following blocking:
1077  /// 0 1 2 3 4
1078  /// | a | b | c | d | e |
1079  /// --|---|---|---|---|---|
1080  /// 0 a | | | | | |
1081  /// --|---|---|---|---|---|
1082  /// 1 b | | | | | |
1083  /// --|---|---|---|---|---|
1084  /// 2 c | | | | | |
1085  /// --|---|---|---|---|---|
1086  /// 3 d | | | | | |
1087  /// --|---|---|---|---|---|
1088  /// 4 e | | | | | |
1089  /// --|---|---|---|---|---|
1090  ///
1091  /// If we want a block matrix corresponding to the concatenation of the
1092  /// blocks [(a,d), (a,e)
1093  /// , (b,e)* ]
1094  /// where top left and top right blocks comes from the function
1095  /// get_block(...), the bottom left entry is empty, and the bottom right is
1096  /// a modified block.
1097  ///
1098  /// Then we create a VectorMatrix of size 2 by 2
1099  /// VectorMatrix<BlockSelector> selected_block(2,2);
1100  ///
1101  /// In the [0][0] entry:
1102  /// row index is 0,
1103  /// column index is 3,
1104  /// do we want this block? - yes (true).
1105  /// selected_block[0][0].select_block(0,3,true);
1106  ///
1107  /// In the [0][1] entry:
1108  /// row index is 0,
1109  /// column index is 4,
1110  /// do we want this block? - yes (true).
1111  /// selected_block[0][0].select_block(0,4,true);
1112  ///
1113  /// In the [1][0] entry:
1114  /// row index is 1,
1115  /// column index is 3,
1116  /// do we want this block? - no (false).
1117  /// selected_block[0][0].select_block(1,3,false);
1118  ///
1119  /// In the [1][1] entry:
1120  /// row index is 1,
1121  /// column index is 4,
1122  /// do we want this block? - yes (true).
1123  /// selected_block[0][0].select_block(1,4,true,block_pt);
1124  ///
1125  /// where block_pt is a pointer to the modified block.
1126  ///
1127  /// Then we can call:
1128  ///
1129  /// CRDoubleMatrix my_block = get_concatenated_block(selected_block);
1130  ///
1131  /// NOTE: It is important to set the row and column indices even if you do
1132  /// not want the block. This is because, if we allow the row and column
1133  /// indices to be "not set", then we can have a whole empty block row without
1134  /// any indices. But concatenation of blocks without communication requires
1135  /// both the row and column distributions, so we know how to permute the
1136  /// rows and columns. So in short, we require that the column and row
1137  /// indices to always be set for every entry in the
1138  /// VectorMatrix<BlockSelector> object for convenience and consistency
1139  /// checks.
1141  &selected_block)
1142  {
1143 #ifdef PARANOID
1144 
1145  unsigned const para_selected_block_nrow = selected_block.nrow();
1146  unsigned const para_selected_block_ncol = selected_block.ncol();
1147  unsigned const para_nblocks = this->nblock_types();
1148 
1149  // Check that selected_block size is not 0.
1150  if(para_selected_block_nrow == 0)
1151  {
1152  std::ostringstream error_msg;
1153  error_msg << "selected_block has nrow = 0.\n";
1154  throw OomphLibError(error_msg.str(),
1155  OOMPH_CURRENT_FUNCTION,
1156  OOMPH_EXCEPTION_LOCATION);
1157  }
1158 
1159  // Check that the number of blocks is not outside of the range
1160  // nblock_types(). Since this function builds matrices for block
1161  // preconditioning, it does not make sense for us to concatenate more
1162  // blocks than nblock_types().
1163  if((para_selected_block_nrow > para_nblocks) ||
1164  (para_selected_block_ncol > para_nblocks))
1165  {
1166  std::ostringstream error_msg;
1167  error_msg << "Trying to concatenate a "
1168  << para_selected_block_nrow << " by "
1169  << para_selected_block_ncol
1170  << " block matrix,\n"
1171  << "but there are only "
1172  << para_nblocks << " block types.\n";
1173  throw OomphLibError(error_msg.str(),
1174  OOMPH_CURRENT_FUNCTION,
1175  OOMPH_EXCEPTION_LOCATION);
1176  }
1177 
1178  // Check that selected_block make sense, i.e. the row indices of each row
1179  // are the same, and the column indices of each column are the same.
1180 
1181  // First check if the row indices are consistent.
1182  // For each row, loop through the columns, comparing the row index against
1183  // the first column.
1184  for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
1185  {
1186  const unsigned col_0_row_index = selected_block[row_i][0].row_index();
1187 
1188  for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
1189  {
1190  const unsigned para_b_i = selected_block[row_i][col_i].row_index();
1191  const unsigned para_b_j = selected_block[row_i][col_i].column_index();
1192 
1193  if(col_0_row_index != para_b_i)
1194  {
1195  std::ostringstream err_msg;
1196  err_msg << "Block index for block(" << row_i << "," << col_i <<") "
1197  << "contains block indices (" << para_b_i << ","
1198  << para_b_j << ").\n"
1199  << "But the row index for the first column is "
1200  << col_0_row_index <<", they must be the same!\n";
1201  throw OomphLibError(err_msg.str(),
1202  OOMPH_CURRENT_FUNCTION,
1203  OOMPH_EXCEPTION_LOCATION);
1204  }
1205  }
1206  }
1207 
1208  // Do the same for the column indices, consistency check.
1209  // For each column, loop through the rows, comparing the column index
1210  // against the first row.
1211  for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
1212  {
1213  const unsigned row_0_col_index = selected_block[0][col_i].column_index();
1214 
1215  for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
1216  {
1217  const unsigned para_b_i = selected_block[row_i][col_i].row_index();
1218  const unsigned para_b_j = selected_block[row_i][col_i].column_index();
1219 
1220  if(row_0_col_index != para_b_j)
1221  {
1222  std::ostringstream err_msg;
1223  err_msg << "Block index for block(" << row_i << "," << col_i <<") "
1224  << "contains block indices (" << para_b_i << ","
1225  << para_b_j << ").\n"
1226  << "But the col index for the first row is "
1227  << row_0_col_index <<", they must be the same!\n";
1228  throw OomphLibError(err_msg.str(),
1229  OOMPH_CURRENT_FUNCTION,
1230  OOMPH_EXCEPTION_LOCATION);
1231  }
1232  }
1233  }
1234 
1235  // Check to see if the values in selected_block is within the range
1236  // nblock_types()
1237  //
1238  // Since we know that the column and row indices are consistent (by the
1239  // two paranoia checks above), we only need to check the column indices
1240  // in the first row, and the row indices in the first column.
1241 
1242  // Check that the row indices in the first column are within the range
1243  // nblock_types()
1244  for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
1245  {
1246  const unsigned para_b_i = selected_block[row_i][0].row_index();
1247  const unsigned para_b_j = selected_block[row_i][0].column_index();
1248  if(para_b_i > para_nblocks)
1249  {
1250  std::ostringstream err_msg;
1251  err_msg << "Block index for block(" << row_i << ",0) "
1252  << "contains block indices (" << para_b_i << ","
1253  << para_b_j << ").\n"
1254  << "But there are only " << para_nblocks
1255  << " nblock_types().\n";
1256  throw OomphLibError(err_msg.str(),
1257  OOMPH_CURRENT_FUNCTION,
1258  OOMPH_EXCEPTION_LOCATION);
1259  }
1260  }
1261 
1262  // Check that the col indices in the first row are within the range
1263  // nblock_types()
1264  for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
1265  {
1266  const unsigned para_b_i = selected_block[0][col_i].row_index();
1267  const unsigned para_b_j = selected_block[0][col_i].column_index();
1268  if(para_b_j > para_nblocks)
1269  {
1270  std::ostringstream err_msg;
1271  err_msg << "Block index for block(0," << col_i << ") "
1272  << "contains block indices (" << para_b_i << ","
1273  << para_b_j << ").\n"
1274  << "But there are only " << para_nblocks
1275  << " nblock_types().\n";
1276  throw OomphLibError(err_msg.str(),
1277  OOMPH_CURRENT_FUNCTION,
1278  OOMPH_EXCEPTION_LOCATION);
1279  }
1280  }
1281 
1282  // Stricter test - can be removed is required in the future. For the first
1283  // column, check that the row indices does not repeat.
1284  std::set<unsigned> row_index_set;
1285  for (unsigned row_i = 0; row_i < para_selected_block_nrow; row_i++)
1286  {
1287  std::pair<std::set<unsigned>::iterator,bool> row_index_set_ret;
1288 
1289  const unsigned row_i_index = selected_block[row_i][0].row_index();
1290 
1291  row_index_set_ret
1292  = row_index_set.insert(row_i_index);
1293 
1294  if(!row_index_set_ret.second)
1295  {
1296  std::ostringstream err_msg;
1297  err_msg << "The row " << row_i_index << " is already inserted.\n";
1298  throw OomphLibError(err_msg.str(),
1299  OOMPH_CURRENT_FUNCTION,
1300  OOMPH_EXCEPTION_LOCATION);
1301  }
1302  }
1303 
1304  // Stricter test - can be removed is required in the future. For the first
1305  // row, check that the column indices does not repeat.
1306  std::set<unsigned> col_index_set;
1307  for (unsigned col_i = 0; col_i < para_selected_block_ncol; col_i++)
1308  {
1309  std::pair<std::set<unsigned>::iterator,bool> col_index_set_ret;
1310 
1311  const unsigned col_i_index = selected_block[0][col_i].column_index();
1312 
1313  col_index_set_ret
1314  = col_index_set.insert(col_i_index);
1315 
1316  if(!col_index_set_ret.second)
1317  {
1318  std::ostringstream err_msg;
1319  err_msg << "The col " << col_i_index << " is already inserted.\n";
1320  throw OomphLibError(err_msg.str(),
1321  OOMPH_CURRENT_FUNCTION,
1322  OOMPH_EXCEPTION_LOCATION);
1323  }
1324  }
1325 
1326  // Loop through all the block_pt and check:
1327  // 1) The non-null pointers point to built matrices.
1328  // 2) The distribution matches those defined by selected_block within
1329  // Block_distribution_pt.
1330  for (unsigned block_i = 0; block_i < para_selected_block_nrow; block_i++)
1331  {
1332  for (unsigned block_j = 0; block_j < para_selected_block_ncol; block_j++)
1333  {
1334  const CRDoubleMatrix* tmp_block_pt
1335  = selected_block[block_i][block_j].replacement_block_pt();
1336 
1337  if(tmp_block_pt != 0)
1338  {
1339  if(!tmp_block_pt->built())
1340  {
1341  std::ostringstream err_msg;
1342  err_msg << "The matrix pointed to by block("
1343  << block_i << "," << block_j << ") is not built.\n";
1344  throw OomphLibError(err_msg.str(),
1345  OOMPH_CURRENT_FUNCTION,
1346  OOMPH_EXCEPTION_LOCATION);
1347  }
1348 
1349  const LinearAlgebraDistribution* const tmp_block_dist_pt
1350  = tmp_block_pt->distribution_pt();
1351 
1352  const unsigned row_selected_block
1353  = selected_block[block_i][block_j].row_index();
1354 
1355  const LinearAlgebraDistribution* const another_tmp_block_dist_pt
1356  = block_distribution_pt(row_selected_block);
1357 
1358  if(*tmp_block_dist_pt != *another_tmp_block_dist_pt)
1359  {
1360  std::ostringstream err_msg;
1361  err_msg << "block_distribution_pt " << row_selected_block << "\n"
1362  << "does not match the distribution from the block_pt() "
1363  << " in selected_block["
1364  << block_i << "][" << block_j <<"].\n";
1365  throw OomphLibError(err_msg.str(),
1366  OOMPH_CURRENT_FUNCTION,
1367  OOMPH_EXCEPTION_LOCATION);
1368  }
1369  }
1370  }
1371  }
1372 
1373  // Attempt a similar check for the column index. This is not as rigorous
1374  // since a CRDoubleMatrix does not have a distribution for the columns.
1375  // However, we can check if the ncol of the matrices in block_pt matches
1376  // those in the block_distribution_pt corresponding to the columns.
1377  // (I hope this makes sense... both the row and columns are permuted in
1378  // CRDoubleMatrixHelpers::concatenate_without_communication(...))
1379  //
1380  // The test for the row distributions checks if the nrow_local is correct.
1381  // We do not have the equivalent for columns.
1382  for (unsigned block_i = 0; block_i < para_selected_block_nrow; block_i++)
1383  {
1384  for (unsigned block_j = 0; block_j < para_selected_block_ncol; block_j++)
1385  {
1386  // Cache the block_pt
1387  const CRDoubleMatrix* tmp_block_pt
1388  = selected_block[block_i][block_j].replacement_block_pt();
1389 
1390  if(tmp_block_pt != 0)
1391  {
1392  const unsigned tmp_block_ncol = tmp_block_pt->ncol();
1393 
1394  const unsigned selected_block_col
1395  = selected_block[block_i][block_j].column_index();
1396 
1397  // YES, nrow, this is not incorrect.
1398  const unsigned another_tmp_block_ncol
1399  = block_distribution_pt(selected_block_col)->nrow();
1400 
1401  if(tmp_block_ncol != another_tmp_block_ncol)
1402  {
1403  std::ostringstream err_msg;
1404  err_msg << "block_pt in selected_block["
1405  << block_i << "][" << block_j << "] "
1406  << "has ncol = " << tmp_block_ncol << ".\n"
1407  << "But the corresponding block_distribution_pt("
1408  << selected_block_col << ") has nrow = "
1409  << another_tmp_block_ncol
1410  << ").\n";
1411  throw OomphLibError(err_msg.str(),
1412  OOMPH_CURRENT_FUNCTION,
1413  OOMPH_EXCEPTION_LOCATION);
1414  }
1415  }
1416  }
1417  }
1418 #endif
1419 
1420  // The return matrix.
1421  MATRIX output_matrix;
1422 
1423  // How many sub matrices are there in the row and column?
1424  const unsigned nblock_row = selected_block.nrow();
1425  const unsigned nblock_col = selected_block.ncol();
1426 
1427  // Get the row and col distributions, this is required for concatenation
1428  // without communication.
1429  Vector<LinearAlgebraDistribution*> row_dist_pt(nblock_row,0);
1430  Vector<LinearAlgebraDistribution*> col_dist_pt(nblock_col,0);
1431 
1432  // For the row distributions, use the first column of selected_block
1433  // Also, store the index of the block rows.
1434  Vector<unsigned> block_row_index(nblock_row,0);
1435  for (unsigned row_i = 0; row_i < nblock_row; row_i++)
1436  {
1437  const unsigned selected_row_index
1438  = selected_block[row_i][0].row_index();
1439 
1440  row_dist_pt[row_i] = Block_distribution_pt[selected_row_index];
1441  block_row_index[row_i] = selected_row_index;
1442  }
1443 
1444  // For the col distributions, use the first row of selected_block
1445  Vector<unsigned> block_col_index(nblock_col,0);
1446  for (unsigned col_i = 0; col_i < nblock_col; col_i++)
1447  {
1448  const unsigned selected_col_index
1449  = selected_block[0][col_i].column_index();
1450 
1451  col_dist_pt[col_i] = Block_distribution_pt[selected_col_index];
1452  block_col_index[col_i] = selected_col_index;
1453  }
1454 
1455  // Now build the output matrix. The output_matrix needs a distribution,
1456  // this distribution is a concatenation of the block rows. But because
1457  // concatenation of distributions requires communication, we try to
1458  // minimise this process by creating it once, then store a key to the
1459  // concatenated distribution. First check to see if the block row indices
1460  // is already a key in Auxiliary_block_distribution_pt, if it is in there,
1461  // we use the distribution it corresponds to. Otherwise, we create the
1462  // distribution and store it for possible further use.
1463  std::map<Vector<unsigned>,
1464  LinearAlgebraDistribution* >::const_iterator iter;
1465 
1466  iter = Auxiliary_block_distribution_pt.find(block_row_index);
1467 
1468  if(iter != Auxiliary_block_distribution_pt.end())
1469  {
1470  output_matrix.build(iter->second);
1471  }
1472  else
1473  {
1475  LinearAlgebraDistributionHelpers::concatenate(row_dist_pt,*tmp_dist_pt);
1476  insert_auxiliary_block_distribution(block_row_index,tmp_dist_pt);
1477  output_matrix.build(tmp_dist_pt);
1478  }
1479 
1480  // Do the same for the column dist, since we might need it for the RHS
1481  // vector..
1482  iter = Auxiliary_block_distribution_pt.find(block_col_index);
1483  if(iter == Auxiliary_block_distribution_pt.end())
1484  {
1486  LinearAlgebraDistributionHelpers::concatenate(col_dist_pt,*tmp_dist_pt);
1487  insert_auxiliary_block_distribution(block_col_index,tmp_dist_pt);
1488  }
1489 
1490  // Storage for the pointers to CRDoubleMatrices to concatenate.
1491  DenseMatrix<CRDoubleMatrix*> block_pt(nblock_row,nblock_col,0);
1492 
1493  // Vector of Vectors of unsigns to indicate whether we have created
1494  // CRDoubleMatrices with new or not... so we can delete it later.
1495  // 0 - no new CRDoubleMatrix is created.
1496  // 1 - a new CRDoubleMatrix is created.
1497  // If we ever use C++11, remove this and use smart pointers.
1498  Vector<Vector<unsigned> > new_block(nblock_row,
1499  Vector<unsigned>(nblock_col,0));
1500 
1501  // Get blocks if wanted.
1502  for (unsigned block_i = 0; block_i < nblock_row; block_i++)
1503  {
1504  for (unsigned block_j = 0; block_j < nblock_col; block_j++)
1505  {
1506  const bool block_wanted = selected_block[block_i][block_j].wanted();
1507 
1508  if(block_wanted)
1509  {
1510  CRDoubleMatrix* tmp_block_pt
1511  = selected_block[block_i][block_j].replacement_block_pt();
1512 
1513  if(tmp_block_pt == 0)
1514  {
1515  new_block[block_i][block_j] = 1;
1516 
1517  block_pt(block_i,block_j) = new CRDoubleMatrix;
1518 
1519  // temp variables for readability purposes.
1520  const unsigned tmp_block_i = block_row_index[block_i];
1521  const unsigned tmp_block_j = block_col_index[block_j];
1522 
1523  // Get the block.
1524  this->get_block(tmp_block_i,tmp_block_j,
1525  *block_pt(block_i,block_j));
1526  }
1527  else
1528  {
1529  block_pt(block_i,block_j) = tmp_block_pt;
1530  }
1531  }
1532  }
1533  }
1534 
1535  // Perform the concatenation.
1537  col_dist_pt,
1538  block_pt,
1539  output_matrix);
1540 
1541  // Delete any new CRDoubleMatrices we created.
1542  for (unsigned block_i = 0; block_i < nblock_row; block_i++)
1543  {
1544  for (unsigned block_j = 0; block_j < nblock_col; block_j++)
1545  {
1546  if(new_block[block_i][block_j])
1547  {
1548  delete block_pt(block_i,block_j);
1549  }
1550  }
1551  }
1552 
1553  return output_matrix;
1554  } // EOFunc get_concatenated_block(...)
1555 
1556  /// \short Takes the naturally ordered vector and extracts the blocks
1557  /// indicated by the block number (the values) in the Vector
1558  /// block_vec_number all at once, then concatenates them without
1559  /// communication. Here, the values in block_vec_number is the block number
1560  /// in the current preconditioner.
1561  /// This is a non-const function because distributions may be created
1562  /// and stored in Auxiliary_block_distribution_pt for future use.
1563  void get_concatenated_block_vector(const Vector<unsigned>& block_vec_number,
1564  const DoubleVector& v,
1565  DoubleVector& b);
1566 
1567  /// \short Takes concatenated block ordered vector, b, and copies its
1568  /// entries to the appropriate entries in the naturally ordered vector, v.
1569  /// Here the values in block_vec_number indicates which blocks the vector
1570  /// b is a concatenation of. The block number are those in the current
1571  /// preconditioner. If the preconditioner is a subsidiary block
1572  /// preconditioner the other entries in v that are not associated with it
1573  /// are left alone.
1575  const Vector<unsigned>& block_vec_number,
1576  const DoubleVector& b,
1577  DoubleVector& v) const;
1578 
1579  /// \short Takes the naturally ordered vector and rearranges it into a
1580  /// vector of sub vectors corresponding to the blocks, so s[b][i] contains
1581  /// the i-th entry in the vector associated with block b.
1582  /// Note: If the preconditioner is a subsidiary preconditioner then only the
1583  /// sub-vectors associated with the blocks of the subsidiary preconditioner
1584  /// will be included. Hence the length of v is master_nrow() whereas the
1585  /// total length of the s vectors is the sum of the lengths of the
1586  /// individual block vectors defined in block_vec_number.
1587  void get_block_vectors(const Vector<unsigned> & block_vec_number,
1588  const DoubleVector& v,
1589  Vector<DoubleVector >& s) const;
1590 
1591  /// \short Takes the naturally ordered vector and rearranges it into a
1592  /// vector of sub vectors corresponding to the blocks, so s[b][i] contains
1593  /// the i-th entry in the vector associated with block b.
1594  /// Note: If the preconditioner is a subsidiary preconditioner then only the
1595  /// sub-vectors associated with the blocks of the subsidiary preconditioner
1596  /// will be included. Hence the length of v is master_nrow() whereas the
1597  /// total length of the s vectors is Nrow.
1598  /// This is simply a wrapper around the other get_block_vectors(...) function
1599  /// where the block_vec_number Vector is the identity, i.e.
1600  /// block_vec_number is [0, 1, ..., nblock_types - 1].
1601  void get_block_vectors(const DoubleVector& v,
1602  Vector<DoubleVector >& s) const;
1603 
1604  /// \short Takes the vector of block vectors, s, and copies its entries into
1605  /// the naturally ordered vector, v. If this is a subsidiary block
1606  /// preconditioner only those entries in v that are associated with its
1607  /// blocks are affected. The block_vec_number indicates which block the
1608  /// vectors in s came from. The block number corresponds to the block
1609  /// numbers in this preconditioner.
1610  void return_block_vectors(const Vector<unsigned>& block_vec_number,
1611  const Vector<DoubleVector >& s,
1612  DoubleVector& v) const;
1613 
1614  /// \short Takes the vector of block vectors, s, and copies its entries into
1615  /// the naturally ordered vector, v. If this is a subsidiary block
1616  /// preconditioner only those entries in v that are associated with its
1617  /// blocks are affected. The block_vec_number indicates which block the
1618  /// vectors in s came from. The block number corresponds to the block
1619  /// numbers in this preconditioner.
1620  /// This is simply a wrapper around the other return_block_vectors(...)
1621  /// function where the block_vec_number Vector is the identity, i.e.
1622  /// block_vec_number is [0, 1, ..., nblock_types - 1].
1624  DoubleVector& v) const;
1625 
1626  /// \short Takes the naturally ordered vector, v and returns the n-th
1627  /// block vector, b. Here n is the block number in the current
1628  /// preconditioner.
1629  void get_block_vector(const unsigned& n, const DoubleVector& v,
1630  DoubleVector& b) const;
1631 
1632  /// \short Takes the n-th block ordered vector, b, and copies its entries
1633  /// to the appropriate entries in the naturally ordered vector, v.
1634  /// Here n is the block number in the current block preconditioner.
1635  /// If the preconditioner is a subsidiary block preconditioner
1636  /// the other entries in v that are not associated with it
1637  /// are left alone.
1638  void return_block_vector(const unsigned& n, const DoubleVector& b,
1639  DoubleVector& v) const;
1640 
1641  /// \short Given the naturally ordered vector, v, return
1642  /// the vector rearranged in block order in w. This function calls
1643  /// get_concatenated_block_vector(...) with the identity block mapping.
1644  ///
1645  /// This function has been re-written to work with the new dof type
1646  /// coarsening feature. The old function is kept alive in
1647  /// internal_get_block_ordered_preconditioner_vector(...) and is moved to
1648  /// the private section of the code. The differences between the two are:
1649  ///
1650  /// 1) This function extracts all the block vectors (in one go) via the
1651  /// function internal_get_block_vectors(...), and concatenates them.
1652  ///
1653  /// 2) The old function makes use of the variables ending in "get_ordered",
1654  /// thus is slightly more efficient since it does not have to concatenate
1655  /// any block vectors.
1656  ///
1657  /// 3) The old function no longer respect the new indirections if dof types
1658  /// have been coarsened.
1659  ///
1660  /// 4) This function extracts the most fine grain dof-level vectors and
1661  /// concatenates them. These dof-level vectors respect the re-ordering
1662  /// caused by the coarsening of dof types. The overhead associated with
1663  /// concatenating DoubleVectors without communication is very small.
1664  ///
1665  /// This function should be used.
1667  DoubleVector& w);
1668 
1669  /// \short Takes the block ordered vector, w, and reorders it in natural
1670  /// order. Reordered vector is returned in v. Note: If the preconditioner is
1671  /// a subsidiary preconditioner then only the components of the vector
1672  /// associated with the blocks of the subsidiary preconditioner will be
1673  /// included. Hence the length of v is master_nrow() whereas that of the
1674  /// vector w is of length this->nrow().
1675  ///
1676  /// This is the return function for the function
1677  /// get_block_ordered_preconditioner_vector(...).
1678  ///
1679  /// It calls the function return_concatenated_block_vector(...) with the
1680  /// identity block number ordering.
1682  DoubleVector& v) const;
1683 
1684  /// \short Return the number of block types.
1685  unsigned nblock_types() const
1686  {
1687 #ifdef PARANOID
1688  if(Block_to_dof_map_coarse.size() == 0)
1689  {
1690  std::ostringstream error_msg;
1691  error_msg <<"The Block_to_dof_map_coarse vector is not setup for \n"
1692  << "this block preconditioner.\n\n"
1693 
1694  << "This vector is always set up within block_setup(...).\n"
1695  << "If block_setup() is already called, then perhaps there is\n"
1696  << "something wrong with your block preconditionable elements.\n"
1697  << std::endl;
1698  throw OomphLibError(error_msg.str(),
1699  OOMPH_CURRENT_FUNCTION,
1700  OOMPH_EXCEPTION_LOCATION);
1701  }
1702 #endif
1703 
1704  // Return the number of block types.
1705  return Block_to_dof_map_coarse.size();
1706  } // EOFunc nblock_types(...)
1707 
1708  /// \short Return the total number of DOF types.
1709  unsigned ndof_types() const
1710  {
1711 #ifdef PARANOID
1712  // Subsidiary preconditioners don't really need the meshes
1713  if (this->is_master_block_preconditioner())
1714  {
1715  std::ostringstream err_msg;
1716  unsigned n=nmesh();
1717  if (n==0)
1718  {
1719  err_msg << "No meshes have been set for this block preconditioner!\n"
1720  << "Set one with set_nmesh(...), set_mesh(...)" << std::endl;
1721  throw OomphLibError(err_msg.str(),
1722  OOMPH_CURRENT_FUNCTION,
1723  OOMPH_EXCEPTION_LOCATION);
1724  for (unsigned m=0;m<n;m++)
1725  {
1726  if (Mesh_pt[m]==0)
1727  {
1728  err_msg << "The mesh pointer to mesh " << m << " is null!\n"
1729  << "Set a non-null one with set_mesh(...)" << std::endl;
1730  throw OomphLibError(err_msg.str(),
1731  OOMPH_CURRENT_FUNCTION,
1732  OOMPH_EXCEPTION_LOCATION);
1733 
1734  }
1735  }
1736  }
1737  }
1738 #endif
1739 
1740  // If this is a subsidiary block preconditioner, then the function
1741  // turn_into_subsidiary_block_preconditioner(...) should have been called,
1742  // this function would have set up the look up lists between coarsened
1743  // dof types and the internal dof types. Of coarse, the user (the writer
1744  // of the preconditioners) should not care about the internal dof types
1745  // and what happens under the hood. Thus they should get the number of
1746  // coarsened dof types (i.e. the number of dof types the preconditioner
1747  // above (parent preconditioner) decides to give to this preconditioner).
1749  {
1750 #ifdef PARANOID
1751  if(Doftype_coarsen_map_coarse.size() == 0)
1752  {
1753  std::ostringstream error_msg;
1754  error_msg <<"The Doftype_coarsen_map_coarse vector is not setup for \n"
1755  << "this SUBSIDIARY block preconditioner.\n\n"
1756 
1757  << "For SUBSIDARY block preconditioners at any level, this\n"
1758  << "vector is set up in the function \n"
1759  << "turn_into_subsidiary_block_preconditioner(...).\n\n"
1760 
1761  << "Being a SUBSIDIARY block preconditioner means that \n"
1762  << "(Master_block_preconditioner_pt == 0) is true.\n"
1763  << "The Master_block_preconditioner_pt MUST be set in the "
1764  << "function \n"
1765  << "turn_into_subsidiary_block_preconditioner(...).\n\n"
1766 
1767  << "Somewhere, the Master_block_preconditioner_pt pointer is\n"
1768  << "set but not by the function\n"
1769  << "turn_into_subsidiary_block_preconditioner(...).\n"
1770  << std::endl;
1771  throw OomphLibError(error_msg.str(),
1772  OOMPH_CURRENT_FUNCTION,
1773  OOMPH_EXCEPTION_LOCATION);
1774  }
1775 #endif
1776  // return the number of dof types.
1777  return Doftype_coarsen_map_coarse.size();
1778  }
1779  else
1780  // Otherwise the number of ndof types is the same as the internal number
1781  // of dof types, since no coarsening of the dof types is done at the
1782  // top-most master level.
1783  {
1784  return internal_ndof_types();
1785  }
1786  } // EOFunc ndof_types(...)
1787 
1788 
1789  /// \short Access to i-th mesh (of the various meshes that contain block
1790  /// preconditionable elements of the same number of dof type).
1791  ///
1792  /// WARNING: This should only be used if the derived class is the
1793  /// upper-most master block preconditioner. An error is thrown is
1794  /// this function is called from a subsidiary preconditioner.
1795  /// They (and since every block preconditioner can in principle
1796  /// be used as s subsidiary preconditioner: all block preconditioners)
1797  /// should store local copies of "their meshes" (if they're needed
1798  /// for anything)
1799  const Mesh* mesh_pt(const unsigned& i) const
1800  {
1801 #ifdef PARANOID
1803  {
1804  std::ostringstream error_msg;
1805  error_msg << "The mesh_pt() function should not be called on a\n"
1806  << "subsidiary block preconditioner." << std::endl;
1807  throw OomphLibError(error_msg.str(),
1808  OOMPH_CURRENT_FUNCTION,
1809  OOMPH_EXCEPTION_LOCATION);
1810  }
1811 #endif
1812 
1813  const Mesh* mesh_i_pt = Mesh_pt[i];
1814 
1815 #ifdef PARANOID
1816  if(mesh_i_pt == 0)
1817  {
1818  std::ostringstream error_msg;
1819  error_msg << "Mesh pointer " << mesh_i_pt << " is null.";
1820  throw OomphLibError(error_msg.str(),
1821  OOMPH_CURRENT_FUNCTION,
1822  OOMPH_EXCEPTION_LOCATION);
1823  }
1824 #endif
1825 
1826  return mesh_i_pt;
1827  } // EOFunc mesh_pt(...)
1828 
1829  /// \short Return the number of meshes in Mesh_pt.
1830  ///
1831  /// WARNING: This should only be used if the derived class is the
1832  /// upper-most master block preconditioner. All block preconditioners)
1833  /// should store local copies of "their meshes" (if they're needed
1834  /// for anything)
1835  unsigned nmesh() const
1836  {
1837  return Mesh_pt.size();
1838  } // EOFunc nmesh()
1839 
1840  /// \short Return the block number corresponding to a global index i_dof.
1841  int block_number(const unsigned& i_dof) const
1842  {
1843  int internal_block_number = this->internal_block_number(i_dof);
1844 
1845  if(internal_block_number == -1)
1846  {
1847  return internal_block_number;
1848  }
1849  else
1850  {
1851  // Map the internal block to the "external" block number, i.e. what the
1852  // writer of the preconditioner is expects.
1853  unsigned block_i = 0;
1854  while(std::find(Block_to_dof_map_fine[block_i].begin(),
1855  Block_to_dof_map_fine[block_i].end(),
1856  internal_block_number)
1857  == Block_to_dof_map_fine[block_i].end())
1858  {
1859  block_i++;
1860  }
1861 
1862  return block_i;
1863  }
1864  }
1865 
1866  /// \short Given a global dof number, returns the index in the block it
1867  /// belongs to.
1868  /// This is the overall index, not local block (in parallel).
1869  int index_in_block(const unsigned& i_dof) const
1870  {
1871  // the dof block number
1872  int internal_dof_block_number = this->internal_dof_number(i_dof);
1873 
1874  if(internal_dof_block_number >= 0)
1875  {
1876  // the external block number
1877  unsigned ex_blk_number = this->block_number(i_dof);
1878 
1879  int internal_index_in_dof = this->internal_index_in_dof(i_dof);
1880 
1881  // find the processor which this global index in block belongs to.
1882  unsigned block_proc
1883  = internal_block_distribution_pt(internal_dof_block_number)
1884  ->rank_of_global_row(internal_index_in_dof);
1885 
1886  // Add up all of the first rows.
1887  const unsigned ndof_in_block
1888  = Block_to_dof_map_fine[ex_blk_number].size();
1889 
1890  unsigned index = 0;
1891  for (unsigned dof_i = 0; dof_i < ndof_in_block; dof_i++)
1892  {
1894  Block_to_dof_map_fine[ex_blk_number][dof_i])
1895  ->first_row(block_proc);
1896  }
1897 
1898  // Now add up all the nrow_local up to this dof block.
1899  unsigned j = 0;
1900 
1901  while(int(Block_to_dof_map_fine[ex_blk_number][j])!= internal_dof_block_number)
1902  {
1904  Block_to_dof_map_fine[ex_blk_number][j])->nrow_local(block_proc);
1905  j++;
1906  }
1907 
1908  // Now add the index of this block...
1909  index += (internal_index_in_dof -
1911  internal_dof_block_number)
1912  ->first_row(block_proc));
1913 
1914  return index;
1915  }
1916 
1917  return -1;
1918  }
1919 
1920  /// \short Access function to the block distributions (const version).
1922  block_distribution_pt(const unsigned& b) const
1923  {
1924 #ifdef PARANOID
1925  if(Block_distribution_pt.size() == 0)
1926  {
1927  std::ostringstream error_msg;
1928  error_msg <<"Block distributions are not set up.\n"
1929  << "Have you called block_setup(...)?\n"
1930  << std::endl;
1931  throw OomphLibError(error_msg.str(),
1932  OOMPH_CURRENT_FUNCTION,
1933  OOMPH_EXCEPTION_LOCATION);
1934  }
1935  if(b > nblock_types())
1936  {
1937  std::ostringstream error_msg;
1938  error_msg <<"You requested the distribution for the block "
1939  << b << ".\n"
1940  << "But there are only " << nblock_types() << " block types.\n"
1941  << std::endl;
1942  throw OomphLibError(error_msg.str(),
1943  OOMPH_CURRENT_FUNCTION,
1944  OOMPH_EXCEPTION_LOCATION);
1945  }
1946 #endif
1947 
1948  return Block_distribution_pt[b];
1949  } // EOFunc block_distribution_pt(...)
1950 
1951  /// \short Access function to the block distributions (non-const version).
1953  block_distribution_pt(const unsigned b)
1954  {
1955 #ifdef PARANOID
1956  if(Block_distribution_pt.size() == 0)
1957  {
1958  std::ostringstream error_msg;
1959  error_msg <<"Block distributions are not set up.\n"
1960  << "Have you called block_setup(...)?\n"
1961  << std::endl;
1962  throw OomphLibError(error_msg.str(),
1963  OOMPH_CURRENT_FUNCTION,
1964  OOMPH_EXCEPTION_LOCATION);
1965  }
1966  if(b > nblock_types())
1967  {
1968  std::ostringstream error_msg;
1969  error_msg <<"You requested the distribution for the block "
1970  << b << ".\n"
1971  << "But there are only " << nblock_types() << " block types.\n"
1972  << std::endl;
1973  throw OomphLibError(error_msg.str(),
1974  OOMPH_CURRENT_FUNCTION,
1975  OOMPH_EXCEPTION_LOCATION);
1976  }
1977 #endif
1978 
1979  return Block_distribution_pt[b];
1980  } // EOFunc block_distribution_pt(...)
1981 
1982  /// \short Access function to the dof-level block distributions.
1984  dof_block_distribution_pt(const unsigned& b)
1985  {
1986 #ifdef PARANOID
1987  if(b > ndof_types())
1988  {
1989  std::ostringstream error_msg;
1990  error_msg <<"You requested the distribution for the dof block "
1991  << b << ".\n"
1992  << "But there are only " << ndof_types() << " DOF types.\n"
1993  << std::endl;
1994  throw OomphLibError(error_msg.str(),
1995  OOMPH_CURRENT_FUNCTION,
1996  OOMPH_EXCEPTION_LOCATION);
1997  }
1998 
1999 #endif
2000 
2002  {
2003 #ifdef PARANOID
2004  if(Internal_block_distribution_pt.size() == 0)
2005  {
2006  std::ostringstream error_msg;
2007  error_msg <<"Internal block distributions are not set up.\n"
2008  << "Have you called block_setup(...)?\n"
2009  << std::endl;
2010  throw OomphLibError(error_msg.str(),
2011  OOMPH_CURRENT_FUNCTION,
2012  OOMPH_EXCEPTION_LOCATION);
2013  }
2014 #endif
2015 
2016  // The dof block is distribution is the same as the internal
2017  // block distribution.
2019  }
2020  else
2021  {
2022 #ifdef PARANOID
2023  if(Dof_block_distribution_pt.size() == 0)
2024  {
2025  std::ostringstream error_msg;
2026  error_msg <<"Dof block distributions are not set up.\n"
2027  << "Have you called block_setup(...)?\n"
2028  << std::endl;
2029  throw OomphLibError(error_msg.str(),
2030  OOMPH_CURRENT_FUNCTION,
2031  OOMPH_EXCEPTION_LOCATION);
2032  }
2033 #endif
2034  return Dof_block_distribution_pt[b];
2035  }
2036  } // EOFunc block_distribution_pt(...)
2037 
2038 
2039  /// \short Access function to the distribution of the master
2040  /// preconditioner. If this preconditioner does not have a master
2041  /// preconditioner then the distribution of this preconditioner is returned.
2043  {
2045  {
2046  return this->distribution_pt();
2047  }
2048  else
2049  {
2050  return Master_block_preconditioner_pt->master_distribution_pt();
2051  }
2052  } // EOFunc master_distribution_pt(...)
2053 
2054  /// \short Return the number of DOF types in mesh i.
2055  /// WARNING: This should only be used by the upper-most master block
2056  /// preconditioner. An error is thrown is
2057  /// this function is called from a subsidiary preconditioner.
2058  /// They (and since every block preconditioner can in principle
2059  /// be used as s subsidiary preconditioner: all block preconditioners)
2060  /// should store local copies of "their meshes" (if they're needed
2061  /// for anything)
2062  unsigned ndof_types_in_mesh(const unsigned& i) const
2063  {
2064 #ifdef PARANOID
2066  {
2067  std::ostringstream err_msg;
2068  err_msg << "A subsidiary block preconditioner should not care about\n"
2069  << "anything to do with meshes.";
2070  throw OomphLibError(err_msg.str(), OOMPH_CURRENT_FUNCTION,
2071  OOMPH_EXCEPTION_LOCATION);
2072  }
2073 #endif
2074  if(Ndof_types_in_mesh.size() == 0)
2075  {
2076  return mesh_pt(i)->ndof_types();
2077  }
2078  else
2079  {
2080  return Ndof_types_in_mesh[i];
2081  }
2082  } // EOFunc ndof_types_in_mesh(...)
2083 
2084  /// \short Return true if this preconditioner is a subsidiary
2085  /// preconditioner.
2087  {
2088  return (this->Master_block_preconditioner_pt != 0);
2089  } // EOFunc is_subsidiary_block_preconditioner()
2090 
2091  /// \short Return true if this preconditioner is the master block
2092  /// preconditioner.
2094  {
2095  return (this->Master_block_preconditioner_pt == 0);
2096  } // EOFunc is_master_block_preconditioner()
2097 
2098  /// \short Set the base part of the filename to output blocks to. If it is
2099  /// set then all blocks will be output at the end of block_setup. If it is
2100  /// left empty nothing will be output.
2101  void set_block_output_to_files(const std::string& basefilename)
2102  {
2103  Output_base_filename = basefilename;
2104  } // EOFunc set_block_output_to_files(...)
2105 
2106  /// \short Turn off output of blocks (by clearing the basefilename string).
2108  {
2109  Output_base_filename.clear();
2110  } // EOFunc disable_block_output_to_files()
2111 
2112  /// \short Test if output of blocks is on or not.
2113  bool block_output_on() const
2114  {
2115  return Output_base_filename.size() > 0;
2116  } // EOFunc block_output_on()
2117 
2118  /// Output all blocks to numbered files. Called at the end of get blocks if
2119  /// an output filename has been set.
2120  void output_blocks_to_files(const std::string& basefilename,
2121  const unsigned& precision = 8) const
2122  {
2123  unsigned nblocks = internal_nblock_types();
2124 
2125  for(unsigned i=0; i<nblocks; i++)
2126  {
2127  for(unsigned j=0; j<nblocks; j++)
2128  {
2129  // Construct the filename.
2130  std::string filename(basefilename + "_block_"
2132  + "_" + StringConversion::to_string(j));
2133 
2134  // Write out the block.
2135  get_block(i,j).sparse_indexed_output(filename, precision, true);
2136  }
2137  }
2138  } // EOFunc output_blocks_to_files(...)
2139 
2140  /// \short A helper method to reduce the memory requirements of block
2141  /// preconditioners. Once the methods get_block(...), get_blocks(...)
2142  /// and build_preconditioner_matrix(...) have been called in this and
2143  /// all subsidiary block preconditioners this method can be called to
2144  /// clean up.
2146  {
2148  {
2149  Index_in_dof_block_dense.clear();
2150  Dof_number_dense.clear();
2151 #ifdef OOMPH_HAS_MPI
2152  Index_in_dof_block_sparse.clear();
2153  Dof_number_sparse.clear();
2154  Global_index_sparse.clear();
2155  Index_in_dof_block_sparse.clear();
2156  Dof_number_sparse.clear();
2157 #endif
2158  Dof_dimension.clear();
2159  }
2160  Ndof_in_block.clear();
2163  } // EOFunc post_block_matrix_assembly_partial_clear()
2164 
2165  /// \short Access function to the master block preconditioner pt.
2167  {
2168 #ifdef PARANOID
2170  {
2171  std::ostringstream error_message;
2172  error_message << "This block preconditioner does not have "
2173  << "a master preconditioner.";
2174  throw OomphLibError(error_message.str(), OOMPH_CURRENT_FUNCTION,
2175  OOMPH_EXCEPTION_LOCATION);
2176  }
2177 #endif
2179  } // EOFunc master_block_preconditioner_pt()
2180 
2181  /// \short Clears all BlockPreconditioner data. Called by the destructor
2182  /// and the block_setup(...) methods
2184  {
2185 
2186  Replacement_dof_block_pt.clear();
2187 
2188  // clear the Distributions
2189  this->clear_distribution();
2190  unsigned nblock = Internal_block_distribution_pt.size();
2191  for (unsigned b = 0; b < nblock; b++)
2192  {
2194  }
2196 
2197  // clear the global index
2198  Global_index.clear();
2199 
2200  // call the post block matrix assembly clear
2202 
2203 #ifdef OOMPH_HAS_MPI
2204  // storage if the matrix is distributed
2205  unsigned nr = Rows_to_send_for_get_block.nrow();
2206  unsigned nc = Rows_to_send_for_get_block.ncol();
2207  for (unsigned p = 0; p < nc; p++)
2208  {
2209  delete[] Rows_to_send_for_get_ordered[p];
2210  delete[] Rows_to_recv_for_get_ordered[p];
2211  for (unsigned b = 0; b < nr; b++)
2212  {
2213  delete[] Rows_to_recv_for_get_block(b,p);
2214  delete[] Rows_to_send_for_get_block(b,p);
2215  }
2216  }
2225 
2226 #endif
2227 
2228  // zero
2230  {
2231  Nrow = 0;
2232  Internal_ndof_types = 0;
2234  }
2235 
2236  // delete the prec matrix dist pt
2241 
2242  // Delete any existing (external) block distributions.
2243  const unsigned n_existing_block_dist
2244  = Block_distribution_pt.size();
2245  for (unsigned dist_i = 0; dist_i < n_existing_block_dist; dist_i++)
2246  {
2247  delete Block_distribution_pt[dist_i];
2248  }
2249 
2250  // Clear the vector.
2251  Block_distribution_pt.clear();
2252 
2253 
2254  // Create the identity key.
2255  Vector<unsigned> preconditioner_matrix_key(n_existing_block_dist,0);
2256  for (unsigned i = 0; i < n_existing_block_dist; i++)
2257  {
2258  preconditioner_matrix_key[i] = i;
2259  }
2260 
2261  // Now iterate through Auxiliary_block_distribution_pt
2262  // and delete all distributions, except for the one which corresponds
2263  // to the identity since this is already deleted.
2264  std::map<Vector<unsigned>, LinearAlgebraDistribution*>::iterator iter
2266 
2267  while(iter != Auxiliary_block_distribution_pt.end())
2268  {
2269  if(iter->first != preconditioner_matrix_key)
2270  {
2271  delete iter->second;
2272  iter++;
2273  }
2274  else
2275  {
2276  ++iter;
2277  }
2278  }
2279 
2280  // Now clear it.
2282 
2283  // Delete any dof block distributions
2284  const unsigned ndof_block_dist = Dof_block_distribution_pt.size();
2285  for (unsigned dof_i = 0; dof_i < ndof_block_dist; dof_i++)
2286  {
2287  delete Dof_block_distribution_pt[dof_i];
2288  }
2289  Dof_block_distribution_pt.clear();
2290 
2291  } // EOFunc clear_block_preconditioner_base()
2292 
2293  /// \short debugging method to document the setup.
2294  /// Should only be called after block_setup(...).
2295  void document()
2296  {
2297  oomph_info << std::endl;
2298  oomph_info << "===========================================" << std::endl;
2299  oomph_info << "Block Preconditioner Documentation" << std::endl
2300  << std::endl;
2301  oomph_info << "Number of DOF types: " << internal_ndof_types() << std::endl;
2302  oomph_info << "Number of block types: " << internal_nblock_types()
2303  << std::endl;
2304  oomph_info << std::endl;
2306  {
2307  for (unsigned d = 0; d < Internal_ndof_types; d++)
2308  {
2309  oomph_info << "Master DOF number " << d << " : "
2310  << this->internal_master_dof_number(d) << std::endl;
2311  }
2312  }
2313  oomph_info << std::endl;
2314  for (unsigned b = 0; b < internal_nblock_types(); b++)
2315  {
2316  oomph_info << "Block " << b << " DOF types:";
2317  for (unsigned i = 0; i < Block_number_to_dof_number_lookup[b].size();
2318  i++)
2319  {
2321  }
2322  oomph_info << std::endl;
2323  }
2324  oomph_info << std::endl;
2325  oomph_info << "Master block preconditioner distribution:" << std::endl;
2326  oomph_info << *master_distribution_pt() << std::endl;
2327  oomph_info << "Internal preconditioner matrix distribution:" << std::endl;
2329  oomph_info << "Preconditioner matrix distribution:" << std::endl;
2331  for (unsigned b = 0; b < Internal_nblock_types; b++)
2332  {
2333  oomph_info << "Internal block " << b << " distribution:" << std::endl;
2334  oomph_info << *Internal_block_distribution_pt[b] << std::endl;
2335  }
2336  for (unsigned b = 0; b < nblock_types(); b++)
2337  {
2338  oomph_info << "Block " << b << " distribution:" << std::endl;
2339  oomph_info << *Block_distribution_pt[b] << std::endl;
2340  }
2341 
2342  // DS: the functions called here no longer exist and this function is
2343  // never used as far as I can tell, so it should be fine to comment this
2344  // bit out:
2345  // if (is_master_block_preconditioner())
2346  // {
2347  // oomph_info << "First look-up row: " << this->first_lookup_row()
2348  // << std::endl;
2349  // oomph_info << "Number of look-up rows: "
2350  // << this->nlookup_rows() << std::endl;
2351  // }
2352  oomph_info << "===========================================" << std::endl;
2353  oomph_info << std::endl;
2354  } // EOFunc document()
2355 
2356  /// \short Access function for the Doftype_coarsen_map_fine
2357  /// variable.
2359  {
2360  return Doftype_coarsen_map_fine;
2361  }
2362 
2363  /// \short Returns the most fine grain dof types in a (possibly coarsened)
2364  /// dof type.
2366  {
2367 #ifdef PARANOID
2368  const unsigned n_dof_types = ndof_types();
2369 
2370  if(i >= n_dof_types)
2371  {
2372  std::ostringstream err_msg;
2373  err_msg << "Trying to get the most fine grain dof types in dof type "
2374  << i << ",\nbut there are only " << n_dof_types
2375  << " number of dof types.\n";
2376  throw OomphLibError(err_msg.str(),
2377  OOMPH_CURRENT_FUNCTION,
2378  OOMPH_EXCEPTION_LOCATION);
2379  }
2380 #endif
2381  return Doftype_coarsen_map_fine[i];
2382  }
2383 
2384  /// \short Access function for the number of most fine grain dof types in
2385  /// a (possibly coarsened) dof type.
2386  unsigned nfine_grain_dof_types_in(const unsigned& i) const
2387  {
2388 #ifdef PARANOID
2389  const unsigned n_dof_types = ndof_types();
2390 
2391  if(i >= n_dof_types)
2392  {
2393  std::ostringstream err_msg;
2394  err_msg << "Trying to get the number of most fine grain dof types "
2395  << "in dof type " << i << ",\n"
2396  << "but there are only " << n_dof_types
2397  << " number of dof types.\n";
2398  throw OomphLibError(err_msg.str(),
2399  OOMPH_CURRENT_FUNCTION,
2400  OOMPH_EXCEPTION_LOCATION);
2401  }
2402 #endif
2403  return Doftype_coarsen_map_fine[i].size();
2404  }
2405 
2406  /// \short Access function to the replaced dof-level blocks.
2408  {
2409  return Replacement_dof_block_pt;
2410  } // EOFunc replacement_block_pt()
2411 
2412  /// \short Setup a matrix vector product.
2413  /// matvec_prod_pt is a pointer to the MatrixVectorProduct,
2414  /// block_pt is a pointer to the block matrix,
2415  /// block_col_indices is a vector indicating which block indices does the
2416  /// RHS vector we want to multiply the matrix by.
2417  ///
2418  /// The distribution of the block column must be the same as the
2419  /// RHS vector being solved. By default, OOMPH-LIB's uniform row distribution
2420  /// is employed. When block matrices are concatenated without communication,
2421  /// the columns are permuted, as a result, the distribution of the columns
2422  /// may no longer be uniform.
2424  CRDoubleMatrix* block_pt,
2425  const Vector<unsigned>& block_col_indices)
2426  {
2427  const unsigned nblock = block_col_indices.size();
2428 
2429  if(nblock == 1)
2430  {
2431  const unsigned col_index = block_col_indices[0];
2432  matvec_prod_pt->setup(block_pt,
2433  Block_distribution_pt[col_index]);
2434  }
2435  else
2436  {
2437  std::map<Vector<unsigned>,
2438  LinearAlgebraDistribution* >::const_iterator iter;
2439 
2440  iter = Auxiliary_block_distribution_pt.find(block_col_indices);
2441  if(iter != Auxiliary_block_distribution_pt.end())
2442  {
2443  matvec_prod_pt->setup(block_pt,iter->second);
2444  }
2445  else
2446  {
2447  Vector<LinearAlgebraDistribution*> tmp_vec_dist_pt(nblock,0);
2448  for (unsigned b = 0; b < nblock; b++)
2449  {
2450  tmp_vec_dist_pt[b] = Block_distribution_pt[block_col_indices[b]];
2451  }
2452 
2455  *tmp_dist_pt);
2456  insert_auxiliary_block_distribution(block_col_indices,tmp_dist_pt);
2457  matvec_prod_pt->setup(block_pt,tmp_dist_pt);
2458  }
2459  }
2460  } // EOFunc setup_matrix_vector_product(...)
2461 
2462  /// \short Setup matrix vector product. This is simply a wrapper
2463  /// around the other setup_matrix_vector_product function.
2465  CRDoubleMatrix* block_pt,
2466  const unsigned& block_col_index)
2467  {
2468  Vector<unsigned> col_index_vector(1,block_col_index);
2469  setup_matrix_vector_product(matvec_prod_pt,block_pt,col_index_vector);
2470  } // EOFunc setup_matrix_vector_product(...)
2471 
2472 // private:
2473 
2474  /// \short Given the naturally ordered vector, v, return
2475  /// the vector rearranged in block order in w. This is a legacy function
2476  /// from the old block preconditioning framework. Kept alive in case it may
2477  /// be needed again.
2478  ///
2479  /// This uses the variables ending in "get_ordered". We no longer use this
2480  /// type of method. This function copy values from v and re-order them
2481  /// in "block order" and place them in w. Block order means that the
2482  /// values in w are the same as the concatenated block vectors.
2483  ///
2484  /// I.e. - v is naturally ordered.
2485  /// v -> s_b, v is ordered into blocks vectors
2486  /// (requires communication)
2487  /// concatenate_without_communication(s_{0,...,nblocks},w) gives w.
2488  ///
2489  /// But this function skips out the concatenation part and builds w directly
2490  /// from v.
2491  ///
2492  /// This is nice but the function is implemented in such a way that it
2493  /// always use all the (internal) blocks and concatenated with the
2494  /// identity ordering. I.e. if this preconditioner has 3 block types, then
2495  /// w will always be:
2496  /// concatenate_without_communication([s_0, s_1, s_2], w). There is easy
2497  /// way to change this.
2498  ///
2499  /// Furthermore, it does not take into account the new dof type coarsening
2500  /// feature. So this function will most likely produce the incorrect vector
2501  /// w from what the user intended. It still works, but w will be the
2502  /// concatenation of the most fine grain dof block vectors with the
2503  /// "natural" dof type ordering.
2504  ///
2505  /// This has been superseded by the function
2506  /// get_block_ordered_preconditioner_vector(...) which does the correct
2507  /// thing.
2509  DoubleVector& w) const;
2510 
2511  /// \short Takes the block ordered vector, w, and reorders it in the natural
2512  /// order. Reordered vector is returned in v. Note: If the preconditioner is
2513  /// a subsidiary preconditioner then only the components of the vector
2514  /// associated with the blocks of the subsidiary preconditioner will be
2515  /// included. Hence the length of v is master_nrow() whereas that of the
2516  /// vector w is of length this->nrow().
2517  ///
2518  /// This is the return function for the function
2519  /// internal_get_block_ordered_preconditioner_vector(...).
2520  /// Both internal_get_block_ordered_preconditioner_vector(...) and
2521  /// internal_return_block_ordered_preconditioner_vector(...) has been
2522  /// superseded by the functions
2523  ///
2524  /// get_block_ordered_preconditioner_vector(...) and
2525  /// return_block_ordered_preconditioner_vector(...),
2526  ///
2527  /// Thus this function is moved to the private section of the code.
2529  const DoubleVector& w, DoubleVector& v) const;
2530 
2531  /// \short Return the number internal blocks. This should be the same
2532  /// as the number of internal dof types. Internally, the block
2533  /// preconditioning framework always work with the most fine grain
2534  /// blocks. I.e. it always deal with the most fine grain dof-level blocks.
2535  /// This allows for coarsening of dof types. When we extract a block,
2536  /// we look at the Block_to_dof_map_fine vector to find out which most fine
2537  /// grain dof types belongs to this block.
2538  ///
2539  /// The preconditioner writer should not have to deal with internal
2540  /// dof/block types and thus this function has been moved to private.
2541  ///
2542  /// This is legacy code from before the coarsening dof type functionality
2543  /// was added. This is kept alive because it is still used in the
2544  /// internal workings of the block preconditioning framework.
2545  ///
2546  /// The function nblock_types(...) should be used if the number of block
2547  /// types is required.
2548  unsigned internal_nblock_types() const
2549  {
2550 #ifdef PARANOID
2551  if(Internal_nblock_types == 0)
2552  {
2553  std::ostringstream err_msg;
2554  err_msg <<"(Internal_nblock_types == 0) is true. \n"
2555  << "Did you remember to call the function block_setup(...)?\n\n"
2556 
2557  << "This variable is always set up within block_setup(...).\n"
2558  << "If block_setup() is already called, then perhaps there is\n"
2559  << "something wrong with your block preconditionable elements.\n"
2560  << std::endl;
2561  throw OomphLibError(err_msg.str(),
2562  OOMPH_CURRENT_FUNCTION,
2563  OOMPH_EXCEPTION_LOCATION);
2564  }
2565 
2567  {
2568  std::ostringstream err_msg;
2569  err_msg << "The number of internal block types and "
2570  << "internal dof types does not match... \n\n"
2571  << "Internally, the number of block types and the number of dof "
2572  << "types must be the same.\n"
2573  << std::endl;
2574  throw OomphLibError(err_msg.str(),
2575  OOMPH_CURRENT_FUNCTION,
2576  OOMPH_EXCEPTION_LOCATION);
2577  }
2578 #endif
2579 
2580  // return the number of internal block types.
2581  return Internal_nblock_types;
2582  } // EOFunc internal_nblock_types(...)
2583 
2584  /// \short Return the number of internal dof types. This is the number of
2585  /// most fine grain dof types. The preconditioner writer should not have to
2586  /// concern him/her-self with the internal dof/block types. Thus this fuction
2587  /// is moved to private.
2588  /// We have kept this function alive since it it still used deep within
2589  /// the inner workings of the block preconditioning framework.
2590  unsigned internal_ndof_types() const
2591  {
2593  // If this is a subsidiary block preconditioner, then the variable
2594  // Internal_ndof_types must always be set up.
2595  {
2596 #ifdef PARANOID
2597  if(Internal_ndof_types == 0)
2598  {
2599  std::ostringstream error_msg;
2600  error_msg <<"(Internal_ndof_types == 0) is true.\n"
2601  << "This means that the Master_block_preconditioner_pt pointer is\n"
2602  << "set but possibly not by the function\n"
2603  << "turn_into_subsidiary_block_preconditioner(...).\n\n"
2604 
2605  << "This goes against the block preconditioning framework "
2606  << "methodology.\n"
2607  << "Many machinery relies on the look up lists set up by the \n"
2608  << "function turn_into_subsidiary_block_preconditioner(...) \n"
2609  << "between the parent and child block preconditioners.\n"
2610  << std::endl;
2611  throw OomphLibError(error_msg.str(),
2612  OOMPH_CURRENT_FUNCTION,
2613  OOMPH_EXCEPTION_LOCATION);
2614  }
2615 #endif
2616  return Internal_ndof_types;
2617  }
2618  else
2619  // Else, this is a master block preconditioner, calculate the number of
2620  // dof types from the meshes.
2621  {
2622  unsigned ndof = 0;
2623  for (unsigned i = 0; i < nmesh(); i++)
2624  {ndof += ndof_types_in_mesh(i);}
2625  return ndof;
2626  }
2627  } // EOFunc internal_ndof_types(...)
2628 
2629  /// \short Takes the n-th block ordered vector, b, and copies its entries
2630  /// to the appropriate entries in the naturally ordered vector, v.
2631  /// Here n is the block number in the current block preconditioner.
2632  /// If the preconditioner is a subsidiary block preconditioner
2633  /// the other entries in v that are not associated with it
2634  /// are left alone.
2635  ///
2636  /// This version works with the internal block types. This is legacy code
2637  /// but is kept alive, hence moved to private. Please use the
2638  /// function "return_block_vector(...)".
2639  void internal_return_block_vector(const unsigned& n,
2640  const DoubleVector& b,
2641  DoubleVector& v) const;
2642 
2643  /// \short A helper function, takes the naturally ordered vector, v,
2644  /// and extracts the n-th block vector, b.
2645  /// Here n is the block number in the current preconditioner.
2646  /// NOTE: The ordering of the vector b is the same as the
2647  /// ordering of the block matrix from internal_get_block(...).
2649  const unsigned& n, const DoubleVector& v, DoubleVector& b) const;
2650 
2651 
2652  /// \short Takes the naturally ordered vector and
2653  /// rearranges it into a vector of sub vectors corresponding to the blocks,
2654  /// so s[b][i] contains the i-th entry in the vector associated with block b.
2655  /// The block_vec_number indicates which blocks we want.
2656  /// These blocks and vectors are those corresponding to the internal blocks.
2657  /// Note: If the preconditioner is a subsidiary preconditioner then only the
2658  /// sub-vectors associated with the blocks of the subsidiary preconditioner
2659  /// will be included. Hence the length of v is master_nrow() whereas the
2660  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
2662  const Vector<unsigned>& block_vec_number,
2663  const DoubleVector& v, Vector<DoubleVector >& s) const;
2664 
2665  /// \short A helper function, takes the naturally ordered vector and
2666  /// rearranges it into a vector of sub vectors corresponding to the blocks,
2667  /// so s[b][i] contains the i-th entry in the vector associated with block b.
2668  /// The block_vec_number indicates which blocks we want.
2669  /// These blocks and vectors are those corresponding to the internal blocks.
2670  /// Note: If the preconditioner is a subsidiary preconditioner then only the
2671  /// sub-vectors associated with the blocks of the subsidiary preconditioner
2672  /// will be included. Hence the length of v is master_nrow() whereas the
2673  /// total length of the s vectors is the sum of the Nrow of the sub vectors.
2674  /// This is simply a wrapper around the other internal_get_block_vectors(...)
2675  /// function with the identity block_vec_number vector.
2677  const DoubleVector& v, Vector<DoubleVector >& s) const;
2678 
2679  /// \short A helper function, takes the vector of block vectors, s, and
2680  /// copies its entries into the naturally ordered vector, v.
2681  /// If this is a subsidiary block preconditioner only those entries in v
2682  /// that are associated with its blocks are affected.
2684  const Vector<unsigned>& block_vec_number,
2685  const Vector<DoubleVector >& s, DoubleVector& v) const;
2686 
2687  /// \short A helper function, takes the vector of block vectors, s, and
2688  /// copies its entries into the naturally ordered vector, v.
2689  /// If this is a subsidiary block preconditioner only those entries in v
2690  /// that are associated with its blocks are affected.
2691  /// This is simple a wrapper around the other
2692  /// internal_return_block_vectors(...) function with the identity
2693  /// block_vec_number vector.
2695  const Vector<DoubleVector >& s, DoubleVector& v) const;
2696 
2697  /// \short Gets block (i,j) from the matrix pointed to by
2698  /// Matrix_pt and returns it in output_block. This is associated with the
2699  /// internal blocks. Please use the other get_block(...) function.
2700  void internal_get_block(const unsigned& i, const unsigned& j,
2701  MATRIX& output_block) const;
2702 
2703  /// \short Return the block number corresponding to a global index i_dof.
2704  /// This returns the block number corresponding to the internal blocks.
2705  /// What this means is that this returns the most fine grain dof-block
2706  /// number which this global index i_dof corresponds to. Since the writer
2707  /// of the preconditioner does not need to care about the internal block
2708  /// types, this function should not be used and thus moved to private.
2709  /// This function should not be removed since it is still used deep within
2710  /// the inner workings of the block preconditioning framework.
2711  int internal_block_number(const unsigned& i_dof) const
2712  {
2713  int dn = internal_dof_number(i_dof);
2714  if (dn == -1)
2715  {
2716  return dn;
2717  }
2718  else
2719  {
2721  }
2722  } // EOFunc internal_block_number(...)
2723 
2724  /// \short Return the index in the block corresponding to a global block
2725  /// number i_dof. The index returned corresponds to the internal blocks,
2726  /// which is the most fine grain dof blocks.
2727  int internal_index_in_block(const unsigned& i_dof) const
2728  {
2729  // the index in the dof block
2730  unsigned index = internal_index_in_dof(i_dof);
2731 
2732  // the dof block number
2733  int internal_dof_block_number = internal_dof_number(i_dof);
2734  if (internal_dof_block_number >= 0)
2735  {
2736 
2737  // the 'actual' block number
2738  unsigned blk_number = internal_block_number(i_dof);
2739 
2740  // compute the index in the block
2741  unsigned j = 0;
2742  while (int(Block_number_to_dof_number_lookup[blk_number][j]) !=
2743  internal_dof_block_number)
2744  {
2745  index +=
2747  (Block_number_to_dof_number_lookup[blk_number][j]);
2748  j++;
2749  }
2750 
2751  // and return
2752  return index;
2753  }
2754  return -1;
2755  } // EOFunc internal_index_in_block(...)
2756 
2757  /// \short Access function to the internal block distributions.
2759  internal_block_distribution_pt(const unsigned& b) const
2760  {
2761 #ifdef PARANOID
2762  if(Internal_block_distribution_pt.size() == 0)
2763  {
2764  std::ostringstream error_msg;
2765  error_msg <<"Internal block distributions are not set up.\n"
2766  << "Have you called block_setup(...)?\n"
2767  << std::endl;
2768  throw OomphLibError(error_msg.str(),
2769  OOMPH_CURRENT_FUNCTION,
2770  OOMPH_EXCEPTION_LOCATION);
2771  }
2772  if(b > internal_nblock_types())
2773  {
2774  std::ostringstream error_msg;
2775  error_msg <<"You requested the distribution for the internal block "
2776  << b << ".\n" << "But there are only "
2778  << " block types.\n" << std::endl;
2779  throw OomphLibError(error_msg.str(),
2780  OOMPH_CURRENT_FUNCTION,
2781  OOMPH_EXCEPTION_LOCATION);
2782  }
2783 #endif
2785  } // EOFunc internal_block_distribution_pt(...)
2786 
2787  /// \short insert a Vector<unsigned> and LinearAlgebraDistribution* pair
2788  /// into Auxiliary_block_distribution_pt. The
2789  /// Auxiliary_block_distribution_pt should only contain pointers to
2790  /// distributions concatenated at this block level. We try to ensure this by
2791  /// checking if the block_vec_number vector is within the range
2792  /// nblock_types(). Of course, this does not guarantee correctness, but this
2793  /// is the least we can do.
2795  const Vector<unsigned>& block_vec_number,
2796  LinearAlgebraDistribution* dist_pt)
2797  {
2798  #ifdef PARANOID
2799  const unsigned max_block_number
2800  = *std::max_element(block_vec_number.begin(),
2801  block_vec_number.end());
2802 
2803  const unsigned nblocks = nblock_types();
2804  if(max_block_number >= nblocks)
2805  {
2806  std::ostringstream err_msg;
2807  err_msg << "Cannot insert into Auxiliary_block_distribution_pt\n"
2808  << "because " << max_block_number << " is equal to or \n"
2809  << "greater than " << nblocks << ".\n";
2810  throw OomphLibError(err_msg.str(),
2811  OOMPH_CURRENT_FUNCTION,
2812  OOMPH_EXCEPTION_LOCATION);
2813  }
2814 
2815  // Now check if the pair already exists in Auxiliary_block_distribution_pt.
2816  // This is a stricter test and can be removed if required.
2817 
2818  // Attempt to get an iterator pointing to the pair with the value
2819  // block_vec_number.
2820  std::map<Vector<unsigned>,
2821  LinearAlgebraDistribution* >::const_iterator iter
2822  = Auxiliary_block_distribution_pt.find(block_vec_number);
2823 
2824  if(iter != Auxiliary_block_distribution_pt.end())
2825  // If it exists, we throw an error
2826  {
2827  std::ostringstream err_msg;
2828  err_msg << "Cannot insert into Auxiliary_block_distribution_pt\n"
2829  << "because the first in the pair already exists.\n";
2830  throw OomphLibError(err_msg.str(),
2831  OOMPH_CURRENT_FUNCTION,
2832  OOMPH_EXCEPTION_LOCATION);
2833  }
2834  #endif
2835 
2837  std::make_pair(block_vec_number,
2838  dist_pt));
2839  } // insert_auxiliary_block_distribution(...)
2840 
2841  /// \short Private helper function to check that every element in the block
2842  /// matrix (i,j) matches the corresponding element in the original matrix
2843  void block_matrix_test(const unsigned& i,
2844  const unsigned& j,
2845  const MATRIX* block_matrix_pt) const;
2846 
2847  /// \short Get the index of first occurrence of value in a vector.
2848  /// If the element does not exist, -1 is returned.
2849  /// The optional parameter indicates of the Vector is sorted or not.
2850  /// Complexity: if the Vector is sorted, then on average, logarithmic in the
2851  /// distance between first and last: Performs approximately log2(N)+2
2852  /// element comparisons.
2853  /// Otherwise, up to linear in the distance between first and last:
2854  /// Compares elements until a match is found.
2855  template<typename myType>
2856  inline int get_index_of_value(const Vector<myType>& vec,
2857  const myType val,
2858  const bool sorted = false) const
2859  {
2860  if(sorted)
2861  {
2862  typename Vector<myType>::const_iterator low
2863  = std::lower_bound(vec.begin(),vec.end(),val);
2864 
2865  return (low == vec.end() || *low != val) ? -1 : (low - vec.begin());
2866  }
2867  else
2868  {
2869  int pos = std::find(vec.begin(),vec.end(),val) - vec.begin();
2870  return (pos < int(vec.size()) && pos >=0) ? pos : -1;
2871  }
2872  }
2873 
2874  private:
2875 
2876  protected:
2877 
2878  /// \short Specify the number of meshes required by this block
2879  /// preconditioner.
2880  /// Note: elements in different meshes correspond to different types
2881  /// of DOF.
2882  void set_nmesh(const unsigned& n)
2883  {
2884  Mesh_pt.resize(n,0);
2886  } // EOFunc set_nmesh(...)
2887 
2888 
2889  /// \short Set the i-th mesh for this block preconditioner.
2890  /// Note:
2891  /// The method set_nmesh(...) must be called before this method
2892  /// to specify the number of meshes.
2893  /// By default, it is assumed that each mesh only contains elements of the
2894  /// same type. This condition may be relaxed by setting the boolean
2895  /// allow_multiple_element_type_in_mesh to true, however, each mesh must only
2896  /// contain elements with the same number of dof types.
2897  void set_mesh(const unsigned& i, const Mesh* const mesh_pt,
2898  const bool &allow_multiple_element_type_in_mesh = false)
2899  {
2900 #ifdef PARANOID
2901  // paranoid check that mesh i can be set
2902  if (i >= nmesh())
2903  {
2904  std::ostringstream err_msg;
2905  err_msg
2906  << "The mesh pointer has space for " << nmesh()
2907  << " meshes.\n" << "Cannot store a mesh at entry " << i << "\n"
2908  << "Has set_nmesh(...) been called?";
2909  throw OomphLibError(err_msg.str(),
2910  OOMPH_CURRENT_FUNCTION,
2911  OOMPH_EXCEPTION_LOCATION);
2912  }
2913 
2914  // Check that the mesh pointer is not null.
2915  if(mesh_pt == 0)
2916  {
2917  std::ostringstream err_msg;
2918  err_msg
2919  << "Tried to set the " << i << "-th mesh pointer, but it is null.";
2920  throw OomphLibError(err_msg.str(),
2921  OOMPH_CURRENT_FUNCTION,
2922  OOMPH_EXCEPTION_LOCATION);
2923  }
2924 #endif
2925 
2926  // store the mesh pt and n dof types
2927  Mesh_pt[i]=mesh_pt;
2928 
2929  // Does this mesh contain multiple element types?
2931  = unsigned(allow_multiple_element_type_in_mesh);
2932  } // EOFunc set_mesh(...)
2933 
2934 
2935  /// \short Set replacement dof-level blocks.
2936  /// Only dof-level blocks can be set. This is important due to how the
2937  /// dof type coarsening feature operates.
2938  ///
2939  /// IMPORTANT: The block indices (block_i, block_j) is the dof-level
2940  /// ordering, NOT the block-ordering. The block-ordering is determined by
2941  /// the parameters given to block_setup(...).
2942  /// The DOF-ordering is determined by the two-level ordering scheme of
2943  /// first the elements, then the meshes.
2944  void set_replacement_dof_block(const unsigned &block_i,
2945  const unsigned &block_j,
2947  {
2948 #ifdef PARANOID
2949  // Check if block_setup(...) has been called.
2950  if(nblock_types() == 0)
2951  {
2952  std::ostringstream err_msg;
2953  err_msg << "nblock_types() is 0, has block_setup(...) been called?\n";
2954  throw OomphLibError(err_msg.str(),
2955  OOMPH_CURRENT_FUNCTION,
2956  OOMPH_EXCEPTION_LOCATION);
2957  }
2958 
2959 
2960  // Range checking for replacement dof block.
2961  unsigned para_ndof_types = this->ndof_types();
2962 
2963  if((block_i >= para_ndof_types) ||
2964  (block_j >= para_ndof_types))
2965  {
2966  std::ostringstream err_msg;
2967  err_msg << "Replacement dof block ("
2968  << block_i << "," << block_j << ") is outside of range:\n"
2969  << para_ndof_types;
2970  throw OomphLibError(err_msg.str(),
2971  OOMPH_CURRENT_FUNCTION,
2972  OOMPH_EXCEPTION_LOCATION);
2973  }
2974 
2975 
2976 // // Check that the most fine grain mapping has been used in block_setup(...)
2977 // // i.e. nblock_types() == ndof_types()
2978 // if(ndof_types() != nblock_types())
2979 // {
2980 // std::ostringstream err_msg;
2981 // err_msg << "ndof_types() != nblock_types()\n"
2982 // << "Only the dof-level blocks can be replaced.\n"
2983 // << "Please re-think your blocking scheme.\n";
2984 // throw OomphLibError(err_msg.str(),
2985 // OOMPH_CURRENT_FUNCTION,
2986 // OOMPH_EXCEPTION_LOCATION);
2987 // }
2988 
2989  // Check that the replacement block pt is not null
2990  if(replacement_dof_block_pt == 0)
2991  {
2992  std::ostringstream err_msg;
2993  err_msg << "Replacing block(" << block_i << "," << block_i << ")\n"
2994  << " but the pointer is NULL." << std::endl;
2995  throw OomphLibError(err_msg.str(),
2996  OOMPH_CURRENT_FUNCTION,
2997  OOMPH_EXCEPTION_LOCATION);
2998  }
2999 
3000  // Check that the replacement block has been built
3001  if(!replacement_dof_block_pt->built())
3002  {
3003  std::ostringstream err_msg;
3004  err_msg << "Replacement block(" << block_i << "," << block_i << ")"
3005  << " is not built." << std::endl;
3006  throw OomphLibError(err_msg.str(),
3007  OOMPH_CURRENT_FUNCTION,
3008  OOMPH_EXCEPTION_LOCATION);
3009  }
3010 
3011  // Check if the distribution matches. Determine which natural ordering dof
3012  // this should go to. I.e. we convert from dof-block index to dof index.
3013  // Luckily, this is stored in Block_to_dof_map_coarse.
3014 // const unsigned para_dof_block_i = Block_to_dof_map_coarse[block_i][0];
3015  const unsigned para_dof_block_i = block_i;
3016 
3017  if(*dof_block_distribution_pt(para_dof_block_i) !=
3018  *replacement_dof_block_pt->distribution_pt() )
3019  {
3020  std::ostringstream err_msg;
3021  err_msg << "The distribution of the replacement dof_block_pt\n"
3022  << "is different from the Dof_block_distribution_pt["
3023  << para_dof_block_i<<"].\n";
3024  throw OomphLibError(err_msg.str(),
3025  OOMPH_CURRENT_FUNCTION,
3026  OOMPH_EXCEPTION_LOCATION);
3027  }
3028 
3029  // Now that we know the distribution of the replacement block is
3030  // correct, we check the number of columns.
3031  const unsigned para_dof_block_j = block_j;
3032  unsigned para_replacement_block_ncol = replacement_dof_block_pt->ncol();
3033  unsigned para_required_ncol
3034  = dof_block_distribution_pt(para_dof_block_j)->nrow();
3035  if(para_replacement_block_ncol != para_required_ncol)
3036  {
3037  std::ostringstream err_msg;
3038  err_msg << "Replacement dof block has ncol = "
3039  << para_replacement_block_ncol << ".\n"
3040  << "But required ncol is " << para_required_ncol << ".\n";
3041  throw OomphLibError(err_msg.str(),
3042  OOMPH_CURRENT_FUNCTION,
3043  OOMPH_EXCEPTION_LOCATION);
3044  }
3045 #endif
3046 
3047  // Block_to_dof_map_coarse[x][0] sense because we only can use this if
3048  // nblock_types() == ndof_types(), i.e. each sub-vector is of length 1.
3049  //
3050  // We use this indirection so that the placement of the pointer is
3051  // consistent with internal_get_block(...).
3052 // const unsigned dof_block_i = Block_to_dof_map_coarse[block_i][0];
3053 // const unsigned dof_block_j = Block_to_dof_map_coarse[block_j][0];
3054 
3055 // Replacement_dof_block_pt(dof_block_i,dof_block_j)
3056 // = replacement_dof_block_pt;
3057 
3058  Replacement_dof_block_pt(block_i,block_j)
3060  }
3061 
3062  /// \short Check if any of the meshes are distributed. This is equivalent
3063  /// to problem.distributed() and is used as a replacement.
3065  {
3066 #ifdef OOMPH_HAS_MPI
3067  // is_mesh_distributed() is only available with MPI
3068  for(unsigned i=0, n=nmesh(); i<n; i++)
3069  {
3070  if(mesh_pt(i)->is_mesh_distributed()) { return true; }
3071  }
3072 #endif
3073  return false;
3074  }
3075 
3076  /// \short Return the number of the block associated with global unknown
3077  /// i_dof. If this preconditioner is a subsidiary block preconditioner then
3078  /// the block number in the subsidiary block preconditioner is returned. If
3079  /// a particular global DOF is not associated with this preconditioner then
3080  /// -1 is returned
3081  int internal_dof_number(const unsigned& i_dof) const
3082  {
3083 
3085  {
3086 #ifdef OOMPH_HAS_MPI
3087  unsigned first_row = this->distribution_pt()->first_row();
3088  unsigned nrow_local = this->distribution_pt()->nrow_local();
3089  unsigned last_row = first_row+nrow_local-1;
3090  if (i_dof >= first_row && i_dof <= last_row)
3091  {
3092  return static_cast<int>(Dof_number_dense[i_dof-first_row]);
3093  }
3094  else
3095  {
3096  //int index = this->get_index_of_element(Global_index_sparse,i_dof);
3097  int index = get_index_of_value<unsigned>(Global_index_sparse,i_dof,true);
3098  if (index >= 0)
3099  {
3100  return Dof_number_sparse[index];
3101  }
3102  }
3103  // if we here we couldn't find the i_dof
3104 #ifdef PARANOID
3105  unsigned my_rank = comm_pt()->my_rank();
3106  std::ostringstream error_message;
3107  error_message
3108  << "Proc " << my_rank<<": Requested internal_dof_number(...) for global DOF "
3109  << i_dof << "\n"
3110  << "cannot be found.\n";
3111  throw OomphLibError(
3112  error_message.str(),
3113  OOMPH_CURRENT_FUNCTION,
3114  OOMPH_EXCEPTION_LOCATION);
3115 #endif
3116 #else
3117  return static_cast<int>(Dof_number_dense[i_dof]);
3118 #endif
3119  }
3120  // else this preconditioner is a subsidiary one, and its Block_number
3121  // lookup schemes etc haven't been set up
3122  else
3123  {
3124  // Block number in master prec
3125  unsigned blk_num = Master_block_preconditioner_pt->internal_dof_number(i_dof);
3126 
3127  // Search through the Block_number_in_master_preconditioner for master
3128  // block blk_num and return the block number in this preconditioner
3129  for (unsigned i = 0; i < this->internal_ndof_types(); i++)
3130  {
3131  if (Doftype_in_master_preconditioner_fine[i] == blk_num)
3132  {return static_cast<int>(i);}
3133  }
3134  // if the master block preconditioner number is not found return -1
3135  return -1;
3136  }
3137 
3138  // Shouldn't get here
3139  throw OomphLibError("Never get here\n",
3140  OOMPH_CURRENT_FUNCTION,
3141  OOMPH_EXCEPTION_LOCATION);
3142  // Dummy return
3143  return -1;
3144  }
3145 
3146  /// \short Return the row/column number of global unknown i_dof within it's
3147  /// block.
3148  unsigned internal_index_in_dof(const unsigned& i_dof) const
3149  {
3151  {
3152 #ifdef OOMPH_HAS_MPI
3153  unsigned first_row = this->distribution_pt()->first_row();
3154  unsigned nrow_local = this->distribution_pt()->nrow_local();
3155  unsigned last_row = first_row+nrow_local-1;
3156  if (i_dof >= first_row && i_dof <= last_row)
3157  {
3158  return static_cast<int>(Index_in_dof_block_dense[i_dof-first_row]);
3159  }
3160  else
3161  {
3162  //int index = this->get_index_of_element(Global_index_sparse,i_dof);
3163  int index = get_index_of_value<unsigned>(Global_index_sparse,i_dof,true);
3164  if (index >= 0)
3165  {
3166  return Index_in_dof_block_sparse[index];
3167  }
3168  }
3169  // if we here we couldn't find the i_dof
3170 #ifdef PARANOID
3171  std::ostringstream error_message;
3172  error_message
3173  << "Requested internal_index_in_dof(...) for global DOF " << i_dof << "\n"
3174  << "cannot be found.\n";
3175  throw OomphLibError(
3176  error_message.str(),
3177  OOMPH_CURRENT_FUNCTION,
3178  OOMPH_EXCEPTION_LOCATION);
3179 #endif
3180 #else
3181  return Index_in_dof_block_dense[i_dof];
3182 #endif
3183  }
3184  else
3185  {
3186  return Master_block_preconditioner_pt->internal_index_in_dof(i_dof);
3187  }
3188 
3189  // Shouldn't get here
3190  throw OomphLibError("Never get here\n",
3191  OOMPH_CURRENT_FUNCTION,
3192  OOMPH_EXCEPTION_LOCATION);
3193  // Dummy return
3194  return -1;
3195  }
3196 
3197  /// \short Return the number of degrees of freedom in block b. Note that if
3198  /// this preconditioner acts as a subsidiary preconditioner then b refers
3199  /// to the block number in the subsidiary preconditioner not the master
3200  /// block preconditioner.
3201  unsigned internal_block_dimension(const unsigned& b) const
3202  {
3203 #ifdef PARANOID
3204  const unsigned i_nblock_types = internal_nblock_types();
3205  if(b >= i_nblock_types)
3206  {
3207  std::ostringstream err_msg;
3208  err_msg << "Trying to get internal block dimension for \n"
3209  << "internal block " << b <<".\n"
3210  << "But there are only " << i_nblock_types
3211  << " internal dof types.\n";
3212  throw OomphLibError(err_msg.str(),
3213  OOMPH_CURRENT_FUNCTION,
3214  OOMPH_EXCEPTION_LOCATION);
3215  }
3216 #endif
3217  return Internal_block_distribution_pt[b]->nrow();
3218  }
3219 
3220  /// \short Return the size of the dof "block" i, i.e. how many degrees of
3221  /// freedom are associated with it. Note that if this preconditioner acts as
3222  /// a subsidiary preconditioner, then i refers to the block number in the
3223  /// subsidiary preconditioner not the master block preconditioner
3224  unsigned internal_dof_block_dimension(const unsigned& i) const
3225  {
3226 #ifdef PARANOID
3227  const unsigned i_n_dof_types = internal_ndof_types();
3228  if(i >= i_n_dof_types)
3229  {
3230  std::ostringstream err_msg;
3231  err_msg << "Trying to get internal dof block dimension for \n"
3232  << "internal dof block " << i <<".\n"
3233  << "But there are only " << i_n_dof_types
3234  << " internal dof types.\n";
3235  throw OomphLibError(err_msg.str(),
3236  OOMPH_CURRENT_FUNCTION,
3237  OOMPH_EXCEPTION_LOCATION);
3238  }
3239 #endif
3240  // I don't understand the difference between this function and
3241  // block_dimension(...) but I'm not going to mess with it... David
3242 
3244  {
3245  return Dof_dimension[i];
3246  }
3247  else
3248  {
3249  unsigned master_i = internal_master_dof_number(i);
3250  return Master_block_preconditioner_pt->internal_dof_block_dimension(master_i);
3251  }
3252  }
3253 
3254  /// \short Return the number of dofs (number of rows or columns) in the
3255  /// overall problem. The prefix "master_" is sort of redundant when used as
3256  /// a stand-alone block preconditioner but is required to avoid ambiguities.
3257  /// The latter is stored (and maintained) separately for each specific block
3258  /// preconditioner regardless of its role.
3259  unsigned master_nrow() const
3260  {
3262  {
3263  return Nrow;
3264  }
3265  else
3266  {
3267  return (this->Master_block_preconditioner_pt->master_nrow());
3268  }
3269  }
3270 
3271  /// \short Takes the block number within this preconditioner and returns the
3272  /// corresponding block number in the master preconditioner. If this
3273  /// preconditioner does not have a master block preconditioner then the
3274  /// block number passed is returned
3275  unsigned internal_master_dof_number(const unsigned& b) const
3276  {
3278  return b;
3279  else
3281  }
3282 
3283  /// \short access function to the internal
3284  /// preconditioner matrix distribution pt.
3285  /// preconditioner_matrix_distribution_pt always returns the concatenation
3286  /// of the internal block distributions. Since the writer of the
3287  /// preconditioner does not need to concern themselves with the internal
3288  /// dof/block, please use preconditioner_matrix_distribution_pt().
3291  {
3294  else
3295  return this->distribution_pt();
3296  }
3297 
3298  /// \short Access function to the preconditioner matrix distribution pointer.
3299  /// This is the concatenation of the block distributions with the identity
3300  /// ordering. I.e. if this preconditioner has three block types, with the
3301  /// three associated block distributions dist_b0, dist_b1 and dist_b2, then
3302  /// this distribution is:
3303  /// LinearAlgebraDistributionHelpers::concatenate(dist_b0, dist_b1, dist_b2).
3306  {
3308  }
3309 
3310  /// \short The replacement dof-level blocks.
3312 
3313  /// \short The distribution for the blocks.
3315 
3316  /// \short Mapping for block types to dof types. These are the dof types
3317  /// the writer of the preconditioner expects. For the upper-most master
3318  /// block preconditioner, this would be the sum of the dof types in the
3319  /// meshes. For subsidiary block preconditioners, this is determined by
3320  /// the parent preconditioner when passing in the doftype_coarsen_map_coarse
3321  /// vector in turn_into_subsidiary_block_preconditioner(...).
3323 
3324  /// \short Mapping for the block types to the most fine grain dof types.
3326 
3327  /// \short Mapping for dof types within THIS precondition. This is usually
3328  /// passed down from the parent preconditioner.
3329  /// This list is used to tell which does types should
3330  /// be considered as a single dof type within this preconditioner. I.e. we
3331  /// "coarsen" the dof types. The values are local to this preconditioner,
3332  /// for example, even if the
3333  /// Doftype_in_master_preconditioner_coarse = [2,3,4], the vector
3334  /// Doftype_coarsen_map_coarse = [[0],[1,2]], saying your local dof types
3335  /// 0 should be considered as dof type 0 and dof types 1 and 2 are considered
3336  /// as dof type 1.
3337  ///
3338  /// Furthermore, the dof types are that the preconditioner above this one
3339  /// knows; these dof types may or may not be coarsened. For example, say that
3340  /// this preconditioner expects two dof types, 0 and 1.
3341  /// The preconditioner above this one wishes to use this preconditioner to
3342  /// solve the block associated with it's dof types 2, 3 and 4. It passes the
3343  /// Vector [2,3,4] to this preconditioner via the function
3344  /// turn_into_subsidiary_block_preconditioner(...), this list is to be stored
3345  /// in Doftype_in_master_preconditioner_coarse. It also passes in the
3346  /// 2D vector [[0][1,2]] (as described above), this list is to be stored in
3347  /// Doftype_coarsen_map_coarse. BUT, the master's preconditioner dof types
3348  /// may also be coarsened. I.e. the underlying dof types of the master block
3349  /// preconditioner may be [0,1,2,3,4,5,6,7], for which it may have the
3350  /// Doftype_coarsen_map_coarse = [[0,1][2,3][4,5][6,7]].
3351  ///
3352  /// An additional list has to be kept for the most fine grain dof type
3353  /// mapping. This is stored in Doftype_coarsen_map_fine, in this case it
3354  /// would be:
3355  ///
3356  /// Doftype_coarsen_map_fine = [[0,1][2,3,4,5]], since the dof types passed
3357  /// to this preconditioner is [2, 3, 4] from the master preconditioner, but
3358  /// it actually refers to the underlying dof types [2,3,4,5,6,7].
3359  ///
3360  /// In the case of the top most master block
3361  /// preconditioner, the block_setup(...) function fills the vector with the
3362  /// identity mapping.
3364 
3365  /// \short Mapping the dof types within this preconditioner. The values in
3366  /// here refers to the most grain dof types. This list is automatically
3367  /// generated either in block_setup(...) (for the top-most preconditioner)
3368  /// or the turn_into_subsidiary_block_preconditioner(...) function.
3369  /// Please refer to the comment above Doftype_coarsen_map_coarse for more
3370  /// details.
3372 
3373  /// \short Storage for the default distribution for each internal block.
3375 
3376  /// \short Storage for the default distribution for each dof block at
3377  /// this level.
3379 
3380  /// \short Vector of unsigned to indicate which meshes contain multiple
3381  /// element types.
3383 
3384  /// \short Vector of pointers to the meshes containing the elements used in
3385  /// the block preconditioner. Const pointers to prevent modification of the
3386  /// mesh by the preconditioner (this could be relaxed if needed). If this is
3387  /// a subsidiary preconditioner, then the information is looked up in the
3388  /// master preconditioner.
3390 
3391  /// \short Storage for number of types of degree of freedom of the elements
3392  /// in each mesh.
3394 
3395  /// \short Number of different block types in this preconditioner. Note that
3396  /// this information is maintained if used as a subsidiary or stand-alone
3397  /// block preconditioner, in the latter case it stores the number of blocks
3398  /// within the subsidiary preconditioner.
3400 
3401  ///\short Number of different DOF types in this preconditioner. Note that
3402  /// this information is maintained if used as a subsidiary or stand-alone
3403  /// block preconditioner, in the latter case it stores the number of dofs
3404  /// within the subsidiary preconditioner.
3406 
3407  private:
3408 
3409  /// \short Debugging variable. Set true or false via the access functions
3410  /// turn_on_recursive_debug_flag(...)
3411  /// turn_off_recursive_debug_flag(...)
3412  /// These will turn on/off the debug flag up the hierarchy.
3414 
3415  /// \short Debugging variable. Set true or false via the access functions
3416  /// turn_on_debug_flag(...)
3417  /// turn_off_debug_flag(...)
3419 
3420  /// \short Stores any block-level distributions / concatenation of
3421  /// block-level distributions required. The first in the pair
3422  /// (Vector<unsigned>) represents the block numbers of the distributions
3423  /// concatenated to get the second in the pair (LinearAlgebraDistribution*).
3424  std::map<Vector<unsigned>,LinearAlgebraDistribution*>
3426 
3427  /// \short Number of DOFs (# of rows or columns in the matrix) in this
3428  /// preconditioner. Note that this information is maintained if used as a
3429  /// subsidiary or stand-alone block preconditioner, in the latter case it
3430  /// stores the number of rows within the subsidiary preconditioner.
3431  unsigned Nrow;
3432 
3433  /// \short If the block preconditioner is acting a subsidiary block
3434  /// preconditioner then a pointer to the master preconditioner is stored
3435  /// here. If the preconditioner does not have a master block preconditioner
3436  /// then this pointer remains null.
3438 
3439  /// \short The map between the dof types in this preconditioner and the master
3440  /// preconditioner. If there is no master preconditioner it remains empty.
3441  /// This list contains the mapping for the underlying dof types.
3443 
3444  /// \short The map between the dof types in this preconditioner and the
3445  /// master preconditioner. If there is no master preconditioner, it remains
3446  /// empty. This is the version for which the master preconditioner expects.
3447  /// The dof types in here may or may not be coarsened in the preconditioner
3448  /// above this one.
3450 
3451  /// \short **This was uncommented** Presumably a non-distributed analogue of
3452  /// Index_in_dof_block_sparse.
3454 
3455  /// \short Vector to store the mapping from the global DOF number to its
3456  /// block. Empty if this preconditioner has a master preconditioner, in this
3457  /// case the information is obtained from the master preconditioner.
3459 
3460 #ifdef OOMPH_HAS_MPI
3461 
3462  // The following three vectors store data on the matrix rows/matrix
3463  // columns/dofs (the three are equivalent) that are not on this processor.
3464 
3465  /// \short For global indices outside of the range this->first_row()
3466  /// to this->first_row()+this->nrow_local(), the Index_in_dof_block
3467  /// and Dof_number are stored sparsely in the vectors:
3468  /// + Index_in_dof_block_sparse;
3469  /// + Dof_number_sparse;
3470  /// The corresponding global indices are stored in this vector.
3472 
3473  /// \short Vector to store the mapping from the global DOF number to the
3474  /// index (row/column number) within its block (empty if this preconditioner
3475  /// has a master preconditioner as this information is obtained from the
3476  /// master preconditioner). Sparse version: for global indices outside of the
3477  /// range this->first_row() to this->first_row()+this->nrow_local(). The
3478  /// global index of an element in this vector is defined in
3479  /// Global_index_sparse.
3481 
3482  /// \short Vector to store the mapping from the global DOF number to its
3483  /// block (empty if this preconditioner has a master preconditioner as this
3484  /// information is obtained from the master preconditioner). Sparse
3485  /// version: for global indices outside of the range this->first_row() to
3486  /// this->first_row()+this->nrow_local(). The global index of an element in
3487  /// this vector is defined in Global_index_sparse.
3489 #endif
3490 
3491  /// \short Vector containing the size of each block, i.e. the number of
3492  /// global DOFs associated with it. (Empty if this preconditioner has a
3493  /// master preconditioner as this information is obtain from the master
3494  /// preconditioner.)
3496 
3497  /// \short Vectors of vectors for the mapping from block number and block
3498  /// row to global row number. Empty if this preconditioner has a master
3499  /// preconditioner as this information is obtain from the master
3500  /// preconditioner.
3502 
3503  /// \short Vector of vectors to store the mapping from block number to the
3504  /// DOF number (each element could be a vector because we allow multiple
3505  /// DOFs types in a single block).
3507 
3508  /// \short Vector to the mapping from DOF number to block number.
3510 
3511  /// \short Number of types of degree of freedom associated with each block.
3513 
3514 
3515 
3516 #ifdef OOMPH_HAS_MPI
3517  /// \short The global rows to be sent of block b to processor p (matrix
3518  /// indexed [b][p]).
3520 
3521  /// \short The number of global rows to be sent of block b to processor p
3522  /// (matrix indexed [b][p]).
3524 
3525  /// \short The block rows to be received from processor p for block b
3526  /// (matrix indexed [b][p]).
3528 
3529  /// \short The number of block rows to be received from processor p for
3530  /// block b (matrix indexed [b][p]).
3532 
3533  /// \short The global rows to be sent to processor p for
3534  /// get_block_ordered_... type methods.
3536 
3537  /// \short The number global rows to be sent to processor p for
3538  /// get_block_ordered_... type methods.
3540 
3541  /// \short The preconditioner rows to be received from processor p for
3542  /// get_block_ordered_... type methods.
3544 
3545  /// \short The number of preconditioner rows to be received from processor
3546  /// p for get_block_ordered_... type methods.
3548 #endif
3549 
3550  /// \short The distribution of the (internal) preconditioner matrix. This is
3551  /// formed by concatenating the distribution of the internal blocks.
3552  /// This is obsolete code, maintained for backwards compatibility.
3553  /// Below is the old comment:
3554  ///
3555  /// - only used if this preconditioner is a master preconditioner.
3556  /// Warning: always use the access function
3557  /// internal_preconditioner_matrix_distribution_pt().
3559 
3560  /// \short The distribution of the preconditioner matrix. This is the
3561  /// concatenation of the block distribution.
3563 
3564  /// \short Static boolean to allow block_matrix_test(...) to be run.
3565  /// Defaults to false.
3567 
3568  /// \short String giving the base of the files to write block data into. If
3569  /// empty then do not output blocks. Default is empty.
3571  };
3572 
3573 
3574 }
3575 #endif
void output_blocks_to_files(const std::string &basefilename, const unsigned &precision=8) const
virtual void set_matrix_pt(DoubleMatrixBase *matrix_pt)
Set the matrix pointer.
void set_replacement_block_pt(CRDoubleMatrix *replacement_block_pt)
set Replacement_block_pt.
unsigned nfine_grain_dof_types_in(const unsigned &i) const
Access function for the number of most fine grain dof types in a (possibly coarsened) dof type...
unsigned Column_index
Column index of the block.
Vector< unsigned > Global_index_sparse
For global indices outside of the range this->first_row() to this->first_row()+this->nrow_local(), the Index_in_dof_block and Dof_number are stored sparsely in the vectors:
unsigned nrow_local() const
access function for the num of local rows on this processor.
Vector< int * > Rows_to_send_for_get_ordered
The global rows to be sent to processor p for get_block_ordered_... type methods. ...
void clear_block_preconditioner_base()
Clears all BlockPreconditioner data. Called by the destructor and the block_setup(...) methods.
Vector< const Mesh * > Mesh_pt
Vector of pointers to the meshes containing the elements used in the block preconditioner. Const pointers to prevent modification of the mesh by the preconditioner (this could be relaxed if needed). If this is a subsidiary preconditioner, then the information is looked up in the master preconditioner.
void broken_copy(const std::string &class_name)
Issue error message and terminate execution.
Vector< unsigned > Nrows_to_recv_for_get_ordered
The number of preconditioner rows to be received from processor p for get_block_ordered_... type methods.
unsigned first_row() const
access function for the first row on this processor. If not distributed then this is just zero...
friend std::ostream & operator<<(std::ostream &o_stream, const BlockSelector &block_selector)
Output function, outputs the Row_index, Column_index, Wanted and the address of the Replacement_block...
bool any_mesh_distributed() const
Check if any of the meshes are distributed. This is equivalent to problem.distributed() and is used a...
bool block_output_on() const
Test if output of blocks is on or not.
unsigned first_row() const
access function for the first row on this processor
void set_column_index(const unsigned &column_index)
Set the column index.
Vector< unsigned > Dof_number_to_block_number_lookup
Vector to the mapping from DOF number to block number.
const unsigned & column_index() const
returns the column index.
void insert_auxiliary_block_distribution(const Vector< unsigned > &block_vec_number, LinearAlgebraDistribution *dist_pt)
insert a Vector<unsigned> and LinearAlgebraDistribution* pair into Auxiliary_block_distribution_pt. The Auxiliary_block_distribution_pt should only contain pointers to distributions concatenated at this block level. We try to ensure this by checking if the block_vec_number vector is within the range nblock_types(). Of course, this does not guarantee correctness, but this is the least we can do.
const bool & wanted() const
returns whether the block is wanted or not.
LinearAlgebraDistribution * Internal_preconditioner_matrix_distribution_pt
The distribution of the (internal) preconditioner matrix. This is formed by concatenating the distrib...
void internal_get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
A helper function, takes the naturally ordered vector, v, and extracts the n-th block vector...
Vector< LinearAlgebraDistribution * > Block_distribution_pt
The distribution for the blocks.
cstr elem_len * i
Definition: cfortran.h:607
LinearAlgebraDistribution * block_distribution_pt(const unsigned b)
Access function to the block distributions (non-const version).
void turn_off_debug_flag()
Toggles off the debug flag.
virtual DoubleMatrixBase * matrix_pt() const
Get function for matrix pointer.
unsigned internal_block_dimension(const unsigned &b) const
Return the number of degrees of freedom in block b. Note that if this preconditioner acts as a subsid...
bool Wanted
Bool to indicate if we require this block.
LinearAlgebraDistribution * Preconditioner_matrix_distribution_pt
The distribution of the preconditioner matrix. This is the concatenation of the block distribution...
bool Debug_flag
Debugging variable. Set true or false via the access functions turn_on_debug_flag(...) turn_off_debug_flag(...)
Vector< unsigned > Allow_multiple_element_type_in_mesh
Vector of unsigned to indicate which meshes contain multiple element types.
unsigned nrow_local() const
access function for the num of local rows on this processor. If no MPI then Nrow is returned...
void do_not_want_block()
Indicate that we do not want the block (set Wanted to false).
void return_concatenated_block_vector(const Vector< unsigned > &block_vec_number, const DoubleVector &b, DoubleVector &v) const
Takes concatenated block ordered vector, b, and copies its entries to the appropriate entries in the ...
BlockSelector()
Default constructor, initialise block index i, j to 0 and bool to false.
bool is_subsidiary_block_preconditioner() const
Return true if this preconditioner is a subsidiary preconditioner.
Vector< unsigned > Ndof_in_block
Number of types of degree of freedom associated with each block.
unsigned ndof_types_in_mesh(const unsigned &i) const
Return the number of DOF types in mesh i. WARNING: This should only be used by the upper-most master ...
Vector< unsigned > Index_in_dof_block_dense
This was uncommented Presumably a non-distributed analogue of Index_in_dof_block_sparse.
Vector< Vector< unsigned > > doftype_coarsen_map_fine() const
Access function for the Doftype_coarsen_map_fine variable.
void turn_on_debug_flag()
Toggles on the debug flag.
Vector< unsigned > Nrows_to_send_for_get_ordered
The number global rows to be sent to processor p for get_block_ordered_... type methods.
MapMatrix< unsigned, CRDoubleMatrix * > Replacement_dof_block_pt
The replacement dof-level blocks.
void return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
Vector< unsigned > Dof_number_dense
Vector to store the mapping from the global DOF number to its block. Empty if this preconditioner has...
Vector< unsigned > Index_in_dof_block_sparse
Vector to store the mapping from the global DOF number to the index (row/column number) within its bl...
OomphInfo oomph_info
BlockPreconditioner< MATRIX > * master_block_preconditioner_pt() const
Access function to the master block preconditioner pt.
void null_replacement_block_pt()
Set Replacement_block_pt to null.
void setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const Vector< unsigned > &block_col_indices)
Setup a matrix vector product. matvec_prod_pt is a pointer to the MatrixVectorProduct, block_pt is a pointer to the block matrix, block_col_indices is a vector indicating which block indices does the RHS vector we want to multiply the matrix by.
Vector< Vector< unsigned > > Doftype_coarsen_map_coarse
Mapping for dof types within THIS precondition. This is usually passed down from the parent precondit...
bool is_master_block_preconditioner() const
Return true if this preconditioner is the master block preconditioner.
Vector< unsigned > Doftype_in_master_preconditioner_fine
The map between the dof types in this preconditioner and the master preconditioner. If there is no master preconditioner it remains empty. This list contains the mapping for the underlying dof types.
void get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
void concatenate(const Vector< LinearAlgebraDistribution * > &in_distribution_pt, LinearAlgebraDistribution &out_distribution)
Takes a vector of LinearAlgebraDistribution objects and concatenates them such that the nrow_local of...
unsigned internal_ndof_types() const
Return the number of internal dof types. This is the number of most fine grain dof types...
void get_block_ordered_preconditioner_vector(const DoubleVector &v, DoubleVector &w)
Given the naturally ordered vector, v, return the vector rearranged in block order in w...
BlockSelector(const unsigned &row_index, const unsigned &column_index, const bool &wanted, CRDoubleMatrix *replacement_block_pt=0)
Constructor, takes the row and column indices and a boolean indicating if the block is required or no...
Vector< Vector< unsigned > > Global_index
Vectors of vectors for the mapping from block number and block row to global row number. Empty if this preconditioner has a master preconditioner as this information is obtain from the master preconditioner.
void operator=(const BlockPreconditioner &)
Broken assignment operator.
static bool Run_block_matrix_test
Static boolean to allow block_matrix_test(...) to be run. Defaults to false.
CRDoubleMatrix * Replacement_block_pt
Pointer to the block.
unsigned long nrow() const
Return the number of rows of the matrix.
Definition: matrices.h:504
LinearAlgebraDistribution * distribution_pt() const
access to the LinearAlgebraDistribution
unsigned Nrow
Number of DOFs (# of rows or columns in the matrix) in this preconditioner. Note that this informatio...
void return_block_ordered_preconditioner_vector(const DoubleVector &w, DoubleVector &v) const
Takes the block ordered vector, w, and reorders it in natural order. Reordered vector is returned in ...
Vector< unsigned > get_fine_grain_dof_types_in(const unsigned &i) const
Returns the most fine grain dof types in a (possibly coarsened) dof type.
const unsigned & row_index() const
returns the row index.
unsigned rank_of_global_row(const unsigned i) const
return the processor rank of the global row number i
void build(const unsigned &row_index, const unsigned &column_index, const bool &wanted, CRDoubleMatrix *replacement_block_pt=0)
void set_replacement_dof_block(const unsigned &block_i, const unsigned &block_j, CRDoubleMatrix *replacement_dof_block_pt)
Set replacement dof-level blocks. Only dof-level blocks can be set. This is important due to how the ...
const unsigned ncol() const
return the number of columns. This is the size of the first inner vectors, or returns 0 if the outer ...
void get_dof_level_block(const unsigned &i, const unsigned &j, MATRIX &output_block, const bool &ignore_replacement_block=false) const
Gets dof-level block (i,j). If Replacement_dof_block_pt(i,j) is not null, then the replacement block ...
void internal_return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
A helper function, takes the vector of block vectors, s, and copies its entries into the naturally or...
Describes the distribution of a distributable linear algebra type object. Typically this is a contain...
void post_block_matrix_assembly_partial_clear()
A helper method to reduce the memory requirements of block preconditioners. Once the methods get_bloc...
void setup_matrix_vector_product(MatrixVectorProduct *matvec_prod_pt, CRDoubleMatrix *block_pt, const unsigned &block_col_index)
Setup matrix vector product. This is simply a wrapper around the other setup_matrix_vector_product fu...
void want_block()
Indicate that we require the block (set Wanted to true).
void turn_on_recursive_debug_flag()
Toggles on the recursive debug flag. The change goes up the block preconditioning hierarchy...
virtual ~BlockSelector()
Default destructor.
void deep_copy(const CRDoubleMatrix *const in_matrix_pt, CRDoubleMatrix &out_matrix)
Create a deep copy of the matrix pointed to by in_matrix_pt.
Definition: matrices.h:3244
void get_concatenated_block_vector(const Vector< unsigned > &block_vec_number, const DoubleVector &v, DoubleVector &b)
Takes the naturally ordered vector and extracts the blocks indicated by the block number (the values)...
bool Recursive_debug_flag
Debugging variable. Set true or false via the access functions turn_on_recursive_debug_flag(...) turn_off_recursive_debug_flag(...) These will turn on/off the debug flag up the hierarchy.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:998
virtual ~BlockPreconditioner()
Destructor.
const Mesh * mesh_pt(const unsigned &i) const
Access to i-th mesh (of the various meshes that contain block preconditionable elements of the same n...
unsigned nmesh() const
Return the number of meshes in Mesh_pt.
virtual void block_setup()
Determine the size of the matrix blocks and setup the lookup schemes relating the global degrees of f...
const LinearAlgebraDistribution * block_distribution_pt(const unsigned &b) const
Access function to the block distributions (const version).
static char t char * s
Definition: cfortran.h:572
unsigned internal_master_dof_number(const unsigned &b) const
Takes the block number within this preconditioner and returns the corresponding block number in the m...
int internal_index_in_block(const unsigned &i_dof) const
Return the index in the block corresponding to a global block number i_dof. The index returned corres...
CRDoubleMatrix * replacement_block_pt() const
Returns Replacement_block_pt.
void get_block_other_matrix(const unsigned &i, const unsigned &j, MATRIX *in_matrix_pt, MATRIX &output_matrix)
Get a block from a different matrix using the blocking scheme that has already been set up...
void concatenate_without_communication(const Vector< LinearAlgebraDistribution * > &row_distribution_pt, const Vector< LinearAlgebraDistribution * > &col_distribution_pt, const DenseMatrix< CRDoubleMatrix * > &matrix_pt, CRDoubleMatrix &result_matrix)
Concatenate CRDoubleMatrix matrices.
Definition: matrices.cc:5164
MapMatrix< unsigned, CRDoubleMatrix * > replacement_dof_block_pt() const
Access function to the replaced dof-level blocks.
unsigned nrow() const
access function to the number of global rows.
unsigned long ncol() const
Return the number of columns of the matrix.
Definition: matrices.h:507
void set_nmesh(const unsigned &n)
Specify the number of meshes required by this block preconditioner. Note: elements in different meshe...
unsigned Internal_nblock_types
Number of different block types in this preconditioner. Note that this information is maintained if u...
const LinearAlgebraDistribution * internal_block_distribution_pt(const unsigned &b) const
Access function to the internal block distributions.
std::string Output_base_filename
String giving the base of the files to write block data into. If empty then do not output blocks...
std::map< Vector< unsigned >, LinearAlgebraDistribution * > Auxiliary_block_distribution_pt
Stores any block-level distributions / concatenation of block-level distributions required...
const unsigned nrow() const
returns the number of rows. This is the outer Vector size.
void setup(CRDoubleMatrix *matrix_pt, const LinearAlgebraDistribution *col_dist_pt=0)
Setup the matrix vector product operator. WARNING: This class is wrapper to Trilinos Epetra matrix ve...
void set_mesh(const unsigned &i, const Mesh *const mesh_pt, const bool &allow_multiple_element_type_in_mesh=false)
Set the i-th mesh for this block preconditioner. Note: The method set_nmesh(...) must be called befor...
Vector< LinearAlgebraDistribution * > Internal_block_distribution_pt
Storage for the default distribution for each internal block.
void set_block_output_to_files(const std::string &basefilename)
Set the base part of the filename to output blocks to. If it is set then all blocks will be output at...
void block_matrix_test(const unsigned &i, const unsigned &j, const MATRIX *block_matrix_pt) const
Private helper function to check that every element in the block matrix (i,j) matches the correspondi...
void get_block(const unsigned &i, const unsigned &j, MATRIX &output_matrix, const bool &ignore_replacement_block=false) const
Put block (i,j) into output_matrix. This block accounts for any coarsening of dof types and any repla...
void select_block(const unsigned &row_index, const unsigned &column_index, const bool &wanted, CRDoubleMatrix *replacement_block_pt=0)
Select a block.
void clear_distribution()
clear the distribution of this distributable linear algebra object
Vector< Vector< unsigned > > Block_to_dof_map_fine
Mapping for the block types to the most fine grain dof types.
virtual const OomphCommunicator * comm_pt() const
Get function for comm pointer.
BlockPreconditioner< MATRIX > * Master_block_preconditioner_pt
If the block preconditioner is acting a subsidiary block preconditioner then a pointer to the master ...
bool built() const
access function to the Built flag - indicates whether the matrix has been build - i...
Definition: matrices.h:1177
unsigned ndof_types() const
Return the total number of DOF types.
unsigned internal_nblock_types() const
Return the number internal blocks. This should be the same as the number of internal dof types...
DenseMatrix< unsigned > Nrows_to_send_for_get_block
The number of global rows to be sent of block b to processor p (matrix indexed [b][p]).
DenseMatrix< int * > Rows_to_recv_for_get_block
The block rows to be received from processor p for block b (matrix indexed [b][p]).
Preconditioner base class. Gives an interface to call all other preconditioners through and stores th...
MATRIX get_block(const unsigned &i, const unsigned &j, const bool &ignore_replacement_block=false) const
Return block (i,j). If the optional argument ignore_replacement_block is true, then any blocks in Rep...
LinearAlgebraDistribution * dof_block_distribution_pt(const unsigned &b)
Access function to the dof-level block distributions.
Vector< unsigned > Dof_number_sparse
Vector to store the mapping from the global DOF number to its block (empty if this preconditioner has...
Vector< int * > Rows_to_recv_for_get_ordered
The preconditioner rows to be received from processor p for get_block_ordered_... type methods...
void broken_assign(const std::string &class_name)
Issue error message and terminate execution.
void document()
debugging method to document the setup. Should only be called after block_setup(...).
unsigned Internal_ndof_types
Number of different DOF types in this preconditioner. Note that this information is maintained if use...
void turn_into_subsidiary_block_preconditioner(BlockPreconditioner< MATRIX > *master_block_prec_pt, const Vector< unsigned > &doftype_in_master_preconditioner_coarse)
Function to turn this preconditioner into a subsidiary preconditioner that operates within a bigger "...
unsigned internal_dof_block_dimension(const unsigned &i) const
Return the size of the dof "block" i, i.e. how many degrees of freedom are associated with it...
int internal_dof_number(const unsigned &i_dof) const
Return the number of the block associated with global unknown i_dof. If this preconditioner is a subs...
std::string string(const unsigned &i)
Return the i-th string or "" if the relevant string hasn't been defined.
void get_blocks(DenseMatrix< bool > &required_blocks, DenseMatrix< MATRIX * > &block_matrix_pt) const
Get all the block matrices required by the block preconditioner. Takes a pointer to a matrix of bools...
void return_block_vectors(const Vector< unsigned > &block_vec_number, const Vector< DoubleVector > &s, DoubleVector &v) const
Takes the vector of block vectors, s, and copies its entries into the naturally ordered vector...
Matrix vector product helper class - primarily a wrapper to Trilinos's Epetra matrix vector product m...
unsigned internal_index_in_dof(const unsigned &i_dof) const
Return the row/column number of global unknown i_dof within it's block.
Vector< Vector< unsigned > > Block_to_dof_map_coarse
Mapping for block types to dof types. These are the dof types the writer of the preconditioner expect...
const LinearAlgebraDistribution * preconditioner_matrix_distribution_pt() const
Access function to the preconditioner matrix distribution pointer. This is the concatenation of the b...
Class for dense matrices, storing all the values of the matrix as a pointer to a pointer with assorte...
Definition: communicator.h:50
Vector< LinearAlgebraDistribution * > Dof_block_distribution_pt
Storage for the default distribution for each dof block at this level.
int internal_block_number(const unsigned &i_dof) const
Return the block number corresponding to a global index i_dof. This returns the block number correspo...
Vector< Vector< unsigned > > Block_number_to_dof_number_lookup
Vector of vectors to store the mapping from block number to the DOF number (each element could be a v...
Vector< Vector< unsigned > > Doftype_coarsen_map_fine
Mapping the dof types within this preconditioner. The values in here refers to the most grain dof typ...
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
DenseMatrix< unsigned > Nrows_to_recv_for_get_block
The number of block rows to be received from processor p for block b (matrix indexed [b][p])...
void internal_return_block_ordered_preconditioner_vector(const DoubleVector &w, DoubleVector &v) const
Takes the block ordered vector, w, and reorders it in the natural order. Reordered vector is returned...
unsigned master_nrow() const
Return the number of dofs (number of rows or columns) in the overall problem. The prefix "master_" is...
Vector< unsigned > Doftype_in_master_preconditioner_coarse
The map between the dof types in this preconditioner and the master preconditioner. If there is no master preconditioner, it remains empty. This is the version for which the master preconditioner expects. The dof types in here may or may not be coarsened in the preconditioner above this one.
const LinearAlgebraDistribution * internal_preconditioner_matrix_distribution_pt() const
access function to the internal preconditioner matrix distribution pt. preconditioner_matrix_distribu...
unsigned Row_index
Row index of the block.
void internal_return_block_vector(const unsigned &n, const DoubleVector &b, DoubleVector &v) const
Takes the n-th block ordered vector, b, and copies its entries to the appropriate entries in the natu...
Vector< unsigned > Dof_dimension
Vector containing the size of each block, i.e. the number of global DOFs associated with it...
MATRIX * matrix_pt() const
Access function to matrix_pt. If this is the master then cast the matrix pointer to MATRIX*...
Data structure to store information about a certain "block" or sub-matrix from the overall matrix in ...
void set_master_matrix_pt(MATRIX *in_matrix_pt)
Set the matrix_pt in the upper-most master preconditioner.
void get_block_vector(const unsigned &n, const DoubleVector &v, DoubleVector &b) const
Takes the naturally ordered vector, v and returns the n-th block vector, b. Here n is the block numbe...
unsigned nblock_types() const
Return the number of block types.
void internal_get_block_ordered_preconditioner_vector(const DoubleVector &v, DoubleVector &w) const
Given the naturally ordered vector, v, return the vector rearranged in block order in w...
BlockPreconditioner(const BlockPreconditioner &)
Broken copy constructor.
unsigned ndof_types() const
Return number of dof types in mesh.
Definition: mesh.cc:8415
A class for compressed row matrices. This is a distributable object.
Definition: matrices.h:872
void internal_get_block_vectors(const Vector< unsigned > &block_vec_number, const DoubleVector &v, Vector< DoubleVector > &s) const
Takes the naturally ordered vector and rearranges it into a vector of sub vectors corresponding to th...
int block_number(const unsigned &i_dof) const
Return the block number corresponding to a global index i_dof.
void disable_block_output_to_files()
Turn off output of blocks (by clearing the basefilename string).
Vector< unsigned > Ndof_types_in_mesh
Storage for number of types of degree of freedom of the elements in each mesh.
A general mesh class.
Definition: mesh.h:74
int get_index_of_value(const Vector< myType > &vec, const myType val, const bool sorted=false) const
Get the index of first occurrence of value in a vector. If the element does not exist, -1 is returned. The optional parameter indicates of the Vector is sorted or not. Complexity: if the Vector is sorted, then on average, logarithmic in the distance between first and last: Performs approximately log2(N)+2 element comparisons. Otherwise, up to linear in the distance between first and last: Compares elements until a match is found.
DenseMatrix< int * > Rows_to_send_for_get_block
The global rows to be sent of block b to processor p (matrix indexed [b][p]).
std::string to_string(T object, unsigned float_precision=8)
Conversion function that should work for anything with operator<< defined (at least all basic types)...
MATRIX get_concatenated_block(const VectorMatrix< BlockSelector > &selected_block)
Returns a concatenation of the block matrices specified by the argument selected_block. The VectorMatrix selected_block must be correctly sized as it is used to determine the number of sub block matrices to concatenate.
const LinearAlgebraDistribution * master_distribution_pt() const
Access function to the distribution of the master preconditioner. If this preconditioner does not hav...
int index_in_block(const unsigned &i_dof) const
Given a global dof number, returns the index in the block it belongs to. This is the overall index...
void set_row_index(const unsigned &row_index)
Set the row index.
void turn_off_recursive_debug_flag()
Toggles off the recursive debug flag. The change goes up the block preconditioning hierarchy...
void internal_get_block(const unsigned &i, const unsigned &j, MATRIX &output_block) const
Gets block (i,j) from the matrix pointed to by Matrix_pt and returns it in output_block. This is associated with the internal blocks. Please use the other get_block(...) function.
void resize(const unsigned long &n)
Definition: matrices.h:511