ParallelDofMapping.h 21.5 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ==  http://www.amdis-fem.org                                              ==
// ==                                                                        ==
// ============================================================================
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.



/** \file FeSpaceMapping.h */

23
24
25
26
#ifndef AMDIS_FE_SPACE_MAPPING_H
#define AMDIS_FE_SPACE_MAPPING_H

#include <mpi.h>
27
#include <vector>
28
#include <map>
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
29
#include <set>
30
31
#include <petsc.h>
#include <petscis.h>
32
33
#include <boost/container/flat_map.hpp>
#include <boost/container/flat_set.hpp>
34
35

#include "AMDiS_fwd.h"
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
36
#include "parallel/DofComm.h"
37
38
#include "parallel/MpiHelper.h"
#include "parallel/ParallelTypes.h"
39
#include "parallel/StdMpi.h"
40
41

namespace AMDiS {
Thomas Witkowski's avatar
Thomas Witkowski committed
42

43
44
  using namespace std;

45
  /** \brief
46
47
48
   * Is used to store matrix indices to all DOFs in rank's subdomain. Thus, the
   * class defines a mapping from component number and DOF index to a global
   * matrix index. This class does not calculate the indices by itself!
49
50
51
52
53
54
55
56
57
58
59
60
   */
  class DofToMatIndex
  {
  public:
    DofToMatIndex() {}

    /// Reset the data structure.
    inline void clear()
    {
      data.clear();
    }

61
62
63
64
65
66
67
    /** Add a new mapping for a given DOF.
     * 
     * \param[in]   component       Component number for which the mapping 
     *                              is defined.
     * \param[in]   dof             DOF index
     * \param[in]   globalMatIndex  Global matrix index.
     */
68
69
70
71
72
    inline void add(int component, DegreeOfFreedom dof, int globalMatIndex)
    {
      data[component][dof] = globalMatIndex;
    }

73
    /// Maps a global DOF index to the global matrix index for a specific 
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
    /// system component number.
    inline int get(int component, DegreeOfFreedom dof)
    {
      FUNCNAME("DofToMatIndex::get()");

      TEST_EXIT_DBG(data.count(component))
	("No mapping data for component %d available!\n", component);

      TEST_EXIT_DBG(data[component].count(dof))
	("Mapping for DOF %d in component %d does not exists!\n",
	 dof, component);

      return data[component][dof];
    }

89
90
91
92
93
    /// Returns for a given matrix index the component and (local or global) DOF
    /// index. As the data structure is not made for this kind of reverse
    /// search, this is very slow and should be only used for debugging.
    void getReverse(int rowIndex, int &component, int &dofIndex);

94
95
96
  private:
    /// The mapping data. For each system component there is a specific map that
    /// maps global DOF indices to global matrix indices.
97
    map<int, boost::container::flat_map<DegreeOfFreedom, int> > data;
98
99
100
  };


101
102
103
  /// This class defines the parallel mapping of DOFs for one FE space. It is
  /// used by the class \ref ParallelDofMapping to specifiy the mapping for a
  /// set of  FE spaces.
104
  class ComponentDofMap
105
106
  {
  public:
107
108
    /// This constructor exists only to create std::map of this class and make
    /// use of the operator [] for read access. Should never be called.
109
    ComponentDofMap() 
110
111
112
    {
      ERROR_EXIT("Should not be called!\n");
    }
113
114
     
    /// This is the only valid constructur to be used. 
115
    ComponentDofMap(MeshLevelData* ld);
116
117

    /// Clears all data of the mapping.
118
    void clear();
119

120
121
    /// Maps a DOF index to both, the local and global index of the mapping. The
    /// global index must not be set.
Thomas Witkowski's avatar
Thomas Witkowski committed
122
    MultiIndex& operator[](DegreeOfFreedom d)
123
    {
124
      TEST_EXIT_DBG(dofMap.count(d))("DOF %d is not in map!\n", d);
125

126
      return dofMap[d];
127
    }
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146

    /** \brief
     * Searchs the map for a given DOF. It does not fail, if the DOF is not 
     * mapped by this mapping. In this case, it returns false. If the DOF is 
     * mapped, the result is stored and the function returns true.
     *
     * \param[in]    dof     DOF index for which a mapping is searched.
     * \param[out]   index   In the case that the DOF is mapped, the result
     *                       is stored here.
     */
    inline bool find(DegreeOfFreedom dof, MultiIndex& index)
    {
      DofMap::iterator it = dofMap.find(dof);
      if (it == dofMap.end())
	return false;

      index = it->second;
      return true;
    }
147
    
148
149
    /// Inserts a new DOF to rank's mapping. The DOF is assumed to be owend by
    /// the rank.
150
    void insertRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
151
    {
152
      FUNCNAME("ComponentDofMap::insertRankDof()");
153
      
154
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
155
      
156
157
      dofMap[dof0].local = dof1;
      nLocalDofs++;
Thomas Witkowski's avatar
Thomas Witkowski committed
158
      if (dof1 != -1)
159
	nRankDofs++;
160
161
    }
    
162
163
    /// Inserts a new DOF to rank's mapping. The DOF exists in rank's subdomain
    /// but is owned by a different rank, thus it is part of an interior boundary.
164
    void insertNonRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
165
    {
166
      FUNCNAME("ComponentDofMap::insertNonRankDof()");
167
      
168
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
169
      
170
171
172
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      nonRankDofs.insert(dof0);
173
174
    }
    
175
    /// Checks if a given DOF is in the DOF mapping.
176
    bool isSet(DegreeOfFreedom dof)
177
    {
178
      return static_cast<bool>(dofMap.count(dof));
179
    }
180

Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
181
182
183
    /// Checks if a given DOF is a rank owned DOF of the DOF mapping. The DOF
    /// must be a DOF of the mapping (this is not checked here), otherwise the
    /// result is meaningsless.
184
    bool isRankDof(DegreeOfFreedom dof)
185
    {
186
      return !(static_cast<bool>(nonRankDofs.count(dof)));
187
    }
Thomas Witkowski's avatar
Blbu  
Thomas Witkowski committed
188
189
190
191
192

    bool isRankGlobalDof(int dof)
    {
      return (dof >= rStartDofs && dof < rStartDofs + nRankDofs);
    }
193
    
194
    /// Returns number of DOFs in the mapping.
195
    unsigned int size()
196
    {
197
      return dofMap.size();
198
199
    }
    
200
    /// Returns the raw data of the mapping.
201
    DofMap& getMap()
202
    {
203
      return dofMap;
204
    }
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
205

206
207
208
209
210
211
212
213
214
215
    DofMap::iterator begin()
    {
      return dofMap.begin();
    }

    DofMap::iterator end()
    {
      return dofMap.end();
    }

216
217
    /// Recomputes the mapping.
    void update();
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
218

219
    /// Sets the FE space this mapping corresponds to.
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
220
221
222
223
224
    void setFeSpace(const FiniteElemSpace *fe)
    {
      feSpace = fe;
    }

225
226
    /// Informs the mapping whether the mapping will include DOFs that are not
    /// owned by the rank.
227
    void setNonLocal(bool b)
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
228
    {
229
      isNonLocal = b;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
230
231
    }

232
    /// Informs the mapping whether a global index must be computed.
233
    void setNeedGlobalMapping(bool b)
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
234
    {
235
      needGlobalMapping = b;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
236
237
    }

238
    /// Sets the DOF communicator.
239
    void setDofComm(DofComm &dc)
Thomas Witkowski's avatar
Thomas Witkowski committed
240
    {
241
      dofComm = &dc;
Thomas Witkowski's avatar
Thomas Witkowski committed
242
243
    }

244
245
246
247
248
249
    void setMpiComm(MPI::Intracomm &m, int l)
    {
      mpiComm = m;
      meshLevel = l;
    }

250
  private:
251
    /// Computes a global mapping from the local one.
252
    void computeGlobalMapping();
253
254
255

    /// Computes the global indices of all DOFs in the mapping that are not owned
    /// by the rank.
256
    void computeNonLocalIndices();
257
258

  private:
259
    MeshLevelData *levelData;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
260

261
262
    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;
263

264
265
266
267
    MPI::Intracomm mpiComm;

    int meshLevel;

268
    /// The FE space this mapping belongs to. This is used only the get the
269
    /// correct DOF communicator in \ref dofComm.
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
270
271
    const FiniteElemSpace *feSpace;

272
    /// Mapping data from DOF indices to local and global indices.
273
    DofMap dofMap;
274

275
    /// Set of all DOFs that are in mapping but are not owned by the rank.
276
    boost::container::flat_set<DegreeOfFreedom> nonRankDofs;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
277

278
    /// If true, a global index mapping will be computed for all DOFs.
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    bool needGlobalMapping;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
280

281
282
283
    /// Is true if there are DOFs in at least one subdomain that are not owned
    /// by the rank. If the value is false, each rank contains only DOFs that
    /// are also owned by this rank.
284
    bool isNonLocal;
285

286
287
  public:
    /// 
288
    int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;
289
  };
290

291

292
293
294
295
  class ComponentIterator {
  public:
    virtual ComponentDofMap& operator*() = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
    virtual ComponentDofMap* operator->() = 0;

298
299
300
301
302
303
304
    virtual bool end() = 0;
    
    virtual void next() = 0;

    virtual void reset() = 0;
  };

305

306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
  class ComponentDataInterface
  {
  public:
    virtual ComponentDofMap& operator[](int compNumber) = 0;

    virtual ComponentDofMap& operator[](const FiniteElemSpace *feSpace) = 0;

    virtual bool isDefinedFor(int compNumber) const = 0;

    virtual ComponentIterator& getIteratorData() = 0;
    
    virtual ComponentIterator& getIteratorComponent() = 0;

    virtual void init(vector<const FiniteElemSpace*> &f0,
		      vector<const FiniteElemSpace*> &f1,
321
322
		      bool isNonLocal,
		      MeshLevelData &levelData) = 0;
323

Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
327
328
    vector<const FiniteElemSpace*>& getFeSpaces()
    {
      return feSpaces;
    }

329
330
331
332
333
334
335
336
337
338
  protected:
    /// The FE spaces for all components.
    vector<const FiniteElemSpace*> feSpaces;

    /// The set of all FE spaces. It uniquly contains all different FE spaces
    /// from \ref feSpaces.
    vector<const FiniteElemSpace*> feSpacesUnique;
  };


339
  class ComponentDataEqFeSpace : public ComponentDataInterface
340
341
  {
  public:
342
343
344
345
346
    ComponentDataEqFeSpace()
      : iterData(this),
	iterComponent(this)
    {}

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
    ComponentDofMap& operator[](int compNumber)
    {
      const FiniteElemSpace *feSpace = feSpaces[compNumber];
      return componentData.find(feSpace)->second;
    }

    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
    {
      TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");;
      
      return componentData.find(feSpace)->second;
    }

    bool isDefinedFor(int compNumber) const
    {
      const FiniteElemSpace *feSpace = feSpaces[compNumber];
      return static_cast<bool>(componentData.count(feSpace));
    }

    ComponentIterator& getIteratorData()
    {
      iterData.reset();
      return iterData;
    }
    
    ComponentIterator& getIteratorComponent()
    {
      iterComponent.reset();
      return iterComponent;
    }

    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
380
381
	      bool isNonLocal,
	      MeshLevelData &levelData);
382
383


384
  protected:
385

386
387
    void addFeSpace(const FiniteElemSpace* feSpace,
		    MeshLevelData &levelData);
388
389
390

    class IteratorData : public ComponentIterator {
    public:
391
392
393
394
      IteratorData(ComponentDataEqFeSpace *d)
	: data(d)
      {}

395
396
      ComponentDofMap& operator*()
      {
397
	(*data)[*it];
398
399
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
      ComponentDofMap* operator->()
      {
402
	&((*data)[*it]);
Thomas Witkowski's avatar
Thomas Witkowski committed
403
404
      }

405
406
      bool end()
      {
407
	return (it == data->feSpacesUnique.end());
408
409
410
411
      }

      void next()
      {
412
	++it;
413
414
415
416
      }

      void reset()
      {
417
	it = data->feSpacesUnique.begin();
418
      }
419
420
421
422
423

    protected:
      ComponentDataEqFeSpace *data;

      vector<const FiniteElemSpace*>::iterator it;
424
425
426
427
    };

    class IteratorComponent : public ComponentIterator {
    public:
428
429
430
431
      IteratorComponent(ComponentDataEqFeSpace *d)
	: data(d)
      {}

432
433
      ComponentDofMap& operator*()
      {
434
	(*data)[*it];
435
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438

      ComponentDofMap* operator->()
      {
439
	&((*data)[*it]);
Thomas Witkowski's avatar
Thomas Witkowski committed
440
441
      }

442
443
      bool end()
      {
444
	return (it == data->feSpaces.end());
445
446
447
448
      }

      void next()
      {
449
	++it;
450
451
452
453
      }

      void reset()
      {
454
	it = data->feSpaces.begin();
455
      }
456
457
458
459
460

    protected:
      ComponentDataEqFeSpace *data;

      vector<const FiniteElemSpace*>::iterator it;
461
462
463
464
465
466
467
468
    };
    

    map<const FiniteElemSpace*, ComponentDofMap> componentData;

    IteratorData iterData;

    IteratorComponent iterComponent;
469
470
471
472

    friend class IteratorData;
    
    friend class IteratorComponent;
473
474
475
  };


476
  class ComponentDataDiffFeSpace : public ComponentDataInterface
477
478
  {
  public:
479
480
481
482
    ComponentDataDiffFeSpace()
      : iter(this)
    {}

483
484
485
486
487
488
489
490
491
    ComponentDofMap& operator[](int compNumber)
    {
      TEST_EXIT_DBG(componentData.count(compNumber))("No data for component %d!\n", compNumber);

      return componentData.find(compNumber)->second;
    }

    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
    {
492
      ERROR_EXIT("FE Space acces is not possible for component wise defined DOF mappings\n");
493
494
495
496
    }
    
    ComponentIterator& getIteratorData()
    {
497
498
      iter.reset();
      return iter;
499
500
501
502
    }
    
    ComponentIterator& getIteratorComponent()
    {
503
504
      iter.reset();
      return iter;
505
506
507
508
509
510
511
    }

    bool isDefinedFor(int compNumber) const
    {
      return (static_cast<unsigned int>(compNumber) < componentData.size());
    }

512
513
514
515
516
    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
	      bool isNonLocal,
	      MeshLevelData &levelData);

517
  protected:
518
519
520
521
522
    void addComponent(unsigned int component,
		      const FiniteElemSpace* feSpace,
		      MeshLevelData &levelData);

    class Iterator : public ComponentIterator {
523
    public:
524
525
526
527
528
      Iterator(ComponentDataDiffFeSpace *d)
	: data(d),
	  componentCounter(-1)
      {}

529
530
      ComponentDofMap& operator*()
      {
531
	(*data)[componentCounter];
532
533
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
534
535
      ComponentDofMap* operator->()
      {
536
	&((*data)[componentCounter]);
Thomas Witkowski's avatar
Thomas Witkowski committed
537
538
539
      }


540
541
      bool end()
      {
542
	return (it == data->feSpaces.end());
543
544
545
546
      }

      void next()
      {
547
548
	++it;
	++componentCounter;
549

550
551
	if (it == data->feSpaces.end())
	  componentCounter = -1;
552
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
553

554
      void reset()
Thomas Witkowski's avatar
Thomas Witkowski committed
555
      {
556
557
	it = data->feSpaces.begin();
	componentCounter = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
558
559
      }

560
561
    protected:
      ComponentDataDiffFeSpace *data;
562

563
      vector<const FiniteElemSpace*>::iterator it;
564

565
      int componentCounter;
566
567
568
569
    };

    map<unsigned int, ComponentDofMap> componentData;

570
    Iterator iter;
571

572
    friend class Iterator;
573
574
  };

575
576
577
578
579
580
  /// Used to specify whether a parallel DOF mapping is defined for each
  /// specific component or for each FE space.
  enum DofMappingMode {
    COMPONENT_WISE,
    FESPACE_WISE
  };
581

582
583
584
585
586
  /**
   * Implements the mapping from sets of distributed DOF indices to local and
   * global indices. The mapping works for a given set of FE spaces. Furthermore,
   * this class may compute the matrix indices of the set of DOF indices.
   */
Thomas Witkowski's avatar
Thomas Witkowski committed
587
  class ParallelDofMapping
588
589
  {
  public:
590
    ParallelDofMapping(DofMappingMode mode);
591

592
593
594
595
    /** \brief Initialize the parallel DOF mapping.
     *
     * \param[in]  m                  MPI communicator.
     * \param[in]  fe                 The FE spaces of all components of the 
596
597
598
599
     *                                PDE to be solved.
     * \param[in]  uniqueFe           Unique list of FE spaces. Thus, two
     *                                arbitrary elements of this list are always
     *                                different.
600
     * \param[in]  isNonLocal         If true, at least one rank's mapping con-
601
602
     *                                taines DOFs that are not owend by the rank.
     */
603
    void init(MeshLevelData& mld,
604
	      vector<const FiniteElemSpace*> &fe,
605
	      vector<const FiniteElemSpace*> &uniqueFe,
606
607
608
609
610
611
612
613
614
615
616
617
	      bool isNonLocal = true);

    /// In the case of having only one FE space, this init function can be used.
    void init(MeshLevelData& mld,
	      const FiniteElemSpace *feSpace,
	      bool isNonLocal = true)
    {
      vector<const FiniteElemSpace*> feSpaces;
      feSpaces.push_back(feSpace);
      init(mld, feSpaces, feSpaces, isNonLocal);
    }

618
    void setMpiComm(MPI::Intracomm &m, int l);
619
620
621
    
    /// Clear all data.
    void clear();
622
623
624

    /// Set the DOF communicator objects that are required to exchange information
    /// about DOFs that are on interior boundaries.
625
    void setDofComm(DofComm &dofComm);
626

627
628
629
    /// Returns the DOF communicator.
    DofComm& getDofComm()
    {
630
631
632
      FUNCNAME("ParallelDofMapping::getDofComm()");

      TEST_EXIT_DBG(dofComm)("No DOF communicator object defined!\n");
633
634
635
636

      return *dofComm;
    }

637
638
639
    /// Changes the computation of matrix indices based of either local or
    /// global DOF indices, see \ref needMatIndexFromGlobal
    void setComputeMatIndex(bool global)
640
641
642
643
    {
      needMatIndexFromGlobal = global;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
644
645
646
647
648
    inline bool isMatIndexFromGlobal()
    {
      return needMatIndexFromGlobal;
    }

649

650
651
652
653
    /// Access the DOF mapping for a given component number.
    inline ComponentDofMap& operator[](int compNumber)
    {
      return (*data)[compNumber];
654
655
    }

656
657
    /// Access the DOF mapping for a given FE space    
    inline ComponentDofMap& operator[](const FiniteElemSpace *feSpace) 
Thomas Witkowski's avatar
Thomas Witkowski committed
658
    {
659
      return (*data)[feSpace];
Thomas Witkowski's avatar
Thomas Witkowski committed
660
661
    }

662
    inline bool isDefinedFor(int compNumber) const
Thomas Witkowski's avatar
Thomas Witkowski committed
663
    {
664
      return data->isDefinedFor(compNumber);
Thomas Witkowski's avatar
Thomas Witkowski committed
665
666
    }

667
    /// Returns the number of solution components the mapping is defined on.
Thomas Witkowski's avatar
Thomas Witkowski committed
668
669
670
671
    inline int getNumberOfComponents() const
    {      
      return static_cast<int>(data->getFeSpaces().size());
    }
672

673
    /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank.
674
    inline int getRankDofs()
675
    {
676
      TEST_EXIT_DBG(nRankDofs >= 0)("Should not happen!\n");
677

678
      return nRankDofs;
679
680
    }

681
    /// Returns \ref nLocalDofs, thus the number of DOFs in ranks subdomain.
682
    inline int getLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
683
    {
684
      TEST_EXIT_DBG(nLocalDofs >= 0)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
685

686
      return nLocalDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
687
688
    }

689
    /// Returns \ref nOverallDofs, thus the number of all DOFs in this mapping.
690
    inline int getOverallDofs()
691
    {
692
      TEST_EXIT_DBG(nOverallDofs >= 0)("Should not happen!\n");
693

694
      return nOverallDofs;
695
696
    }

697
698
    /// Returns \ref rStartDofs, thus the smallest global index of a DOF that is
    /// owned by the rank.
699
    inline int getStartDofs()
700
    {
701
      TEST_EXIT_DBG(rStartDofs >= 0)("Should not happen!\n");
702

703
      return rStartDofs;
704
705
    }

706
    /// Update the mapping.
Thomas Witkowski's avatar
Thomas Witkowski committed
707
    void update();
708

709
710
    /// Returns the global matrix index of a given DOF for a given 
    /// component number.
711
    inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
Thomas Witkowski's avatar
Thomas Witkowski committed
712
    {
713
      return dofToMatIndex.get(ithComponent, d);
Thomas Witkowski's avatar
Thomas Witkowski committed
714
715
    }

716
717
718
719
720
721
722
    /// Returns the component number and local/global DOF index for a given
    /// matrix row index. Should be used for debugging only!
    inline void getReverseMatIndex(int index, int &component, int &dofIndex)
    {
      dofToMatIndex.getReverse(index, component, dofIndex);
    }

723
724
    /// Returns the local matrix index of a given DOF for a given 
    /// component number.
725
    inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d)
726
    {
727
      return dofToMatIndex.get(ithComponent, d) - rStartDofs;
728
729
    }

730
731
732
733
734
    // inline const FiniteElemSpace* getFeSpace(int i = 0)
    // {      
    //   TEST_EXIT_DBG(i < feSpacesUnique.size())("Wrong component number!\n");
    //   return feSpacesUnique[i];
    // }
Thomas Witkowski's avatar
Thomas Witkowski committed
735

Thomas Witkowski's avatar
Thomas Witkowski committed
736
737
738
739
740
741
742
743
744
745
746
747
    // Writes all data of this object to an output stream.
    void serialize(ostream &out)
    {
      ERROR_EXIT("MUST BE IMPLEMENTED!\n");
    }

    // Reads the object data from an input stream.
    void deserialize(istream &in)
    {
      ERROR_EXIT("MUST BE IMPLEMENTED!\n");
    }

748
    /// Compute local and global matrix indices.
749
    void computeMatIndex(bool globalIndex);
750

751
752
753
754
    void createIndexSet(IS &is, 
			int firstComponent, 
			int nComponents);

755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
    /// Create a parallel distributed PETSc vector based on this mapping.
    inline void createVec(Vec &vec)
    {
      VecCreateMPI(mpiComm, getRankDofs(), getOverallDofs(), &vec);
    }

    /// Create a parallel distributed PETsc vector based on this mapping but
    /// with a different (larger) global size. This is used in multi-level 
    /// method to embed a local vector into a subdomain spaned by several
    /// ranks.
    inline void createVec(Vec &vec, int nGlobalRows)
    {
      VecCreateMPI(mpiComm, getRankDofs(), nGlobalRows, &vec);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
770
771
772
773
774
    inline void createLocalVec(Vec &vec)
    {
      VecCreateSeq(PETSC_COMM_SELF, getRankDofs(), &vec);
    }

775
776
  protected:
    /// Compute \ref nRankDofs.
777
    int computeRankDofs();
778
779

    /// Compute \ref nLocalDofs.
780
    int computeLocalDofs();
781
782

    /// Compute \ref nOverallDofs.
783
    int computeOverallDofs();
784
785

    /// Compute \ref rStartDofs.
786
    int computeStartDofs();
787

788
  private:
789
790
791
792
    MPI::Intracomm mpiComm;

    int meshLevel;

793
    MeshLevelData *levelData;
794

795
796
    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;
797

798
799
800
    /// Is true if there are DOFs in at least one subdomain that are not owned
    /// by the rank. If the value is false, each rank contains only DOFs that
    /// are also owned by this rank.
801
    bool isNonLocal;
802
803
804
805
806

    /// If matrix indices should be computed, this variable defines if the
    /// mapping from DOF indices to matrix row indices is defined on local
    /// or global DOF indices. If true, the mapping is to specify and to use
    /// on global ones, otherwise on local DOF indices.
Thomas Witkowski's avatar
Thomas Witkowski committed
807
808
809
810
811
    /// In most scenarios the mapping stored on local DOF indices is what we
    /// want to have. Only when periodic boundary conditions are used together
    /// with some global matrix approache, the matrix indices must be stored
    /// also for DOFs that are not part of the local subdomain. Thus, the
    /// mapping will be stored on global DOF indices.
812
813
    bool needMatIndexFromGlobal;

814
815
    /// Maps from components to DOF mappings.
    ComponentDataInterface *data;
Thomas Witkowski's avatar
Thomas Witkowski committed
816

817
    /// Number of DOFs owned by rank.
818
    int nRankDofs;
819
820

    /// Number of DOFs in rank's subdomain.
821
    int nLocalDofs;
822
823

    /// Number of global DOFs (this value is thus the same on all ranks).
824
    int nOverallDofs;
825
826

    /// Smallest global index of a DOF owned by the rank.
827
    int rStartDofs;
828
829
830

    /// Mapping from global DOF indices to global matrix indices under 
    /// consideration of possibly multiple components.
831
    DofToMatIndex dofToMatIndex;
832
833
834
835
  };
}

#endif