ParallelDofMapping.h 23.8 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
    /// Informs the mapping whether a global index must be computed.
226
    void setGlobalMapping(bool b)
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
227
    {
228
      globalMapping = b;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
229
230
    }

231
    /// Sets the DOF communicator.
232
    void setDofComm(DofComm &dc)
Thomas Witkowski's avatar
Thomas Witkowski committed
233
    {
234
      dofComm = &dc;
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
    }

237
238
239
240
241
242
    void setMpiComm(MPI::Intracomm &m, int l)
    {
      mpiComm = m;
      meshLevel = l;
    }

243
  private:
244
    /// Computes a global mapping from the local one.
245
    void computeGlobalMapping();
246
247
248

    /// Computes the global indices of all DOFs in the mapping that are not owned
    /// by the rank.
249
    void computeNonLocalIndices();
250
251

  private:
252
    MeshLevelData *levelData;
Thomas Witkowski's avatar
FETI-DP  
Thomas Witkowski committed
253

254
255
    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;
256

257
258
259
260
    MPI::Intracomm mpiComm;

    int meshLevel;

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

265
    /// Mapping data from DOF indices to local and global indices.
266
    DofMap dofMap;
267

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

271
    /// If true, a global index mapping will be computed for all DOFs.
272
    bool globalMapping;
273

274
275
  public:
    /// 
276
    int nRankDofs, nLocalDofs, nOverallDofs, rStartDofs;
277
  };
278

279

280
281
  /// Abstract iterator interface to access containrs containing values of
  /// type \ref ComponentDofMap.
282
283
284
285
  class ComponentIterator {
  public:
    virtual ComponentDofMap& operator*() = 0;

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

288
289
290
291
292
293
294
    virtual bool end() = 0;
    
    virtual void next() = 0;

    virtual void reset() = 0;
  };

295

296
297
298
  /// Abstract interface to acces DOF mapping data for each component. Allows
  /// to hide specific implementations, which allow, e.g., to efficiently map
  /// all components having the same FE space to the same DOF mapping.
299
300
301
  class ComponentDataInterface
  {
  public:
302
    /// Access via component number
303
304
    virtual ComponentDofMap& operator[](int compNumber) = 0;

305
    /// Acess via FE space pointer
306
307
    virtual ComponentDofMap& operator[](const FiniteElemSpace *feSpace) = 0;

308
309
    /// Checks whether the DOF mapping is defined for a specific 
    /// component number.
310
311
    virtual bool isDefinedFor(int compNumber) const = 0;

312
    /// Returns iterator which iterates over all DOF mappings.
313
314
    virtual ComponentIterator& getIteratorData() = 0;
    
315
316
317
318
    /// Returns iterator which iterates over the DOF mappings of all 
    /// components. If the data is defined for each FE space and more than
    /// one commponent is defined on the same FE space, the iterative will
    /// also iterative multple times over the same DOF mapping object.
319
320
    virtual ComponentIterator& getIteratorComponent() = 0;

321
322
323
324
325
326
327
328
329
330
331
332
    /** \brief
     * Initialization of the object.
     *
     * \param[in]   componentSpaces   Set of the FE spaces for each component.
     * \param[in]   feSpaces          Set of all different FE spaces.
     * \param[in]   globalMapping     Mapping is parallel (true) or only 
     *                                local (false).
     * \param[in]   levelData         Data for multi level method.
     */
    virtual void init(vector<const FiniteElemSpace*> &componentSpaces,
		      vector<const FiniteElemSpace*> &feSpaces,
		      bool globalMapping,
333
		      MeshLevelData &levelData) = 0;
334
335
336

  protected:
    /// The FE spaces for all components.
337
    vector<const FiniteElemSpace*> componentSpaces;
338
339
340

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


345
346
347
  /// This class concretizes the interface class \ref ComponentDataInterface. A
  /// DOF mapping is implemented for each component.
  class ComponentData : public ComponentDataInterface
348
349
  {
  public:
350
351
    ComponentData()
      : iter(this)
352
353
    {}

354
    /// Returns DOF mapping for a given component number.
355
356
    ComponentDofMap& operator[](int compNumber)
    {
357
358
      TEST_EXIT_DBG(componentData.count(compNumber))
	("No data for component %d!\n", compNumber);
359

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

363
364
365
366
367
    /// Just to implement the corresponding virtual function in \ref 
    /// ComponentDataInterface. Of course it does not work as we have data for
    /// each component. Thus there may be different mappings for the same
    /// FE space.
    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
368
    {
369
      ERROR_EXIT("FE Space acces is not possible for component wise defined DOF mappings\n");
370
    }
371
372
    
    /// Return data iterator.
373
374
    ComponentIterator& getIteratorData()
    {
375
376
      iter.reset();
      return iter;
377
378
    }
    
379
    /// Return component iterator.
380
381
    ComponentIterator& getIteratorComponent()
    {
382
383
      iter.reset();
      return iter;
384
385
    }

386
387
388
389
390
391
392
393
    /// Checks whether the DOF mapping is defined for a specific 
    /// component number.
    bool isDefinedFor(int compNumber) const
    {
      return (static_cast<unsigned int>(compNumber) < componentData.size());
    }

    /// Initialization
394
395
    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
396
	      bool globalMapping,
397
	      MeshLevelData &levelData);
398

399
  protected:
400

401
402
    /// Iterator class to iterate over all parallel DOF mappings.
    class Iterator : public ComponentIterator {
403
    public:
404
405
406
      Iterator(ComponentData *d)
	: data(d),
	  componentCounter(-1)
407
408
      {}

409
410
      ComponentDofMap& operator*()
      {
411
	(*data)[componentCounter];
412
413
      }

Thomas Witkowski's avatar
Thomas Witkowski committed
414
415
      ComponentDofMap* operator->()
      {
416
	&((*data)[componentCounter]);
Thomas Witkowski's avatar
Thomas Witkowski committed
417
418
      }

419

420
421
      bool end()
      {
422
	return (it == data->componentSpaces.end());
423
424
425
426
      }

      void next()
      {
427
	++it;
428
429
430
431
	++componentCounter;

	if (it == data->componentSpaces.end())
	  componentCounter = -1;
432
433
434
435
      }

      void reset()
      {
436
437
	it = data->componentSpaces.begin();
	componentCounter = 0;
438
      }
439
440

    protected:
441
442
      /// Pointer to data class over which the iterator must iterate.
      ComponentData *data;
443

444
      /// Internal iterator of the internal data from \ref ComponentData.
445
      vector<const FiniteElemSpace*>::iterator it;
446

447
448
449
      /// Component number of current iteration. 
      int componentCounter;
    };
450
451


452
453
    /// Data mapping from component numbers to DOF mapping objects.
    map<unsigned int, ComponentDofMap> componentData;
454

455
456
    /// Iterator object.
    Iterator iter;
457
458


459

460
    friend class Iterator;
461
462
463
  };


464
465
466
467
468

  /// This class concretizes the interface class \ref ComponentDataInterface. A
  /// DOF mapping is implemented for each finite element space. Thus, different
  /// components sharing the same FE space are handled by the same DOF mapping.
  class FeSpaceData : public ComponentDataInterface
469
470
  {
  public:
471
472
473
    FeSpaceData()
      : iterData(this),
	iterComponent(this)
474
475
    {}

476
    /// Returns DOF mapping for a given component number.
477
478
    ComponentDofMap& operator[](int compNumber)
    {
479
480
      const FiniteElemSpace *feSpace = componentSpaces[compNumber];
      return componentData.find(feSpace)->second;
481
482
    }

483
    /// Returns DOF mapping for a given FE space.
484
485
    ComponentDofMap& operator[](const FiniteElemSpace *feSpace)
    {
486
487
488
      TEST_EXIT_DBG(componentData.count(feSpace))("No data for FE space!\n");;
      
      return componentData.find(feSpace)->second;
489
    }
490
491
492
493
494
495
496
497
498
499

    /// Checks whether the DOF mapping is defined for a specific 
    /// component number.
    bool isDefinedFor(int compNumber) const
    {
      const FiniteElemSpace *feSpace = componentSpaces[compNumber];
      return static_cast<bool>(componentData.count(feSpace));
    }

    /// Return data iterator.
500
501
    ComponentIterator& getIteratorData()
    {
502
503
      iterData.reset();
      return iterData;
504
505
    }
    
506
    /// Return component iterator.
507
508
    ComponentIterator& getIteratorComponent()
    {
509
510
      iterComponent.reset();
      return iterComponent;
511
512
    }

513
    /// Initialization
514
515
    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
516
	      bool globalMapping,
517
518
	      MeshLevelData &levelData);

519

520
  protected:
521

522
523
    /// Iterator class to iterate over all parallel DOF mappings.
    class IteratorData : public ComponentIterator {
524
    public:
525
526
      IteratorData(FeSpaceData *d)
	: data(d)
527
528
      {}

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

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

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

      void next()
      {
546
	++it;
547
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
548

549
      void reset()
Thomas Witkowski's avatar
Thomas Witkowski committed
550
      {
551
	it = data->feSpaces.begin();
Thomas Witkowski's avatar
Thomas Witkowski committed
552
553
      }

554
    protected:
555
      FeSpaceData *data;
556

557
      vector<const FiniteElemSpace*>::iterator it;
558
    };
559

560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596

    /// Iterator class to iterate over all component DOF mappings.
    class IteratorComponent : public ComponentIterator {
    public:
      IteratorComponent(FeSpaceData *d)
	: data(d)
      {}

      ComponentDofMap& operator*()
      {
	(*data)[*it];
      }

      ComponentDofMap* operator->()
      {
	&((*data)[*it]);
      }

      bool end()
      {
	return (it == data->componentSpaces.end());
      }

      void next()
      {
	++it;
      }

      void reset()
      {
	it = data->componentSpaces.begin();
      }

    protected:
      FeSpaceData *data;

      vector<const FiniteElemSpace*>::iterator it;
597
    };
598
    
599

600
    map<const FiniteElemSpace*, ComponentDofMap> componentData;
601

602
    IteratorData iterData;
603

604
605
606
607
608
    IteratorComponent iterComponent;

    friend class IteratorData;
    
    friend class IteratorComponent;
609
610
  };

611
612


613
614
615
616
617
618
  /// 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
  };
619

620
621
622
623
624
  /**
   * 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
625
  class ParallelDofMapping
626
627
  {
  public:
628
629
630
631
632
633
    /** \brief 
     * Constructur for parallel DOF mapping.
     *
     * \param[in]   mode  Defines if DOF mapping is defined either per 
     *                    component or per FE space.
     */
634
    ParallelDofMapping(DofMappingMode mode);
635

636
637
    /** \brief 
     * Initialize the parallel DOF mapping.
638
639
640
     *
     * \param[in]  m                  MPI communicator.
     * \param[in]  fe                 The FE spaces of all components of the 
641
642
643
644
     *                                PDE to be solved.
     * \param[in]  uniqueFe           Unique list of FE spaces. Thus, two
     *                                arbitrary elements of this list are always
     *                                different.
645
     * \param[in]  globalMapping      If true, at least one rank's mapping con-
646
647
     *                                taines DOFs that are not owend by the rank.
     */
648
    void init(MeshLevelData& mld,
649
	      vector<const FiniteElemSpace*> &fe,
650
	      vector<const FiniteElemSpace*> &uniqueFe,
651
	      bool globalMapping = true);
652
653
654
655

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

663
    void setMpiComm(MPI::Intracomm &m, int l);
664
665
666
    
    /// Clear all data.
    void clear();
667
668
669

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

672
673
674
    /// Returns the DOF communicator.
    DofComm& getDofComm()
    {
675
676
677
      FUNCNAME("ParallelDofMapping::getDofComm()");

      TEST_EXIT_DBG(dofComm)("No DOF communicator object defined!\n");
678
679
680
681

      return *dofComm;
    }

682
683
684
    /// Changes the computation of matrix indices based of either local or
    /// global DOF indices, see \ref needMatIndexFromGlobal
    void setComputeMatIndex(bool global)
685
686
687
688
    {
      needMatIndexFromGlobal = global;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
689
690
691
692
693
    inline bool isMatIndexFromGlobal()
    {
      return needMatIndexFromGlobal;
    }

694

695
696
697
698
    /// Access the DOF mapping for a given component number.
    inline ComponentDofMap& operator[](int compNumber)
    {
      return (*data)[compNumber];
699
700
    }

701
702
    /// Access the DOF mapping for a given FE space    
    inline ComponentDofMap& operator[](const FiniteElemSpace *feSpace) 
Thomas Witkowski's avatar
Thomas Witkowski committed
703
    {
704
      return (*data)[feSpace];
Thomas Witkowski's avatar
Thomas Witkowski committed
705
706
    }

707
708
    /// Checks whether the DOF mapping is defined for a specific 
    /// component number.
709
    inline bool isDefinedFor(int compNumber) const
Thomas Witkowski's avatar
Thomas Witkowski committed
710
    {
711
      return data->isDefinedFor(compNumber);
Thomas Witkowski's avatar
Thomas Witkowski committed
712
713
    }

714
    /// Returns the number of solution components the mapping is defined on.
Thomas Witkowski's avatar
Thomas Witkowski committed
715
716
    inline int getNumberOfComponents() const
    {      
717
      return static_cast<int>(componentSpaces.size());
Thomas Witkowski's avatar
Thomas Witkowski committed
718
    }
719

720
    /// Returns \ref nRankDofs, thus the number of DOFs owned by the rank.
721
    inline int getRankDofs()
722
    {
723
      TEST_EXIT_DBG(nRankDofs >= 0)("Should not happen!\n");
724

725
      return nRankDofs;
726
727
    }

728
    /// Returns \ref nLocalDofs, thus the number of DOFs in ranks subdomain.
729
    inline int getLocalDofs()
Thomas Witkowski's avatar
Thomas Witkowski committed
730
    {
731
      TEST_EXIT_DBG(nLocalDofs >= 0)("Should not happen!\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
732

733
      return nLocalDofs;
Thomas Witkowski's avatar
Thomas Witkowski committed
734
735
    }

736
    /// Returns \ref nOverallDofs, thus the number of all DOFs in this mapping.
737
    inline int getOverallDofs()
738
    {
739
      TEST_EXIT_DBG(nOverallDofs >= 0)("Should not happen!\n");
740

741
      return nOverallDofs;
742
743
    }

744
745
    /// Returns \ref rStartDofs, thus the smallest global index of a DOF that is
    /// owned by the rank.
746
    inline int getStartDofs()
747
    {
748
      TEST_EXIT_DBG(rStartDofs >= 0)("Should not happen!\n");
749

750
      return rStartDofs;
751
752
    }

753
    /// Update the mapping.
Thomas Witkowski's avatar
Thomas Witkowski committed
754
    void update();
755

756
757
    /// Returns the global matrix index of a given DOF for a given 
    /// component number.
758
    inline int getMatIndex(int ithComponent, DegreeOfFreedom d)
Thomas Witkowski's avatar
Thomas Witkowski committed
759
    {
760
      return dofToMatIndex.get(ithComponent, d);
Thomas Witkowski's avatar
Thomas Witkowski committed
761
762
    }

763
764
765
766
767
768
769
    /// 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);
    }

770
771
    /// Returns the local matrix index of a given DOF for a given 
    /// component number.
772
    inline int getLocalMatIndex(int ithComponent, DegreeOfFreedom d)
773
    {
774
      return dofToMatIndex.get(ithComponent, d) - rStartDofs;
775
776
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
777
778
779
780
781
782
783
784
785
786
787
788
    // 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");
    }

789
    /// Compute local and global matrix indices.
790
    void computeMatIndex(bool globalIndex);
791

792
793
794
795
    void createIndexSet(IS &is, 
			int firstComponent, 
			int nComponents);

796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
    /// 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
811
812
813
814
815
    inline void createLocalVec(Vec &vec)
    {
      VecCreateSeq(PETSC_COMM_SELF, getRankDofs(), &vec);
    }

816
817
  protected:
    /// Compute \ref nRankDofs.
818
    int computeRankDofs();
819
820

    /// Compute \ref nLocalDofs.
821
    int computeLocalDofs();
822
823

    /// Compute \ref nOverallDofs.
824
    int computeOverallDofs();
825
826

    /// Compute \ref rStartDofs.
827
    int computeStartDofs();
828

829
  private:
830
831
832
833
    MPI::Intracomm mpiComm;

    int meshLevel;

834
    MeshLevelData *levelData;
835

836
837
    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;
838

839
840
841
    /// 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.
842
    bool globalMapping;
843
844
845
846
847

    /// 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
848
849
850
851
852
    /// 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.
853
854
    bool needMatIndexFromGlobal;

855
856
    /// Maps from components to DOF mappings.
    ComponentDataInterface *data;
Thomas Witkowski's avatar
Thomas Witkowski committed
857

858
    /// Number of DOFs owned by rank.
859
    int nRankDofs;
860
861

    /// Number of DOFs in rank's subdomain.
862
    int nLocalDofs;
863
864

    /// Number of global DOFs (this value is thus the same on all ranks).
865
    int nOverallDofs;
866
867

    /// Smallest global index of a DOF owned by the rank.
868
    int rStartDofs;
869
870
871

    /// Mapping from global DOF indices to global matrix indices under 
    /// consideration of possibly multiple components.
872
    DofToMatIndex dofToMatIndex;
873
874
875

    /// FE spaces of all components
    vector<const FiniteElemSpace*> componentSpaces;
876
877
878
879
  };
}

#endif