ZeroOrderAssembler.cc 12.9 KB
Newer Older
1
#include <vector>
2
3
#include <boost/numeric/mtl/mtl.hpp>

4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
#include "Assembler.h"
#include "ZeroOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "ElementMatrix.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)
  {}

  ZeroOrderAssembler* ZeroOrderAssembler::getSubAssembler(Operator* op,
							  Assembler *assembler,
							  Quadrature *quad,
							  bool optimized)
  {
    // check if a assembler is needed at all
32
33
    if (op->zeroOrder.size() == 0)
      return NULL;   
34
35
36
37

    ZeroOrderAssembler *newAssembler;

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

    int myRank = omp_get_thread_num();
41
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68

    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());

	if ((opTerms == assTerms) && 
	    ((*subAssemblers)[i]->getQuadrature() == quad)) {
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
	}
      }
    }
 
    // check if all terms are pw_const
    bool pwConst = true;
    for (int i = 0; i < static_cast<int>( op->zeroOrder[myRank].size()); i++) {
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

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

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

  StandardZOA::StandardZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, false)
Thomas Witkowski's avatar
Thomas Witkowski committed
85
  {}
86
87
88
89
90
91
92

  void StandardZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

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

    int myRank = omp_get_thread_num();
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, nPoints, c);
    }
      
    if (symmetric) {
103
104
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
            
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

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

	for (int i = 0; i < nRow; i++) {
	  double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      //  non symmetric assembling 
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

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

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

  void StandardZOA::calculateElementMatrix(const ElInfo *rowElInfo,
					   const ElInfo *colElInfo,
144
145
					   const ElInfo *smallElInfo,
					   const ElInfo *largeElInfo,
146
147
					   ElementMatrix *mat) 
  {
148
149
150
151
    FUNCNAME("StandardZOA::calculateElementMatrix()");

    TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");

152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
    calculateElementMatrix(smallElInfo, mat);

    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();

    ElementMatrix tmpMat(nRow, nRow);

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nRow; j++) {
	tmpMat[i][j] = 0.0;

	for (int k = 0; k < nRow; k++) {
	  tmpMat[i][j] += m[i][j] * (*mat)[j][i];
	}
      }
    }

    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nRow; j++) {
	(*mat)[i][j] = tmpMat[i][j];
      }
    }

    /*
175
176
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
177
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
    int myRank = omp_get_thread_num();
180
    int nPoints = quadrature->getNumPoints();
181
    std::vector<double> c(nPoints, 0.0);
182

183
    std::vector<OperatorTerm*>::iterator termIt;
184
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt)
Thomas Witkowski's avatar
Thomas Witkowski committed
185
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(smallElInfo, nPoints, c);
186

Thomas Witkowski's avatar
Thomas Witkowski committed
187
188
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
189

Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
192
193
194
	for (int iq = 0; iq < nPoints; iq++) {
	  double val0 = c[iq] * smallElInfo->getDet() * quadrature->getWeight(iq) * 
	    (*(phi->getPhi(i)))(quadrature->getLambda(iq));

	  double val1 = 0.0;
195
196
	  for (int k = 0; k < nCol; k++)
	    val1 += m[j][k] * (*(psi->getPhi(k)))(quadrature->getLambda(iq));
197

Thomas Witkowski's avatar
Thomas Witkowski committed
198
199
200
	  val0 *= val1;

	  (*mat)[i][j] += val0;
201
202
	}
      }
203
    }
204
205
    */

206
207
208
209
210
  }

  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    int nPoints = quadrature->getNumPoints();
211
    std::vector<double> c(nPoints, 0.0);
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231

    int myRank = omp_get_thread_num();
    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();

      for (int i = 0; i < nRow; i++) {
	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
	  (quadrature->getLambda(iq));
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
      }
    }
  }

  FastQuadZOA::FastQuadZOA(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
232
  {}
233
234
235
236
237
238
239

  void FastQuadZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
240
241
242
243
244
245
246
247
248
249
#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_PHI);
	firstCall = false;
      }
250
251
    }

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

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

261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	c[iq] *= elInfo->getDet();

	const double *psi = psiFast->getPhi(iq);
	const double *phi = phiFast->getPhi(iq);
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
	  }
	}
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
292
293
294
295
  void FastQuadZOA::calculateElementMatrix(const ElInfo *rowElInfo,
					   const ElInfo *colElInfo,
					   const ElInfo *smallElInfo,
					   const ElInfo *largeElInfo,
					   ElementMatrix *mat) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
296
297
298
299
    FUNCNAME("FastQuadZOA::calculateElementMatrix()");

    ERROR_EXIT("CHECK FIRST, IF IT WILL REALLY WORK!\n");

Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
302
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
307
308
309
310
311
312
313
314
315
316

    if (firstCall) {
#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_PHI);
	firstCall = false;
      }
    }

317
    std::vector<double> c(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(rowElInfo, nPoints, c);
    }

    for (int iq = 0; iq < nPoints; iq++) {
      c[iq] *= smallElInfo->getDet();
      
      const double *psi = psiFast->getPhi(iq);
      const double *phi = phiFast->getPhi(iq);

      for (int i = 0; i < nCol; i++) { 
	for (int j = 0; j < nRow; j++) { 
	  double val = quadrature->getWeight(iq) * c[iq] * psi[i];

	  double tmpval = 0.0;
	  for (int k = 0; k < nCol; k++) {
335
	    tmpval += m[j][k] * phi[k];
Thomas Witkowski's avatar
Thomas Witkowski committed
336
337
338
339
340
341
342
343
344
345
	  }
	  val *= tmpval;

	  (*mat)[j][i] += val;
	}
      }

    }    
  }

346
347
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
348
349
350
    int myRank = omp_get_thread_num();
    int nPoints = quadrature->getNumPoints();

351
    if (firstCall) {
352
353
354
355
356
357
358
359
360
361
#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_PHI);
	firstCall = false;
      }
362
363
    }

364
    std::vector<double> c(nPoints, 0.0);
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
    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++) {
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
      }
    }
  }

  PrecalcZOA::PrecalcZOA(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

  void PrecalcZOA::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    if (firstCall) {
388
389
390
391
392
393
394
395
396
397
398
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);
	firstCall = false;
      }
399
400
    }

401
    std::vector<double> c(1, 0.0);
402
403
404
405
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

    for (int i = 0; i < size; i++) {
406
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, c);
407
408
    }

409
    c[0] *= elInfo->getDet();
410
411
412

    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
413
	(*mat)[i][i] += c[0] * q00->getValue(i,i);
414
	for (int j = i + 1; j < nCol; j++) {
415
	  double val = c[0] * q00->getValue(i, j);
416
417
418
419
420
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
    } else {
421
      for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
422
	for (int j = 0; j < nCol; j++) {
423
	  (*mat)[i][j] += c[0] * q00->getValue(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
424
	}
425
      }
426
427
428
429
430
431
    }
  }

  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    if (firstCall) {
432
433
434
435
436
437
438
439
440
441
442
#ifdef _OPENMP
#pragma omp critical
#endif 
      {
	q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					  owner->getColFESpace()->getBasisFcts(), 
					  quadrature);
	q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
				 quadrature);	
	firstCall = false;
      }
443
444
445
446
447
    }

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

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

452
    c[0] *= elInfo->getDet();
453
454

    for (int i = 0; i < nRow; i++)
455
      (*vec)[i] += c[0] * q0->getValue(i);
456
457
458
  }

}