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