FirstOrderAssembler.cc 13.9 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
24
25
26
#include "Assembler.h"
#include "FirstOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "OpenMP.h"

namespace AMDiS {

  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;

  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
  std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;


  FirstOrderAssembler::FirstOrderAssembler(Operator *op,
					   Assembler *assembler,
					   Quadrature *quad,
					   bool optimized,
					   FirstOrderType type)
    : SubAssembler(op, assembler, quad, 1, optimized, type)
27
28
29
30
31
  {
    VectorOfFixVecs<DimVec<double> > 
      Lb(assembler->getRowFESpace()->getMesh()->getDim(), 1, NO_INIT);
    tmpLb.resize(omp_get_overall_max_threads(), Lb);
  }
32
33
34
35
36
37
38
39
40

  FirstOrderAssembler* 
  FirstOrderAssembler::getSubAssembler(Operator* op,
				       Assembler *assembler,
				       Quadrature *quad,
				       FirstOrderType type,
				       bool optimized)
  {
    std::vector<SubAssembler*> *subAssemblers =
41
42
43
44
45
46
47
48
      optimized ?
      (type == GRD_PSI ? 
       &optimizedSubAssemblersGrdPsi : 
       &optimizedSubAssemblersGrdPhi) :
      (type == GRD_PSI ? 
       &standardSubAssemblersGrdPsi :
       &standardSubAssemblersGrdPhi);
    
49
50
    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*> opTerms 
51
      = (type == GRD_PSI) ? op->firstOrderGrdPsi[myRank] : op->firstOrderGrdPhi[myRank];
52
53

    // check if a assembler is needed at all
54
    if (opTerms.size() == 0)
55
56
57
58
59
60
61
      return NULL;

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

    FirstOrderAssembler *newAssembler;

    // check if a new assembler is needed
62
    for (int i = 0; i < static_cast<int>(subAssemblers->size()); i++) {
63
64
65
66
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());   

67
      if (opTerms == assTerms && (*subAssemblers)[i]->getQuadrature() == quad)
68
69
70
71
72
	return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
    }

    // check if all terms are pw_const
    bool pwConst = true;
73
    for (int i = 0; i < static_cast<int>(opTerms.size()); i++) {
74
75
76
77
78
79
80
      if (!(opTerms[i])->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
81
    if (!optimized) {
82
83
      newAssembler = 
	(type == GRD_PSI) ?
Thomas Witkowski's avatar
Thomas Witkowski committed
84
85
	dynamic_cast<FirstOrderAssembler*>(new Stand10(op, assembler, quad)) :
	dynamic_cast<FirstOrderAssembler*>(new Stand01(op, assembler, quad));    
86
87
88
89
    } else {
       if (pwConst) {
 	newAssembler = 
 	  (type == GRD_PSI) ?
Thomas Witkowski's avatar
Thomas Witkowski committed
90
91
 	  dynamic_cast<FirstOrderAssembler*>(new Pre10(op, assembler, quad)) :
 	  dynamic_cast<FirstOrderAssembler*>(new Pre01(op, assembler, quad));
92
93
94
       } else {
 	newAssembler = 
 	  (type == GRD_PSI) ?
Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
 	  dynamic_cast<FirstOrderAssembler*>(new Quad10(op, assembler, quad)) :
 	  dynamic_cast<FirstOrderAssembler*>(new Quad01(op, assembler, quad));
97
98
       }
    }
99
100
101

    subAssemblers->push_back(newAssembler);
    return newAssembler;
102
  }
103
104
105
106
107

  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
  {}

108
  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
109
110
111
112
113
114
115
  {
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();
116
117
118
119
    int myRank = omp_get_thread_num();
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);
    std::vector<double> phival(nCol);
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nCol; i++) {
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
      }

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	for (int j = 0; j < nCol; j++) {
138
	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
139
140
141
142
143
	}
      }
    }
  }

144
  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
  {
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
165
	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
166
167
168
169
      }
    }
  }

170
171
172

  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
173
  {}
174
175


176
  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
177
178
179
180
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

    if (firstCall) {
181
182
183
184
185
186
187
188
189
190
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
	basFcts = owner->getColFESpace()->getBasisFcts();
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
191
192
193
    }

    int nPoints = quadrature->getNumPoints();
194
195
196
    int myRank = omp_get_thread_num();
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);
197
198
199
200
201
202
203
204
205
206
207
208

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

209
      const double *phi = phiFast->getPhi(iq);
210
211
212
213
      grdPsi = psiFast->getGradient(iq);

      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
214
	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
215
216
217
218
      }
    }
  }

219
  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

    if (firstCall) {
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
	basFcts = owner->getColFESpace()->getBasisFcts();
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
	firstCall = false;
      }
    }

    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {

      Lb[iq] *= elInfo->getDet();
      grdPsi = psiFast->getGradient(iq);

      for (int i = 0; i < nRow; i++) {
254
	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
255
256
257
      }
    }
  }
258
259
260
261
262
263
264

  Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }


265
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
266
267
268
269
270
  {
    const int *k;
    const double *values;

    if (firstCall) {
271
272
273
274
275
276
277
278
279
280
281
#ifdef _OPENMP
#pragma omp critical
#endif 
      {	
	q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);
	firstCall = false;
      }
282
283
284
285
    }

    const int **nEntries = q10->getNumberEntries();
    int myRank = omp_get_thread_num();
286
287
    // Do not need do resize Lb, because it's size is always at least one.
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
    Lb[0].set(0.0);

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
    }

    Lb[0] *= elInfo->getDet();

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	k = q10->getKVec(i, j);
	values = q10->getValVec(i, j);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][k[m]];
	}
304
	mat[i][j] += val;
305
306
307
308
309
310
311
312
      }
    }
  }

  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
  {}

313
  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
314
315
316
317
318
319
320
321
  {
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
322
323
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      for (int i = 0; i < nCol; i++) {
	(*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
      }

      for (int i = 0; i < nRow; i++) {
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	for (int j = 0; j < nCol; j++)
342
	  mat[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
343
344
345
346
347
348
      }
    } 
  }

  Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
349
  {}
350

351
  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
352
353
354
355
  {
    VectorOfFixVecs<DimVec<double> > *grdPhi;

    if (firstCall) {
356
357
358
359
360
361
362
363
364
365
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
	psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
	basFcts = owner->getColFESpace()->getBasisFcts();
	phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
	firstCall = false;
      }
366
367
368
369
    }

    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
370
371
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387

    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq].set(0.0);
    }
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
    }
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

      const double *psi = psiFast->getPhi(iq);
      grdPhi = phiFast->getGradient(iq);

      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
388
	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
389
390
391
392
393
394
395
396
397
      }
    }
  }

  Pre01::Pre01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
  {
  }

398
  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
399
400
401
402
403
  {
    const int *l;
    const double *values;

    if (firstCall) {
404
405
406
407
408
409
410
411
412
413
414
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q01 = Q01PsiPhi::provideQ01PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);
	firstCall = false;
      }
415
416
417
418
    }

    const int **nEntries = q01->getNumberEntries();
    int myRank = omp_get_thread_num();
419
420
    // Do not need to resize Lb, because it's size is always at least one!
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
    Lb[0].set(0.0);

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
    }

    Lb[0] *= elInfo->getDet();

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	l = q01->getLVec(i, j);
	values = q01->getValVec(i, j);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][l[m]];
	}
437
	mat[i][j] += val;
438
439
440
441
      }
    }
  }

442
  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
443
444
445
446
447
  {
    const int *k;
    const double *values;

    if (firstCall) {
448
449
450
451
452
453
454
455
456
457
458
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);
	firstCall = false;
      }
459
460
461
462
    }

    const int *nEntries = q1->getNumberEntries();
    int myRank = omp_get_thread_num();
463
464
    // Do not need to resize Lb, because it's size is always at least one!
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
    Lb[0].set(0.0);

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>(terms[myRank][i]))->getLb(elInfo, 1, Lb);
    }

    Lb[0] *= elInfo->getDet();

    for (int i = 0; i < nRow; i++) {
      k = q1->getKVec(i);
      values = q1->getValVec(i);
      double val = 0.0;
      for (int m = 0; m < nEntries[i]; m++) {
	val += values[m] * Lb[0][k[m]];
      }
480
      vec[i] += val;
481
482
483
484
    }
  }

}