Operator.cc 4.59 KB
Newer Older
1
2
3
4
5
6
#include "Operator.h"
#include "ElInfo.h"
#include "Assembler.h"
#include "FixVec.h"
#include "DOFVector.h"
#include "Quadrature.h"
7
#include "OpenMP.h"
8
9
10
11
12
13
14
15
16

namespace AMDiS {

  const Flag Operator::MATRIX_OPERATOR = 1;
  const Flag Operator::VECTOR_OPERATOR = 2;


  int Operator::getQuadratureDegree(int order, FirstOrderType firstOrderType) 
  {
17
    std::vector<OperatorTerm*>* terms = NULL;
18
    int myRank = omp_get_thread_num();
19
20
21

    switch(order) {
    case 0:
22
      terms = &zeroOrder[myRank];
23
24
25
      break;
    case 1:
      if (firstOrderType == GRD_PHI)
26
	terms = &firstOrderGrdPhi[myRank];
27
      else 
28
	terms = &firstOrderGrdPsi[myRank];
29
30
      break;
    case 2:
31
      terms = &secondOrder[myRank];
32
33
34
35
36
37
38
39
40
41
      break;
    }

    const BasisFunction *psi = rowFESpace->getBasisFcts();
    const BasisFunction *phi = colFESpace->getBasisFcts();

    int psiDegree = psi->getDegree();
    int phiDegree = phi->getDegree();
    int maxTermDegree = 0;

Thomas Witkowski's avatar
Thomas Witkowski committed
42
    for (int i = 0; i < static_cast<int>(terms->size()); i++) {
43
44
45
      OperatorTerm *term = (*terms)[i];
      maxTermDegree = max(maxTermDegree, term->degree);
    }
46
47
   
    return psiDegree + phiDegree - order + maxTermDegree;
48
49
50
51
  }


  Operator::Operator(Flag operatorType,
Thomas Witkowski's avatar
Thomas Witkowski committed
52
53
54
55
		     const FiniteElemSpace *row,
		     const FiniteElemSpace *col)
    : rowFESpace(row), 
      colFESpace(col ? col : row),
56
57
58
      type(operatorType), 
      fillFlag(Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS | 
	       Mesh::FILL_DET | Mesh::FILL_GRD_LAMBDA),
Thomas Witkowski's avatar
Thomas Witkowski committed
59
      needDualTraverse(false),
60
61
62
      uhOld(NULL),
      optimized(true)
  {
63
    int maxThreads = omp_get_overall_max_threads();
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78

    assembler.resize(maxThreads);
    secondOrder.resize(maxThreads);
    firstOrderGrdPsi.resize(maxThreads);
    firstOrderGrdPhi.resize(maxThreads);
    zeroOrder.resize(maxThreads);

    for (int i = 0; i < maxThreads; i++) {
      assembler[i] = NULL;
      secondOrder[i].resize(0);
      firstOrderGrdPsi[i].resize(0);
      firstOrderGrdPhi[i].resize(0);
      zeroOrder[i].resize(0);
    }

79
80
    nRow = rowFESpace->getBasisFcts()->getNumber();
    nCol = colFESpace->getBasisFcts()->getNumber();
Thomas Witkowski's avatar
Thomas Witkowski committed
81
  }
82

83
84

  void Operator::setUhOld(const DOFVectorBase<double> *vec)
85
  {
86
87
    uhOld = vec;
    auxFeSpaces.insert(vec->getFESpace());
88
89
  }

90

91
  void Operator::getElementMatrix(const ElInfo *elInfo, 
92
				  ElementMatrix& userMat, 
93
94
				  double factor)
  {
95
96
    int myRank = omp_get_thread_num();

97
    if (!assembler[myRank])
98
      initAssembler(myRank, NULL, NULL, NULL, NULL);
99

100
    assembler[myRank]->calculateElementMatrix(elInfo, userMat, factor);
101
102
  }

103

Thomas Witkowski's avatar
Thomas Witkowski committed
104
105
  void Operator::getElementMatrix(const ElInfo *rowElInfo, const ElInfo *colElInfo,
				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
106
				  ElementMatrix& userMat, 
107
108
109
110
				  double factor)
  {
    int myRank = omp_get_thread_num();

111
    if (!assembler[myRank])
112
113
114
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementMatrix(rowElInfo, colElInfo, 
115
					      smallElInfo, largeElInfo,
116
117
118
    					      userMat, factor);
  }

119

120
  void Operator::getElementVector(const ElInfo *elInfo, 
121
				  ElementVector& userVec, 
122
123
				  double factor)
  {
124
125
    int myRank = omp_get_thread_num();

126
    if (!assembler[myRank])
127
      initAssembler(myRank, NULL, NULL, NULL, NULL);
128

129
    assembler[myRank]->calculateElementVector(elInfo, userVec, factor);
130
131
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
132
133
  void Operator::getElementVector(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				  const ElInfo *smallElInfo, const ElInfo *largeElInfo,
134
				  ElementVector& userVec,
Thomas Witkowski's avatar
Thomas Witkowski committed
135
136
137
138
				  double factor)
  {
    int myRank = omp_get_thread_num();

139
    if (!assembler[myRank])
Thomas Witkowski's avatar
Thomas Witkowski committed
140
141
142
143
144
145
146
      initAssembler(myRank, NULL, NULL, NULL, NULL);

    assembler[myRank]->calculateElementVector(mainElInfo, auxElInfo, 
					      smallElInfo, largeElInfo,
					      userVec, factor);
  }

147
148
  void Operator::initAssembler(int rank,
			       Quadrature *quad2,
149
150
151
			       Quadrature *quad1GrdPsi,
			       Quadrature *quad1GrdPhi,
			       Quadrature *quad0) 
152
  {    
153
154
155
156
#ifdef _OPENMP
#pragma omp critical (initAssembler)
#endif
      {
157
158
159
160
161
162
163
164
	if (optimized)
	  assembler[rank] = 
	    new OptimizedAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
				   rowFESpace, colFESpace);
	else
	  assembler[rank] = 
	    new StandardAssembler(this, quad2, quad1GrdPsi, quad1GrdPhi, quad0,
				  rowFESpace, colFESpace);	
165
166
167
      }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
168
169
170
171
172
  void Operator::finishAssembling()
  {
    assembler[omp_get_thread_num()]->finishAssembling();
  }

173

174
  Assembler *Operator::getAssembler(int rank) 
Thomas Witkowski's avatar
Thomas Witkowski committed
175
  {
176
177
    if (!assembler[rank])
      initAssembler(rank, NULL, NULL, NULL, NULL);    
Thomas Witkowski's avatar
Thomas Witkowski committed
178

179
    return assembler[rank];
Thomas Witkowski's avatar
Thomas Witkowski committed
180
181
  }

182
  void Operator::setAssembler(int rank, Assembler *a)
Thomas Witkowski's avatar
Thomas Witkowski committed
183
  {
184
    assembler[rank] = a; 
185
186
187
  }
 
}