ParallelDofMapping.h 24.1 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
155
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)
	("DOF %d is already defined in mapping!\n", dof0);
156
      
157
158
      dofMap[dof0].local = dof1;
      nLocalDofs++;
Thomas Witkowski's avatar
Thomas Witkowski committed
159
      if (dof1 != -1)
160
	nRankDofs++;
161
162
    }
    
163
164
    /// 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.
165
    void insertNonRankDof(DegreeOfFreedom dof0, DegreeOfFreedom dof1 = -1)
166
    {
167
      FUNCNAME("ComponentDofMap::insertNonRankDof()");
168
      
169
      TEST_EXIT_DBG(dofMap.count(dof0) == 0)("Should not happen!\n");
170
      
171
172
173
      dofMap[dof0].local = dof1;
      nLocalDofs++;
      nonRankDofs.insert(dof0);
174
175
    }
    
176
    /// Checks if a given DOF is in the DOF mapping.
177
    bool isSet(DegreeOfFreedom dof)
178
    {
179
      return static_cast<bool>(dofMap.count(dof));
180
    }
181

Thomas Witkowski's avatar
Blbu    
Thomas Witkowski committed
182
183
184
    /// 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.
185
    bool isRankDof(DegreeOfFreedom dof)
186
    {
187
      return !(static_cast<bool>(nonRankDofs.count(dof)));
188
    }
Thomas Witkowski's avatar
Blbu    
Thomas Witkowski committed
189
190
191
192
193

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

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

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

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

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

226
    /// Informs the mapping whether a global index must be computed.
227
    void setGlobalMapping(bool b)
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
228
    {
229
      globalMapping = b;
Thomas Witkowski's avatar
FETI-DP    
Thomas Witkowski committed
230
231
    }

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

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

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

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

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

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

258
259
260
261
    MPI::Intracomm mpiComm;

    int meshLevel;

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

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

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

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

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

280

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

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

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

    virtual void reset() = 0;
  };

296

297
298
299
  /// 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.
300
301
302
  class ComponentDataInterface
  {
  public:
303
    /// Access via component number
304
305
    virtual ComponentDofMap& operator[](int compNumber) = 0;

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

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

313
    /// Returns iterator which iterates over all DOF mappings.
314
315
    virtual ComponentIterator& getIteratorData() = 0;
    
316
317
318
319
    /// 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.
320
321
    virtual ComponentIterator& getIteratorComponent() = 0;

322
323
324
325
326
327
328
329
330
331
332
333
    /** \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,
334
		      MeshLevelData &levelData) = 0;
335
336
337

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

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


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

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

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

364
365
366
367
368
    /// 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)
369
    {
370
      ERROR_EXIT("FE Space access is not possible for component wise defined DOF mappings\n");
371
    }
372
373
    
    /// Return data iterator.
374
375
    ComponentIterator& getIteratorData()
    {
376
377
      iter.reset();
      return iter;
378
379
    }
    
380
    /// Return component iterator.
381
382
    ComponentIterator& getIteratorComponent()
    {
383
384
      iter.reset();
      return iter;
385
386
    }

387
388
389
390
391
392
393
394
    /// 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
395
396
    void init(vector<const FiniteElemSpace*> &f0,
	      vector<const FiniteElemSpace*> &f1,
397
	      bool globalMapping,
398
	      MeshLevelData &levelData);
399

400
  protected:
401

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

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

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

420

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

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

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

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

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

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

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


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

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


460

461
    friend class Iterator;
462
463
464
  };


465
466
467
468
469

  /// 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
470
471
  {
  public:
472
473
474
    FeSpaceData()
      : iterData(this),
	iterComponent(this)
475
476
    {}

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

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

    /// 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.
501
502
    ComponentIterator& getIteratorData()
    {
503
504
      iterData.reset();
      return iterData;
505
506
    }
    
507
    /// Return component iterator.
508
509
    ComponentIterator& getIteratorComponent()
    {
510
511
      iterComponent.reset();
      return iterComponent;
512
513
    }

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

520

521
  protected:
522

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

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

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

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

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

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

555
    protected:
556
      FeSpaceData *data;
557

558
      vector<const FiniteElemSpace*>::iterator it;
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
597

    /// 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;
598
    };
599
    
600

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

603
    IteratorData iterData;
604

605
606
607
608
609
    IteratorComponent iterComponent;

    friend class IteratorData;
    
    friend class IteratorComponent;
610
611
  };

612
613


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

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

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

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

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

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

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

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

      return *dofComm;
    }

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

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

695

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

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

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

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

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

726
      return nRankDofs;
727
728
    }

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

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

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

742
      return nOverallDofs;
743
744
    }

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

751
      return rStartDofs;
752
753
    }

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

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

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

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

778
779
780
781
782
783
    /// Returns the set of unique FE spaces.
    inline vector<const FiniteElemSpace*>& getFeSpaces()
    {
      return feSpaces;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
784
785
786
787
788
789
790
791
792
793
794
795
    // 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");
    }

796
    /// Compute local and global matrix indices.
797
    void computeMatIndex(bool globalIndex);
798

799
800
801
802
    void createIndexSet(IS &is, 
			int firstComponent, 
			int nComponents);

803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
    /// 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
818
819
820
821
822
    inline void createLocalVec(Vec &vec)
    {
      VecCreateSeq(PETSC_COMM_SELF, getRankDofs(), &vec);
    }

823
824
  protected:
    /// Compute \ref nRankDofs.
825
    int computeRankDofs();
826
827

    /// Compute \ref nLocalDofs.
828
    int computeLocalDofs();
829
830

    /// Compute \ref nOverallDofs.
831
    int computeOverallDofs();
832
833

    /// Compute \ref rStartDofs.
834
    int computeStartDofs();
835

836
  private:
837
838
839
840
    MPI::Intracomm mpiComm;

    int meshLevel;

841
    MeshLevelData *levelData;
842

843
844
    /// DOF communicator for all DOFs on interior boundaries.
    DofComm *dofComm;
845

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

    /// 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
855
856
857
858
859
    /// 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.
860
861
    bool needMatIndexFromGlobal;

862
863
    /// Maps from components to DOF mappings.
    ComponentDataInterface *data;
Thomas Witkowski's avatar
Thomas Witkowski committed
864

865
    /// Number of DOFs owned by rank.
866
    int nRankDofs;
867
868

    /// Number of DOFs in rank's subdomain.
869
    int nLocalDofs;
870
871

    /// Number of global DOFs (this value is thus the same on all ranks).
872
    int nOverallDofs;
873
874

    /// Smallest global index of a DOF owned by the rank.
875
    int rStartDofs;
876
877
878

    /// Mapping from global DOF indices to global matrix indices under 
    /// consideration of possibly multiple components.
879
    DofToMatIndex dofToMatIndex;
880

881
    /// FE spaces of all components.
882
    vector<const FiniteElemSpace*> componentSpaces;
883
884
885

    /// Set of unique FE spaces.
    vector<const FiniteElemSpace*> feSpaces;
886
887
888
889
  };
}

#endif