ZeroOrderAssembler.cc 9.64 KB
Newer Older
1
#include <vector>
2
#include <boost/numeric/mtl/mtl.hpp>
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "Assembler.h"
#include "ZeroOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "OpenMP.h"

namespace AMDiS {

  std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;

  ZeroOrderAssembler::ZeroOrderAssembler(Operator *op,
					 Assembler *assembler,
					 Quadrature *quad,
					 bool optimized)
    : SubAssembler(op, assembler, quad, 0, optimized)
  {}

24

25
26
27
28
29
  ZeroOrderAssembler* ZeroOrderAssembler::getSubAssembler(Operator* op,
							  Assembler *assembler,
							  Quadrature *quad,
							  bool optimized)
  {
30
31
    int myRank = omp_get_thread_num();

32
33
    // check if an assembler is needed at all
    if (op->zeroOrder[myRank].size() == 0)
34
      return NULL;   
35
36
37
38

    ZeroOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
39
      optimized ? &optimizedSubAssemblers : &standardSubAssemblers;
40

41
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
42
43
44
45
46
47
48
49
50
51

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    if (quad) {
      for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
	std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());

	sort(assTerms.begin(), assTerms.end());

52
	if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
53
54
55
56
57
58
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
      }
    }
 
    // check if all terms are pw_const
    bool pwConst = true;
59
    for (int i = 0; i < static_cast<int>(op->zeroOrder[myRank].size()); i++) {
60
61
62
63
64
65
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

66
67
68
    newAssembler = new StandardZOA(op, assembler, quad);


69
    // create new assembler
70
    if (!optimized) {
Thomas Witkowski's avatar
Thomas Witkowski committed
71
      newAssembler = new StandardZOA(op, assembler, quad);
72
    } else {
73
      if (pwConst)
74
  	newAssembler = new PrecalcZOA(op, assembler, quad);
75
      else
76
  	newAssembler = new FastQuadZOA(op, assembler, quad);      
77
    }
78
79
80
81
82

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

83

84
  StandardZOA::StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad)
85
86
87
88
89
    : ZeroOrderAssembler(op, assembler, quad, false)      
  {
    name = "standard zero order assembler";
  }

90

91
  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
92
  {
93
94
    const BasisFunction *psi = rowFeSpace->getBasisFcts();
    const BasisFunction *phi = colFeSpace->getBasisFcts();
95
96

    int nPoints = quadrature->getNumPoints();
97
98
    std::vector<double> c(nPoints, 0.0);
    std::vector<double> phival(nCol);
99
100
101

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
102
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
103
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
104
    
105
    if (symmetric) {
106
107
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
            
108
109
110
111
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
112
	for (int i = 0; i < nCol; i++)
113
114
115
116
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
117
	  mat[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
118
119
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
120
121
	    mat[i][j] += val;
	    mat[j][i] += val;
122
123
124
125
126
127
128
129
	  }
	}
      }
    } else {      //  non symmetric assembling 
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
130
	for (int i = 0; i < nCol; i++)
131
132
133
134
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
135
	  for (int j = 0; j < nCol; j++)
136
	    mat[i][j] += quadrature->getWeight(iq) * c[iq] * psival * phival[j];
137
138
139
140
141
	}
      }
    }
  }

142

143
  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
144
145
  {
    int nPoints = quadrature->getNumPoints();
146
    std::vector<double> c(nPoints, 0.0);
147
    int myRank = omp_get_thread_num();
148

149
    std::vector<OperatorTerm*>::iterator termIt;
150
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
151
152
153
154
155
156
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
157
	double psi = (*(rowFeSpace->getBasisFcts()->getPhi(i)))
158
	  (quadrature->getLambda(iq));
159
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi;
160
161
162
163
      }
    }
  }

164

165
166
  FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
167
  {
168
    name = "fast quadrature zero order assembler";
169
170
    tmpC.resize(omp_get_overall_max_threads());
  }
171

172

173
  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
174
175
176
177
178
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
179
180
181
182
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
183
	const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
184
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
185
	basFcts = colFeSpace->getBasisFcts();
186
187
188
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
189
190
    }

191
    std::vector<double> &c = tmpC[myRank];
Thomas Witkowski's avatar
Thomas Witkowski committed
192
    if (static_cast<int>(c.size()) < nPoints)
193
      c.resize(nPoints);
Thomas Witkowski's avatar
Thomas Witkowski committed
194
    for (int i = 0; i < nPoints; i++)
195
196
197
198
      c[i] = 0.0;

    for (std::vector<OperatorTerm*>::iterator termIt = terms[myRank].begin(); 
	 termIt != terms[myRank].end(); ++termIt)
199
200
201
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);

    if (symmetric) {
202
203
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

204
      for (int iq = 0; iq < nPoints; iq++) {
205
	c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
206
207
	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
208

209
	for (int i = 0; i < nRow; i++) {
210
	  mat[i][i] += c[iq] * psi[i] * phi[i];
211
	  for (int j = i + 1; j < nCol; j++) {
212
	    double val = c[iq] * psi[i] * phi[j];
213
214
	    mat[i][j] += val;
	    mat[j][i] += val;
215
216
217
218
219
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
220
	c[iq] *= elInfo->getDet() * quadrature->getWeight(iq);
221
222
223

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
224
225
226
	for (int i = 0; i < nRow; i++)
	  for (int j = 0; j < nCol; j++)
	    mat[i][j] += c[iq] * psi[i] * phi[j];
227
228
229
      }
    }
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
230

231

232
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
233
  {
234
235
236
    int myRank = omp_get_thread_num();
    int nPoints = quadrature->getNumPoints();

237
    if (firstCall) {
238
239
240
241
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
242
	const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
243
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
244
	basFcts = colFeSpace->getBasisFcts();
245
246
247
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
248
249
    }

250
    std::vector<double> c(nPoints, 0.0);
251
252
253
254
255
256
257
258
259
260
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= elInfo->getDet();

      const double *psi = psiFast->getPhi(iq);
      for (int i = 0; i < nRow; i++) {
261
	vec[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
262
263
264
265
      }
    }
  }

266

267
268
269
  PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
270
    name = "precalculated zero order assembler";
271
272
  }

273

274
  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
275
276
  {
    if (firstCall) {
277
278
279
280
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
281
282
	q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), 
					  colFeSpace->getBasisFcts(), 
283
					  quadrature);
284
	q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);
285
286
	firstCall = false;
      }
287
288
    }

289
    std::vector<double> c(1, 0.0);
290
291
292
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

293
294
    for (int i = 0; i < size; i++)
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);    
295

296
    c[0] *= elInfo->getDet();
297
    
298
299
    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
300
	mat[i][i] += c[0] * q00->getValue(i,i);
301
	for (int j = i + 1; j < nCol; j++) {
302
	  double val = c[0] * q00->getValue(i, j);
303
304
	  mat[i][j] += val;
	  mat[j][i] += val;
305
306
307
	}
      }
    } else {
308
      for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
309
	for (int j = 0; j < nCol; j++) {
310
	  mat[i][j] += c[0] * q00->getValue(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
311
	}
312
      }
313
314
315
    }
  }

316

317
  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
318
319
  {
    if (firstCall) {
320
321
322
323
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
324
325
	q00 = Q00PsiPhi::provideQ00PsiPhi(rowFeSpace->getBasisFcts(), 
					  colFeSpace->getBasisFcts(), 
326
					  quadrature);
327
	q0 = Q0Psi::provideQ0Psi(rowFeSpace->getBasisFcts(), quadrature);	
328
329
	firstCall = false;
      }
330
331
332
333
334
    }

    std::vector<OperatorTerm*>::iterator termIt;

    int myRank = omp_get_thread_num();
335
    std::vector<double> c(1, 0.0);
336
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
337
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
338

339
    c[0] *= elInfo->getDet();
340
341

    for (int i = 0; i < nRow; i++)
342
      vec[i] += c[0] * q0->getValue(i);
343
344
345
  }

}