pattern.hpp 109 KB
Newer Older
prudhomm's avatar
prudhomm committed
1
/*
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
 This file is part of the LifeV library
 Copyright (C) 2001,2002,2003,2004 EPFL, INRIA and Politechnico di Milano

 This library is free software; you can redistribute it and/or
 modify it under the terms of the GNU Lesser General Public
 License as published by the Free Software Foundation; either
 version 2.1 of the License, or (at your option) any later version.

 This library is distributed in the hope that it will be useful,
 but WITHOUT ANY WARRANTY; without even the implied warranty of
 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 Lesser General Public License for more details.

 You should have received a copy of the GNU Lesser General Public
 License along with this library; if not, write to the Free Software
 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
prudhomm's avatar
prudhomm committed
18
*/
19
/*----------------------------------------------------------------------*
prudhomm's avatar
prudhomm committed
20
|
21
22
|
|
prudhomm's avatar
prudhomm committed
23
| #Version  0.1 Experimental   07/7/00. Luca Formaggia & Alessandro Veneziani
24
25
26
27
28
29
|
| #Attempt to define a pattern for vectorial problems 12/10/01. Alain Gauthier
|  Construction of a class VBRPatt which is required for the interface
|  with AZTEC.
|
| #Purposes Defines Patterns for normal and mixed matrices
prudhomm's avatar
prudhomm committed
30
31
|  A pattern defines the graph of a sparse matrix. The pattern class builds
|  that graph starting from a Degree of Freedom (DOF) object.
32
|
prudhomm's avatar
prudhomm committed
33
34
35
36
37
38
39
40
|  The patterns are meant to be used also as way of interrogating the local
|  connectivities of degrees of freedom. In order to avoid ambiguities the
|  degrees of freedom ARE ALWAYS numbered from 1 and they are indicated by
|  using the type name ID (in this way we avoid ambiguity with UInt).
|  On the other hand, we may have a FORTRAN or C-like style for the RAW data
|  storing the pattern which may be CSR of MSR format. In this way we hopefully
|  will be able to better interface external Linear Algebra packages
|  Fortran-style (like SPARSEKIT).
41
42
|
*----------------------------------------------------------------------*/
prudhomm's avatar
prudhomm committed
43
#ifndef PATTERN_OFFSET
44
#define PATTERN_OFFSET 0 // for the Fortran (PATTERN_OFFSET=1) vs C (PATTERN_OFFSET=0) numbering
45
46
47
#endif
#ifndef _PATTERN_HH
#define _PATTERN_HH
prudhomm's avatar
prudhomm committed
48
#include <life/lifecore/life.hpp>
49
50
51
52
53
#ifndef INDEX_T
#define INDEX_T UInt
#endif

#ifndef _VEC_UNKNOWN_HH
prudhomm's avatar
prudhomm committed
54
#include <life/lifearray/vecUnknown.hpp>
55
56
57
58
59
60
#endif



#include<set>
#include<algorithm>
prudhomm's avatar
prudhomm committed
61
#include<string>
62
//#include<functional>
prudhomm's avatar
prudhomm committed
63
#include <life/lifemesh/bareItems.hpp>
64

prudhomm's avatar
prudhomm committed
65
66
namespace LifeV
{
67
const INDEX_T PatternOffset = PATTERN_OFFSET;
prudhomm's avatar
prudhomm committed
68

prudhomm's avatar
prudhomm committed
69
70
/*!
  \class PatternDefs
71
72
73
74
75
  This class containes some useful functions and typedefs which will be
  common to ALL patterns (base and mixed ones). I use BareEdges for the
  dynamic pattern, because I have provided the comparison function. Actually,
  just a pair<UInt,UInt> should work all the same.
*/
76
77
78
class PatternDefs
{
public:
prudhomm's avatar
prudhomm committed
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
    //! Type for indices
    typedef INDEX_T Index_t;

    //! Container for the actual (raw) pattern
    typedef std::vector<Index_t> Container;

    //! type for differences (offsets)
    typedef Container::size_type Diff_t;

    //! convert from Identifier (DOF) to Index
    Index_t _d2i( ID const d ) const;

    //! convert from Index to Identifier
    ID _i2d( Index_t const i ) const;

    //! convert from index to differences (offsets)
    Diff_t _i2o( Index_t const i ) const;

    //! convert from identifier to differences
    Diff_t _d2o( ID const d ) const;

    //! convert container of indices to identifiers (in place)
    Container & _i2d( Container & list_of_indices ) const;

    //! convert container of identifiers to indices (in place)
104
    Container & _d2i( Container & list_of_dof ) const;
prudhomm's avatar
prudhomm committed
105
106

    //! convert container of identifiers to differences (in place)
107
108
    Container & _d2o( Container & list_of_dof ) const;
protected:
prudhomm's avatar
prudhomm committed
109
110
    //! Container for the dynamic pattern
    typedef std::set<BareEdge, cmpBareItem<BareEdge> > DynPattern;
111
112
113
};


prudhomm's avatar
prudhomm committed
114
115
/*!
  \class BasePattern
116
117
  This is the building block for the pattern classes.
*/
118
class BasePattern : public PatternDefs
119
120
{
public:
prudhomm's avatar
prudhomm committed
121

122
123
124

  /*!
      \typedef enum
125

126
127
128
129
130
131
132
133
134
135
136
    */
    typedef enum PatternType
    {
        STANDARD_PATTERN,          //!< Standard pattern
        EDGE_COUPLING_PATTERN      //!< Additional coupling trhough edges
    };





prudhomm's avatar
prudhomm committed
137
    //! Default constructor (size zero)
138
    BasePattern();
139

prudhomm's avatar
prudhomm committed
140
141
142
143
144
145
    /*!
      Constructor for given pattern size
      @param ex_nnz number of nonzero entries
      @param ex_nrow number of rows
      @param ex_ncol number of columns
    */
146

prudhomm's avatar
prudhomm committed
147
148
    BasePattern( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol );

149
    //! Default distructor
150
    virtual ~BasePattern() {};
151

prudhomm's avatar
prudhomm committed
152
153
154
155
156
157
158
159
160
161
162
    //! Copy constructor
    BasePattern( BasePattern const &rightHandPattern );

    //! Number of Rows
    inline UInt nRows() const;

    //! Number of Columns
    inline UInt nCols() const;

    //! Total number of nonzero entries
    inline UInt nNz() const;
163

prudhomm's avatar
prudhomm committed
164
165
    //! returns true iff pattern is empty
    inline bool isEmpty() const;
166

prudhomm's avatar
prudhomm committed
167
168
    //! returns true iff the first item furnished by the raw pattern data
    //! is the diagonal entry
169
170
171
172
    bool diagFirst() const
    {
        return _diagfirst;
    }
173

prudhomm's avatar
prudhomm committed
174
175
176
    //! some info for the curious
    void showMe( bool const verbose = false,
                 std::ostream & out = std::cout ) const;
177

178
179
180
    //! Returns true if the indices (i,j) are in the pattern
    virtual bool isThere( Index_t const i, Index_t const j ) const = 0;

181
protected:
prudhomm's avatar
prudhomm committed
182
183
184
185
186
187
    /*!
      Builds the dynamic pattern. It is the standard routine for a
      SYMMETRIC pattern associated to a single degree of freedom object.
      It uses the local pattern provided through the DOF object.
      If nbcomp > 1, it reproduces the pattern for (nbcomp X nbcomp) blocks
    */
188
    template <typename DOF1>
prudhomm's avatar
prudhomm committed
189
190
    bool setpatt( DOF1 const & dof1,
                  DynPattern & dynpatt,
191
192
                  UInt const nbcomp = 1 );

prudhomm's avatar
prudhomm committed
193
194
195
196
197
    /*!
      Builds the dynamic pattern. It uses the local pattern provided through
      the DOF object, augmented by the couplings needed for IP stabilization.
      If nbcomp > 1, it reproduces the pattern for (nbcomp X nbcomp) blocks.
      @author Miguel Fernandez, 12/2003
198
199
200
201
202
     */
    template <typename DOF, typename MESH>
    bool setpatt( const DOF& dof, const MESH& mesh, DynPattern & dynpatt,
                  UInt const nbcomp );

prudhomm's avatar
prudhomm committed
203
204
205
206
207
208
209
210
211
212
213
214
215
    /*!
      Builds the SYMMETRIC pattern associated to StiffDG operators.
      This pattern should work with several penalization terms (IP,
      Bassi et al., ecc...). Boundary integrals involve one element
      and its adjacent neighbours.
      @author Daniele A. di Pietro, 01/2004
    */
    template<typename DOF, typename DOFBYFACE>
    bool setpattDG(DOF const & dof, DOFBYFACE const & dofbyface,
                   DynPattern & dynpatt, UInt const nbcomp = 1);

    /*!
      It builds a mixed pattern. It is the standard routine for a
prudhomm's avatar
prudhomm committed
216
      pattern associated to two degrees of     freedom object.
prudhomm's avatar
prudhomm committed
217
218
219
      it uses the MixedLocalPattern class for the local patterns.
      See the documentation in the implementation part for more details
    */
220
221
    template <typename DOF1, typename DOF2>
    bool setpatt( DOF1 const & dof1, DOF2 const & dof2, DynPattern & dynpatt,
222
                  UInt const bRows, UInt const bCols );
223
224
225
226
227
228
229
230

    static const UInt _defaultsize = 0; // MUST BE 0

    UInt _nnz;
    UInt _nrows;
    UInt _ncols;
    bool _filled; //!< True if the pattern has been filled!
    bool _diagfirst; //!< True if the pattern has diagfirst!
231
232
233
};

///////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
234
235
236
237
//
//      C S R Pattern
//
///////////////////////////////////////////////////
238

239
240
class MSRPatt;

prudhomm's avatar
prudhomm committed
241
class CSRPatt : public BasePattern
242
243
244
{
public:

prudhomm's avatar
prudhomm committed
245
    //! Default constructor
246
    CSRPatt();
prudhomm's avatar
prudhomm committed
247
248
249
250
251
252
253

    /*!
      Constructor
      @param ex_nnz number of nonzero entries
      @param ex_nrow number of rows
      @param ex_ncol number of columns
    */
254
255
    CSRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol );

prudhomm's avatar
prudhomm committed
256
    //! Constructor where raw data are given from outside
257
258
259
    CSRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol,
             const std::vector<Index_t> &ex_ia,
             const std::vector<Index_t> &ex_ja );
prudhomm's avatar
prudhomm committed
260
261
262
263
264
265
266
267
268
269
270
271
272
273
    /*
      TODO: As it
      stands, This is a pretty useless contructor, unless we build a smarter
      version which takes in input raw data also in a other types. This
      indeed will help interfacing with external libraries We will need to
      change it to something of the type

      template<typename Const_Iter>
      CSRPatt(Const_Iter &ex_ia, Const_Iter  &ex_ja )
      where Const_iter is an iterator to a sequence (then  a pointer as well)
      LUCA: I will do it soon!
    */

    //! Copy constructor
274
    CSRPatt( const CSRPatt &RightHandCSRP );
prudhomm's avatar
prudhomm committed
275
276

    //! Constructor for single DOF (square matrix), possibly with more than
277
278
279
    //! one component
    template <typename DOF>
    CSRPatt( DOF const & dof, UInt const nbcomp = 1 );
prudhomm's avatar
prudhomm committed
280
281
282

    //! build function for single DOF (square matrix), possibly with more than
    //! one component
283
284
285
    template <typename DOF>
    bool buildPattern( DOF const & dof, const UInt nbcomp );

prudhomm's avatar
prudhomm committed
286
287
288
289
290
    /*!
      Constructor for single DOF (square matrix), possibly with more than
      one component; version handling patterns coming from IP stabilization.
      @author Christoph Winkelmann, 09/2004
    */
291
292
    template <typename DOF, typename MESH>
    CSRPatt( const DOF& dof, const UInt nbcomp, const MESH& mesh );
prudhomm's avatar
prudhomm committed
293
294
295
296
297

    /* build function for single dof(square matrix), possibly with more than
       one component; version handling patterns coming from IP stabilization.
       @author Christoph Winkelmann, 09/2004
    */
298
299
300
    template <typename DOF, typename MESH>
    bool buildPattern( DOF const& dof, const UInt nbcomp, const MESH& mesh );

prudhomm's avatar
prudhomm committed
301
302
303
304
305
306
    /*!
      Constructors for two DOF (in general non-square matrix)
      This constructors are required if we want to use the pattern in a
      MixedPattern. It may be useful also stand alone. Look at the
      documentation of buildPattern<DOF1,DOF2> to have details.
    */
307
308
309
310
311
312
313
314
    template <typename DOF1, typename DOF2>
    CSRPatt( DOF1 const & dof1, DOF2 const & dof2,
             UInt const bRows = 1, UInt const bCols = 1 );

    template <typename DOF1, typename DOF2>
    bool buildPattern( DOF1 const& dof1, DOF2 const & dof2,
                       UInt const bRows = 1, UInt const bCols = 1 );

prudhomm's avatar
prudhomm committed
315
316
317
318
319
320
321
322
    CSRPatt( const MSRPatt& msrPatt );

    /*!
      Version which construct two patterns: patt and its transpose one,
      in addition an other Container for the access to
      the columns of the transpose matrix is built.
      @author Alain, Nov. 2002
    */
323
    template <typename DOF1, typename DOF2>
324
    friend void
325
326
327
328
329
330
    buildPattTpatt( DOF1 const& dof1, DOF2 const & dof2,
                    CSRPatt &Patt, CSRPatt &Tpatt,
                    UInt const bRows = 1, UInt const bCols = 1 );

    CSRPatt & operator= ( const CSRPatt& RhCsr );

prudhomm's avatar
prudhomm committed
331
332
333
334
335
    //! return a copy of ia as a container
    const Container& ia() const { return _ia; };

    //! return a copy of ja as a container
    const Container& ja() const { return _ja; };
336

prudhomm's avatar
prudhomm committed
337
338
    //! return a copy of jaT as a container
    const Container& jaT() const { return _jaT; };
339

prudhomm's avatar
prudhomm committed
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
    //! Give ia (in a raw form)
    Index_t * giveRawCSR_ia() { return & ( _ia.front() ); }

    //! Give ja (in a raw form)
    Index_t * giveRawCSR_ja() { return & ( _ja.front() ); }

    //! Give jaT (in a raw form)
    Index_t * giveRawCSR_jaT() { return & ( _jaT.front() ); }

    //! Give ia (in a raw form)
    Index_t const * giveRawCSR_ia() const { return & ( _ia.front() ); }

    //! Give ja (in a raw form)
    Index_t const * giveRawCSR_ja() const { return & ( _ja.front() ); }

    //! Give jaT (in a raw form)
    Index_t const * giveRawCSR_jaT() const { return & ( _jaT.front() ); }

    //! Give ia (as container)
    Container & give_ia() { return _ia; }

    //! Give ja (as container)
    Container & give_ja() { return _ja; }

    //! Give jaT (as container)
    Container & give_jaT() { return _jaT; }

    //! Give ia (as container)
    Container const & give_ia() const { return _ia; }

    //! Give ja (as container)
    Container const & give_ja() const { return _ja; }

    //! Give jaT (as container)
    Container const & give_jaT() const { return _jaT; }

    //! Min and Max elements in a row.
    inline std::pair<UInt, UInt> giveMinMax() const;

    //! N of neighbours of the DOF numbeered d . BEWARE d is INCLUDED!
    inline UInt nbNeighbours( ID const d ) const;

    //! the i-th (start from 1) neighbour of dof d. The first is d itself
    inline ID neighbour( ID const i, ID const d ) const;

    //! put neighbours of dof d in a list
    inline void neighbours( ID const d, Container & start ) const;

    /*!
      The following template function extracts a row (useful to implement A*b),
      returning two sequences: One with the column numbering (coldata) of one
      with the corresponding offsets in the vector holding the matrix value.
      coldata is given as a sequence of Index_t (i.e. the content depend by
      the value of PATETRN_OFFSET), while the position sequence are always
      offsets (i.e. starting from 0) Iter is either an iterator to a container
      or a pointer. The function Returns the n. or row elements.

      IMPORTANT: for efficency reason the sequences pointed by coldata and
      position MUST have been dimensioned boforehand in order to have the
      sufficient dimension (use giveMinMax) NO CHECKS ARE MADE!
    */
401
402
403
    template <typename Iter>
    inline UInt row( Diff_t const row, Iter coldata, Iter position ) const;

prudhomm's avatar
prudhomm committed
404
405
406
407
408
409
410
411
412
413
414
    //! Locate position in vector containing DOF ID couple (i,j)
    inline std::pair<Diff_t, bool> locate_dof( ID const i,
                                               ID const j ) const;

    //! Locate position in vector containing index couple (i,j)
    inline std::pair<Diff_t, bool> locate_index( Index_t const i,
                                                 Index_t const j ) const;

    //! pattern visualization
    void showMe( bool const verbose = false,
                 std::ostream& c = std::cout ) const ;
415

prudhomm's avatar
prudhomm committed
416
417
    //! pattern visualization a la Matlab
    void spy( std::string const & filename = "matrice.m" ) const;
418
419
420

    //! column-concatenation of two CSR block patterns
    friend CSRPatt colUnify( CSRPatt const &patt1, CSRPatt const &patt2 );
prudhomm's avatar
prudhomm committed
421
422

    //! column-concatenation of one CSR block pattern and ncolZero null columns
423
    friend CSRPatt colUnify( CSRPatt const &patt1, UInt const ncolZero );
prudhomm's avatar
prudhomm committed
424
425

    //! column-concatenation of one ncolZero null columns and CSR block pattern
426
    friend CSRPatt colUnify( UInt const ncolZero, CSRPatt const &patt1 );
prudhomm's avatar
prudhomm committed
427

428
429
    //! row-concatenation of two CSR block patterns
    friend CSRPatt rowUnify( CSRPatt const &patt1, CSRPatt const &patt2 );
prudhomm's avatar
prudhomm committed
430
431

    //! row-concatenation of one CSR block pattern and nrowZero null rows
432
    friend CSRPatt rowUnify( CSRPatt const &patt1, UInt const nrowZero );
prudhomm's avatar
prudhomm committed
433
434

    //! row-concatenation of nrowZero null rows and one CSR block pattern
435
436
437
438
    friend CSRPatt rowUnify( UInt const nrowZero, CSRPatt const &patt1 );

    //! Return a block diagonal pattern of nblock blocks
    friend CSRPatt diagblockMatrix( CSRPatt const &patt, UInt const nblock );
439

440
441
    //! Returns true if the indices (i,j) are in the pattern
    bool isThere( Index_t const i, Index_t const j ) const;
442
443

protected:
prudhomm's avatar
prudhomm committed
444
445
446
447
448
449
450
451
    //! Row offset for row index i
    Diff_t _row_off( Index_t i ) const { return _i2o( _ia[ _i2o( i ) ] ); }


    //! Locate position in vector containing index couple (i,j)
    std::pair<Diff_t, bool> locate_pattern( Index_t const i,
                                            Index_t const j ) const;

452
453
454
    Container _ia; //!< point to the rows entries
    Container _ja; //!< contain col indices for each row
    Container _jaT; //!< point to the col entries of the transpose pattern
455
456

private:
457
    template <typename DOF>
prudhomm's avatar
prudhomm committed
458
459
    void _buildPattern( DOF const & dof,
                        DynPattern const & dynpatt,
460
                        UInt const nbcomp );
461
462
463
};

///////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
464
465
466
467
//
//      V B R Pattern
// It holds a CSR pattern of variable size blocks
//
468
469
//////////////////////////////////////////////////

prudhomm's avatar
prudhomm committed
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
/*!
  \class VBRPatt
  The block pattern is given by the class CSRPatt:
  _ia : points to the location in _ja of the first block
  entry in each block row.
  _ja : contains the block column indices of the block pattern.

  BEWARE : the interpretation of _nnz, _nrows and _ncols of the class
  BasePattern is valid for the BLOCK pattern, that is
  _nnz : number of non zero BLOCKS,
  _nrows : number of BLOCK rows,
  _ncols : number of BLOCK columns.

  VBRPatt: VBR format accepts blocks of variable size.
  This pattern is more general than needed. In fact we work
  with blocks having a fixed size.
*/
487
488
class VBRPatt: public CSRPatt
{
489
490
491
492
public:
    VBRPatt()
    {}
    VBRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol );
prudhomm's avatar
prudhomm committed
493
494
495
496
497
498
    VBRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol,
             const std::vector<Index_t> &ex_ia,
             const std::vector<Index_t> &ex_ja,
             const std::vector<Index_t> &ex_indx,
             const std::vector<Index_t> &ex_rpntr,
             const std::vector<Index_t> &ex_cpntr );
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515

    VBRPatt( const VBRPatt &RightHandVBRP );

    //! Constructors for single DOF (square matrix), the size of the
    //! blocks is given by GenericVecHdl.
    template <typename DOF>
    VBRPatt( DOF const & dof, UInt const blockSize );
    bool buildPattern( UInt const blockSize );

    //! Constructors for two DOF (in general non-square matrix)
    //! as for a single DOF, the block pattern is built by CSRPatt and
    //! the size of the blocks is given by GenericVecHdl.
    template <typename DOF1, typename DOF2>
    VBRPatt( DOF1 const & dof1, DOF2 const & dof2, UInt const blockSize );

    VBRPatt & operator= ( const VBRPatt& RhVbr );

prudhomm's avatar
prudhomm committed
516
517
    //! return a copy of indx as a container
    inline const Container& indx() const { return _indx; }
518

prudhomm's avatar
prudhomm committed
519
520
    //! return a copy of rpntr as a container
    inline const Container& rpntr() const { return _rpntr; }
521

prudhomm's avatar
prudhomm committed
522
523
    //! return a copy of cpntr as a container
    inline const Container& cpntr() const { return _cpntr; }
524

prudhomm's avatar
prudhomm committed
525
526
527
528
529
    //! Give ia (in a raw form)
    inline Index_t* giveRawVBR_ia() { return giveRawCSR_ia(); }

    //! Give ja (in a raw form)
    inline Index_t* giveRawVBR_ja() { return giveRawCSR_ja(); }
530

prudhomm's avatar
prudhomm committed
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
    //! Give indx (in a raw form)
    inline Index_t* giveRawVBR_indx() { return & ( _indx.front() ); }

    //! Give rpntr (in a raw form)
    inline Index_t* giveRawVBR_rpntr() { return & ( _rpntr.front() ); }

    //! Give cpntr (in a raw form)
    inline Index_t* giveRawVBR_cpntr() { return & ( _cpntr.front() ); }

    //! Give indx (as container)
    inline Container & give_indx() { return _indx; }

    //! Give rpntr (as container)
    inline Container & give_rpntr() { return _rpntr; }

    //! Give cpntr (as container)
    inline Container & give_cpntr() { return _cpntr; }

    //! Give the block row to which row belongs. No check of _rpntr size!
550
    inline Diff_t rbloc( Diff_t const row ) const
551
    {
prudhomm's avatar
prudhomm committed
552
553
554
555
556
        // size of square block
        UInt blsize = _rpntr[ 1 ] - _rpntr[ 0 ];

        // block row offset in which is row
        return row / blsize + PatternOffset;
557
558
    }

prudhomm's avatar
prudhomm committed
559
    //! Give the block column to which col belongs
560
    inline Diff_t cbloc( Diff_t const col ) const
561
    {
prudhomm's avatar
prudhomm committed
562
563
564
565
566
        // size of square block
        UInt blsize = _rpntr[ 1 ] - _rpntr[ 0 ];

        // block row offset in which is row
        return col / blsize + PatternOffset;
567
568
    }

prudhomm's avatar
prudhomm committed
569
570
    //! Give the local row numbering
    //! of row in the corresponding block (given by rbloc)
571
    inline Diff_t locr( Diff_t const row ) const
572
    {
prudhomm's avatar
prudhomm committed
573
574
575
576
577
        // size of square block
        UInt blsize = _rpntr[ 1 ] - _rpntr[ 0 ];

        // local row number into the block
        return row % blsize + PatternOffset;
578
579
    }

prudhomm's avatar
prudhomm committed
580
581
    //! Give the local column numbering
    //! of row in the corresponding block (given by cbloc)
582
583
    inline Diff_t locc( Diff_t const col ) const
    {
prudhomm's avatar
prudhomm committed
584
585
586
587
588
        // size of square block
        UInt blsize = _rpntr[ 1 ] - _rpntr[ 0 ];

        // local col number into the block
        return col % blsize + PatternOffset;
589
590
    }

prudhomm's avatar
prudhomm committed
591
592
593
    /*!
      The following template function extracts a row (useful to implement A*b),
      returning two sequences:
594

prudhomm's avatar
prudhomm committed
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
      coldata : column numbering (of elements and NOT blocks)

      position: coldata corresponding offsets in the vector holding the matrix
      value.

      coldata is given as a sequence of Index_t (i.e. the
      content depends on the value of PATTERN_OFFSET), while
      the position sequence are always offsets (i.e. starting from 0)

      Iter is either an iterator to a container or a pointer. The
      function returns the n. or row elements.

      IMPORTANT: for efficiency reason the sequences pointed by coldata
      and position MUST have been dimensioned before in order to
      have the sufficient dimension (use giveMinMax) NO CHECKS ARE MADE!
    */
    inline UInt row( Diff_t const row,
                     Container & coldata,
                     Container & position ) const;
614
615
616
617
618
619
620
621

    // Here the routines which locate the position in the vector
    // containing the matrix value of an entry (i,j), which may be a
    // couple of indices or of ID's (degree of freedoms identifiers)
    // ALAIN : using locate routines member of class CSRPatt gives
    // information on the block position (element (1,1) of the block).

    //! pattern visualization
prudhomm's avatar
prudhomm committed
622
623
624
625
626
    void showMe( bool const verbose = false,
                 std::ostream& c = std::cout ) const ;

    //! pattern visualization a la Matlab
    void spy( std::string const & filename = "matrice.m" ) const;
627
628

private:
prudhomm's avatar
prudhomm committed
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
    //! _indx(i) points to the location in the pattern
    //! of the (1,1) element of the i-th block entry.
    //! _indx(_nnz+1)=1+number of non zeros elements of the pattern.
    Container _indx;

    //! _rpntr(i)-_rpntr(0) is the row index of the
    //! first point row in the i-th block row.
    //! _rpntr(_nrows+1)=_rpntr(0)+number of element rows.
    Container _rpntr;

    //!  the column index of the first point column in
    //! the i-th block column.
    //! _cpntr=_rpntr FOR SQUARE PATTERN of SQUARE BLOCKS.
    //! _cpntr(_ncols+1)=_cpntr(0)+number of element columns.
    Container _cpntr;
644
645
646
};

//////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
647
648
649
650
//
//      C S R Symmetric Pattern
// It holds only the upper triangular part
//
651
652
//////////////////////////////////////////////////

prudhomm's avatar
prudhomm committed
653
class CSRPattSymm : public BasePattern
654
{
655
656
public:
    CSRPattSymm();
prudhomm's avatar
prudhomm committed
657

658
    CSRPattSymm( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol );
prudhomm's avatar
prudhomm committed
659

660
661
662
    CSRPattSymm( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol,
                 const std::vector<Index_t> &ex_ia,
                 const std::vector<Index_t> &ex_ja );
prudhomm's avatar
prudhomm committed
663

664
    CSRPattSymm( const CSRPattSymm &RightHandCSRP );
prudhomm's avatar
prudhomm committed
665

666
667
    template <typename DOF>
    CSRPattSymm( DOF const & dof );
prudhomm's avatar
prudhomm committed
668

669
670
671
    CSRPattSymm & operator= ( const CSRPattSymm& RhCsr );

    template <typename DOF>
prudhomm's avatar
prudhomm committed
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
    bool buildPattern( DOF const & dof );

    const Container& ia() const { return _ia; };

    const Container& ja() const { return _ja; };

    //! Give ia (in a raw form)
    Index_t * giveRawCSR_ia() { return & ( _ia.front() ); }

    //! Give ja (in a raw form)
    Index_t * giveRawCSR_ja() { return & ( _ja.front() ); }

    //! Give ia (as container)
    Container & give_ia() { return _ia; }

    //! Give ja (as container)
    Container & give_ja() { return _ja; }
689

690
    inline std::pair<UInt, UInt> giveMinMax() const;
691

prudhomm's avatar
prudhomm committed
692
693
    //! N of neighbours of the DOF numbeered d . BEWARE d is INCLUDED!
    //! Note: this function is inefficient
694
    UInt nbNeighbours( ID const d ) const;
prudhomm's avatar
prudhomm committed
695
696
697

    //! the i-th (start from 1) neighbour of dof d. The first is d itself
    //! Note: this function is inefficient
698
    ID neighbour( ID const i, ID const d ) const ;
prudhomm's avatar
prudhomm committed
699
700
701

    //! put neighbours of dof d in a list
    //! Note: this function is inefficient
702
    void neighbours( ID const d, Container & start ) const;
prudhomm's avatar
prudhomm committed
703
704

    //! extracts a row (useful to implement A*b)
705
    template <typename Iter>
prudhomm's avatar
prudhomm committed
706
    inline UInt row( Diff_t const row, Iter coldata, Iter position ) const;
707

708
    inline std::pair<Diff_t, bool> locate_dof( ID const i, ID const j ) const;
709

prudhomm's avatar
prudhomm committed
710
711
712
713
714
    inline std::pair<Diff_t, bool> locate_index( Index_t const i,
                                                 Index_t const j ) const;

    //! pattern visualization
    void showMe( bool verbose = false, std::ostream& c = std::cout ) const;
715

prudhomm's avatar
prudhomm committed
716
717
    //! pattern visualization a la Matlab
    void spy( std::string const & filename = "matrice.m" ) const;
prudhomm's avatar
prudhomm committed
718
719


prudhomm's avatar
prudhomm committed
720
721
722
723
724
    /*
      DOF are ALWAYS numbered from 1: the following definition set the correct
      numbering as a function of current numbering (specified by
      PATTERN_OFFSET)
    */
725

726
    //! Returns true if the indices (i,j) are in the pattern
727
728
    bool isThere( Index_t i, Index_t j ) const;

729
protected:
prudhomm's avatar
prudhomm committed
730
731
732
733
734
735
    //! Row offset for row index i
    Diff_t _row_off( Index_t i ) const { return _i2o( _ia[ _i2o( i ) ] ); }

    std::pair<Diff_t, bool> locate_pattern( Index_t const i,
                                            Index_t const j ) const;

736
737
738
private:
    Container _ia;
    Container _ja;
739
740
741
742
};


///////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
743
744
745
746
//
//      M S R Pattern
//  It implements the MSR format
//
747
748
//////////////////////////////////////////////////

prudhomm's avatar
prudhomm committed
749
class MSRPatt : public BasePattern
750
751
{
public:
prudhomm's avatar
prudhomm committed
752

753
    MSRPatt();
prudhomm's avatar
prudhomm committed
754

755
    MSRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol );
prudhomm's avatar
prudhomm committed
756

757
758
759
    MSRPatt( UInt ex_nnz, UInt ex_nrow, UInt ex_ncol,
             const std::vector<Index_t> &bindx,
             const std::vector<Index_t> &ybind );
prudhomm's avatar
prudhomm committed
760

761
    MSRPatt( const MSRPatt &RightHandMSRP );
prudhomm's avatar
prudhomm committed
762

763
    MSRPatt( const CSRPatt &RightHandCSRP );
prudhomm's avatar
prudhomm committed
764

765
    MSRPatt & operator= ( const MSRPatt& RhMsr );
prudhomm's avatar
prudhomm committed
766
767

    //! Constructor for standard FEM
768
769
770
    template <typename DOF>
    MSRPatt( DOF const & dof, UInt const nbcomp = 1 );

prudhomm's avatar
prudhomm committed
771
772
    //! Constructor for continuous FEM with IP stabilization
    //! @author Miguel Fernandez, 12/2003
773
    template <typename DOF, typename MESH>
774
    MSRPatt( const DOF& dof, const MESH& mesh, const UInt nbcomp, PatternType patternType = EDGE_COUPLING_PATTERN );
775

prudhomm's avatar
prudhomm committed
776
777
778
779
780
781
782
783
784
    //! Constructor for DG FEM
    //! @author Daniele A. Di Pietro, 10/2004
    template<typename DOF, typename DOFBYFACE>
    MSRPatt(DOF const & dof,
            DOFBYFACE const & dofbyface,
            const std::string& type,
            UInt const nbcomp = 1);

    //! build function for standard FEM
785
786
787
    template <typename DOF>
    bool buildPattern( DOF const & dof, UInt const nbcomp );

prudhomm's avatar
prudhomm committed
788
789
    //! build function for continuous FEM with IP stabilization
    //! @author Miguel Fernandez, 12/2003
790
791
792
    template <typename DOF, typename MESH>
    bool buildPattern( const DOF& dof, const MESH& mesh, const UInt nbcomp );

prudhomm's avatar
prudhomm committed
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
    //! build function for DG FEM
    //! @author Daniele A. Di Pietro, 10/2004
    template<typename DOF, typename DOFBYFACE>
    bool buildPattern(DOF const & dof,
                      DOFBYFACE const & dofbyface,
                      const std::string& type,
                      UInt const nbcomp = 1);

    const Container& bindx() const { return _bindx; };

    const Container& ybind() const { return _ybind; };

    //! Give _bindx (in a raw form)
    Index_t * giveRaw_bindx() { return & ( _bindx.front() ); }

    //! Give _ybind (in a raw form)
    Index_t * giveRaw_ybind() { return & ( _ybind.front() ); }

    //! Give _bindx (in a raw form)
    Index_t const * giveRaw_bindx() const { return & ( _bindx.front() ); }

    //! Give _ybind (in a raw form)
    Index_t const * giveRaw_ybind() const { return & ( _ybind.front() ); }

    //! Give _bindx (as container)
    Container & give_bindx() { return _bindx; }

    //! Give _ybind (as container)
    Container & give_ybind() { return _ybind; }

    //! Give _bindx (as container)
    Container const & give_bindx() const { return _bindx; }

    //! Give _ybind (as container)
    Container const & give_ybind() const { return _ybind; }

829
    inline UInt nbNeighbours( ID const d ) const;
prudhomm's avatar
prudhomm committed
830

831
    inline ID neighbour( ID const i, ID const d ) const;
prudhomm's avatar
prudhomm committed
832

833
    void neighbours( ID const d, Container & start ) const;
prudhomm's avatar
prudhomm committed
834
835

    //! extracts a row (useful to implement A*b)
836
    template <typename Iter>
prudhomm's avatar
prudhomm committed
837
    inline UInt row( Diff_t const row, Iter coldata, Iter position ) const;
838

839
    inline std::pair<Diff_t, bool> locate_dof( ID const i, ID const j ) const;
prudhomm's avatar
prudhomm committed
840

841
842
    inline std::pair<UInt, UInt> giveMinMax() const;

prudhomm's avatar
prudhomm committed
843
844
    inline std::pair<Diff_t, bool> locate_index( Index_t const i,
                                                 Index_t const j ) const;
845

prudhomm's avatar
prudhomm committed
846
847
848
849
850
851
852
853
854
855
    //! pattern visualization
    void showMe( bool verbose = false, std::ostream& c = std::cout ) const ;

    //! pattern visualization a la Matlab
    void spy( std::string const & filename = "matrice.m" ) const ;

    //! Construction of a diagonal matrix of n blocks
    friend void diagblockMatrix( MSRPatt &ans,
                                 MSRPatt const &patt,
                                 UInt const nblock );
856

857
    //! Returns true if the indices (i,j) are in the pattern
858
859
    bool isThere( Index_t i, Index_t j ) const;

860
protected:
prudhomm's avatar
prudhomm committed
861
862
863
864
865
866
    //! Row offset for row index i
    Diff_t _row_off( Index_t i ) const { return _i2o( _bindx[ _i2o( i ) ] ); }

    std::pair<Diff_t, bool> locate_pattern( Index_t const i,
                                            Index_t const j ) const;

867
868
private:
    Container _bindx;
prudhomm's avatar
prudhomm committed
869
870
    Container _ybind;
    // WARNING: AV January 2001: ybind helps in reading the matrix by columns;
871
872
873
874
875
876
877
    // a counterpart of ybind should be added also to the CSR format
    // ybind should facilitate also the locate_pattern subroutines
    // IT IS BASED ON THE ASSUMPTION THAT THE PATTERN IS SYMMETRIC

    template <typename DOF>
    void _buildPattern( DOF const & dof, DynPattern const & dynpatt,
                        UInt const nbcomp );
878
879
880
};

///////////////////////////////////////////////////
881
//
prudhomm's avatar
prudhomm committed
882
883
884
885
886
887
888
889
890
891
892
893
//      Mixed Local Patterns
//
//////////////////////////////////////////////////

/*!
  \class MixedLocalPattern
  This class MixedLocalPattern<FE1, FE2> is used to provide the local
  pattern for two basis finite element, it corresponds to the pattern
  required to build element matrices of the type
  \f$ m_{ij}=\int_K \phi_i^1\phi^2_j \f$.
  The number of rows correspond to the number of local
  degrees of freedom for the finite element of type 1, i.e FE1::nbNode,
894
  while the number of colums are FE2::nbNode. The major assumption is that
prudhomm's avatar
prudhomm committed
895
896
897
  \f$ \forall K \f$ and \f$ i,j\in K \f$, the support of
  \f$ \phi_i^1\vert_K \f$ and \f$ \phi_j^2\vert_K \f$ has non-null
  intersection. In other words, the local matrix has a full pattern.
prudhomm's avatar
prudhomm committed
898
  The pattern follows the rule of numbering
prudhomm's avatar
prudhomm committed
899
*/
900
template <typename FE1, typename FE2>
901
902
903
class MixedLocalPattern: public PatternDefs
{
public:
904
905
906
    UInt nRows() const;
    UInt nCols() const;
    UInt nbPattern() const;
prudhomm's avatar
prudhomm committed
907
908
    UInt patternFirst( UInt const i ) const; // Numbering from 0
    UInt patternSecond( UInt const i ) const; // Numbering from 0
909
};
prudhomm's avatar
prudhomm committed
910

911
///////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
912
913
914
//
//      Mixed Pattern
//
915
//////////////////////////////////////////////////
prudhomm's avatar
prudhomm committed
916
917
918
919
920
921
922
923
924

/*!
  \class MixedPattern
  The class MixedPattern is a very general class able to held multiblock
  matrix patterns. The block numbering starts from (0,0). Each block contains
  a pointer to a Pattern class and a couple of offsets, which indicate how the
  LOCAL block row/cols numbering has to be increased to get the GLOBAL
  numbering (i.e. the numbering associated to to the global matrix).

925
  It may be used in two forms
prudhomm's avatar
prudhomm committed
926
927
928
929
930
931
932
933
934
935

  1) as a viewer, then the local patterns are contructued externally and then
  "linked" to the mixed patter object, or

  2) by delegating to contruction of  the local patterns to the mixed pattern
  object.

  In the first case the destruction of the mixed pattern object will NOT imply
  the destruction of the local patterns, while in the other case everything
  is destroyed.
936
*/
937
template <UInt BROWS, UInt BCOLS, typename PATTERN = CSRPatt>
prudhomm's avatar
prudhomm committed
938
class MixedPattern : public PatternDefs
939
940
{
public:
prudhomm's avatar
prudhomm committed
941

prudhomm's avatar
prudhomm committed
942
    //! default constructor
943
944
    MixedPattern();

prudhomm's avatar
prudhomm committed
945
946
947
948
949
950
    /*!
      constructor which makes the pattern construct and link to an external
      pattern, useful to construct a diagonal pattern
      @arg type "full" or "diag"
      @author Miguel Fernandez, 11/2002
    */
951
952
    MixedPattern( PATTERN & ex_patt, const std::string& type = "full" );

prudhomm's avatar
prudhomm committed
953
    //! destructor
954
955
    ~MixedPattern();

prudhomm's avatar
prudhomm committed
956
    //! make a diagonal pattern for vectorial problem
957
958
    void makeDiagPattern( PATTERN & ex_patt );

prudhomm's avatar
prudhomm committed
959
960
961
962
963
    //! Number of blocks (rows and columns)
    inline std::pair<UInt, UInt> nBlocks() const;

    //! Number of rows in block (m,n)
    inline UInt nRows( Diff_t const m, Diff_t const n ) const;
964

965
    inline UInt nCols( Diff_t const m, Diff_t const n ) const;
966

prudhomm's avatar
prudhomm committed
967
968
969
970
971
972
973
974
975
976
977
978
979
    //! Nonzeros on block (m,n)
    inline UInt nNz( Diff_t m, Diff_t n ) const;

    //! Global number of rows
    inline UInt nRows() const;

    //! Global number of cols
    inline UInt nCols() const;

    //! Non zeros in global matrix
    UInt nNz() const;

    //! Number of neighbours at block level. (local ID numbering)
980
    UInt nbNeighbours( Diff_t const m, Diff_t const n, ID const d ) const ;
981

prudhomm's avatar
prudhomm committed
982
983
984
    //! i-th neighbour of dof d at block level. (local ID numbering)
    ID neighbour( Diff_t const m, Diff_t const n,
                  ID const i, ID const d ) const;
985

prudhomm's avatar
prudhomm committed
986
987
988
    //! put neighbours of dof d at block level. (local ID numbering) in a list
    inline void neighbours( Diff_t const m, Diff_t const n,
                            ID const d, Container & neighs ) const;
989

prudhomm's avatar
prudhomm committed
990
991
992
993
994
    //! extracts a row (useful to implement A*b) at block level.
    //! (local ID numbering)
    template <typename Iter>
    inline UInt row( Diff_t const m, Diff_t const n,
                     Diff_t const row, Iter coldata, Iter Position ) const;
995

prudhomm's avatar
prudhomm committed
996
997
    //! Number of neighbours at global level. (global ID numbering)
    UInt nbNeighbours( ID const d_g ) const ;
998

prudhomm's avatar
prudhomm committed
999
1000
    //! i-th neighbour of dof d at global level. (global ID numbering)
    ID neighbour( ID const i_g, ID const d_g ) const;
For faster browsing, not all history is shown. View entire blame