ZeroOrderAssembler.cc 12.4 KB
Newer Older
1
2
3
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
32
33
34
35
36
#include <vector>
#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
    if (op->zeroOrder.size() == 0) {
      return NULL;
    }

    ZeroOrderAssembler *newAssembler;

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

    int myRank = omp_get_thread_num();
40
    std::vector<OperatorTerm*> opTerms = op->zeroOrder[myRank];
41
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
152
153
154
155
    FUNCNAME("StandardZOA::calculateElementMatrix()");

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

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
    DimMat<double> *m = smallElInfo->getSubElemCoordsMat();

Thomas Witkowski's avatar
Thomas Witkowski committed
156
    TEST_EXIT(m)("No subElemCoordsMat!\n");
157

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    int myRank = omp_get_thread_num();
159
    int nPoints = quadrature->getNumPoints();
160
    std::vector<double> c(nPoints, 0.0);
161

162
163
    std::vector<OperatorTerm*>::iterator termIt;
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
Thomas Witkowski's avatar
Thomas Witkowski committed
164
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(smallElInfo, nPoints, c);
165
166
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
167
168
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
169

Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
173
174
	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;
175
	  for (int k = 0; k < nCol; k++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
176
	    val1 += (*m)[j][k] * (*(psi->getPhi(k)))(quadrature->getLambda(iq));
177
	  }
178

Thomas Witkowski's avatar
Thomas Witkowski committed
179
180
181
	  val0 *= val1;

	  (*mat)[i][j] += val0;
182
	}
Thomas Witkowski's avatar
Thomas Witkowski committed
183

184
      }
185
186
187
188
189
190
    }
  }

  void StandardZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    int nPoints = quadrature->getNumPoints();
191
    std::vector<double> c(nPoints, 0.0);
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211

    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)
212
  {}
213
214
215
216
217
218
219

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

    if (firstCall) {
220
221
222
223
224
225
226
227
228
229
#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;
      }
230
231
    }

232
    std::vector<double> c(nPoints, 0.0);
233
234
235
236
237
238
    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) {
239
240
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
      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
270
271
272
273
274
275
  void FastQuadZOA::calculateElementMatrix(const ElInfo *rowElInfo,
					   const ElInfo *colElInfo,
					   const ElInfo *smallElInfo,
					   const ElInfo *largeElInfo,
					   ElementMatrix *mat) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
278
279
    FUNCNAME("FastQuadZOA::calculateElementMatrix()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
    DimMat<double> *m = smallElInfo->getSubElemCoordsMat();

    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;
      }
    }

297
    std::vector<double> c(nPoints, 0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
    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++) {
	    tmpval += (*m)[j][k] * phi[k];
	  }
	  val *= tmpval;

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

    }    
  }

326
327
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
328
329
330
    int myRank = omp_get_thread_num();
    int nPoints = quadrature->getNumPoints();

331
    if (firstCall) {
332
333
334
335
336
337
338
339
340
341
#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;
      }
342
343
    }

344
    std::vector<double> c(nPoints, 0.0);
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
    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) {
368
369
370
371
372
373
374
375
376
377
378
#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;
      }
379
380
    }

381
    std::vector<double> c(1, 0.0);
382
383
384
385
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

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

389
    c[0] *= elInfo->getDet();
390
391
392

    if (symmetric) {
      for (int i = 0; i < nRow; i++) {
393
	(*mat)[i][i] += c[0] * q00->getValue(i,i);
394
	for (int j = i + 1; j < nCol; j++) {
395
	  double val = c[0] * q00->getValue(i, j);
396
397
398
399
400
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
    } else {
401
      for (int i = 0; i < nRow; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
402
	for (int j = 0; j < nCol; j++) {
403
	  (*mat)[i][j] += c[0] * q00->getValue(i, j);
Thomas Witkowski's avatar
Thomas Witkowski committed
404
	}
405
      }
406
407
408
409
410
411
    }
  }

  void PrecalcZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    if (firstCall) {
412
413
414
415
416
417
418
419
420
421
422
#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;
      }
423
424
425
426
427
    }

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

    int myRank = omp_get_thread_num();
428
    std::vector<double> c(1, 0.0);
429
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
430
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
431
432
    }

433
    c[0] *= elInfo->getDet();
434
435

    for (int i = 0; i < nRow; i++)
436
      (*vec)[i] += c[0] * q0->getValue(i);
437
438
439
  }

}