FirstOrderAssembler.cc 13.8 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
    int myRank = omp_get_thread_num();
Thomas Witkowski's avatar
Thomas Witkowski committed
50
51
    std::vector<OperatorTerm*> opTerms = 
      (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

  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
Thomas Witkowski's avatar
Thomas Witkowski committed
106
107
108
109
  {
    psi = owner->getRowFESpace()->getBasisFcts();
    phi = owner->getColFESpace()->getBasisFcts();
  }
110

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

    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
    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);    
125
126
127
128
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

129
130
      for (int i = 0; i < nCol; i++)
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));      
131
132
133

      for (int i = 0; i < nRow; i++) {
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
134
	for (int j = 0; j < nCol; j++)
Thomas Witkowski's avatar
Thomas Witkowski committed
135
	  mat[i][j] += quadrature->getWeight(iq) * phival[j] * (Lb[iq] * grdPsi);
136
137
138
139
      }
    }
  }

140
  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
141
142
143
144
145
146
147
  {
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
    Lb.resize(nPoints);

148
149
150
151
    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);    
152
153
154
155
156
157
  
    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);
158
	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
159
160
161
162
      }
    }
  }

163
164
165

  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
166
  {}
167
168


169
  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
170
171
172
173
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

    if (firstCall) {
174
175
176
177
178
179
180
181
182
183
#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;
      }
184
185
186
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
191
    for (int iq = 0; iq < nPoints; iq++)
192
193
      Lb[iq].set(0.0);

Thomas Witkowski's avatar
Thomas Witkowski committed
194
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
195
196
197
198
199
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, nPoints, Lb);
  
    for (int iq = 0; iq < nPoints; iq++) {
      Lb[iq] *= elInfo->getDet();

200
      const double *phi = phiFast->getPhi(iq);
201
202
      grdPsi = psiFast->getGradient(iq);

Thomas Witkowski's avatar
Thomas Witkowski committed
203
      for (int i = 0; i < nRow; i++)
204
	for (int j = 0; j < nCol; j++)
205
	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
206
207
208
    }
  }

209
  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
  {
    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);

Thomas Witkowski's avatar
Thomas Witkowski committed
231
232
233
    for (int iq = 0; iq < nPoints; iq++)
      Lb[iq].set(0.0);    
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
234
235
236
237
238
239
      (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);

Thomas Witkowski's avatar
Thomas Witkowski committed
240
      for (int i = 0; i < nRow; i++)
241
	vec[i] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]);
242
243
    }
  }
244
245
246
247
248
249
250

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


251
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
252
253
254
255
256
  {
    const int *k;
    const double *values;

    if (firstCall) {
257
258
259
260
261
262
263
264
265
266
267
#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;
      }
268
269
270
271
    }

    const int **nEntries = q10->getNumberEntries();
    int myRank = omp_get_thread_num();
272
273
    // Do not need do resize Lb, because it's size is always at least one.
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
274
275
276
277
278
279
280
281
282
283
284
285
286
    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;
287
288
	for (int m = 0; m < nEntries[i][j]; m++)
	  val += values[m] * Lb[0][k[m]];	
289
	mat[i][j] += val;
290
291
292
293
294
295
      }
    }
  }

  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
296
297
298
299
  {
    VectorOfFixVecs<DimVec<double> > 
      grdPhi(assembler->getRowFESpace()->getMesh()->getDim(), nCol, NO_INIT);
    tmpGrdPhi.resize(omp_get_overall_max_threads(), grdPhi);
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
302

    psi = owner->getRowFESpace()->getBasisFcts();
    phi = owner->getColFESpace()->getBasisFcts();
303
  }
304

305
  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
306
307
308
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
309
310
    VectorOfFixVecs<DimVec<double> >& Lb = tmpLb[myRank];
    VectorOfFixVecs<DimVec<double> >& grdPhi = tmpGrdPhi[myRank];
311
    Lb.resize(nPoints);
312

313
    for (int iq = 0; iq < nPoints; iq++)
314
      Lb[iq].set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
315

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

322
      for (int i = 0; i < nCol; i++)
323
324
325
326
327
	(*(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++)
Thomas Witkowski's avatar
Thomas Witkowski committed
328
	  mat[i][j] += quadrature->getWeight(iq) * psival * (Lb[iq] * grdPhi[j]);
329
330
331
332
333
334
      }
    } 
  }

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

337
  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
338
339
340
341
  {
    VectorOfFixVecs<DimVec<double> > *grdPhi;

    if (firstCall) {
342
343
344
345
346
347
348
349
350
351
#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;
      }
352
353
354
355
    }

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

359
    for (int iq = 0; iq < nPoints; iq++)
360
      Lb[iq].set(0.0);
361
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++)
362
363
364
365
366
367
368
369
      (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);

370
      for (int i = 0; i < nRow; i++)
371
	for (int j = 0; j < nCol; j++)
372
	  mat[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
373
374
375
376
377
378
379
380
    }
  }

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

381
  void Pre01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
382
383
384
385
386
  {
    const int *l;
    const double *values;

    if (firstCall) {
387
388
389
390
391
392
393
394
395
396
397
#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;
      }
398
399
400
401
    }

    const int **nEntries = q01->getNumberEntries();
    int myRank = omp_get_thread_num();
402
403
    // Do not need to resize Lb, because it's size is always at least one!
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
    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]];
	}
420
	mat[i][j] += val;
421
422
423
424
      }
    }
  }

425
  void Pre10::calculateElementVector(const ElInfo *elInfo, ElementVector& vec)
426
427
428
429
430
  {
    const int *k;
    const double *values;

    if (firstCall) {
431
432
433
434
435
436
437
438
439
440
441
#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;
      }
442
443
444
445
    }

    const int *nEntries = q1->getNumberEntries();
    int myRank = omp_get_thread_num();
446
447
    // Do not need to resize Lb, because it's size is always at least one!
    VectorOfFixVecs<DimVec<double> > &Lb = tmpLb[myRank];
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
    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]];
      }
463
      vec[i] += val;
464
465
466
467
    }
  }

}