ZeroOrderAssembler.cc 12.4 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
69

    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
70
71
72
73
74
75
76
77
78
    if (!optimized) {
	newAssembler = NEW StandardZOA(op, assembler, quad);
    } else {
      if (pwConst) {
	  newAssembler = NEW PrecalcZOA(op, assembler, quad);
      } else {
	  newAssembler = NEW FastQuadZOA(op, assembler, quad);
      }
    }
79
80
81
82
83
84
85

    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
86
  {}
87
88
89
90
91
92
93

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

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

    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) {
104
105
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");
            
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
144
      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,
145
146
					   const ElInfo *smallElInfo,
					   const ElInfo *largeElInfo,
147
148
					   ElementMatrix *mat) 
  {
149
150
151
152
153
154
    FUNCNAME("StandardZOA::calculateElementMatrix()");

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

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
155
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
156

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
176
177
178
	  val0 *= val1;

	  (*mat)[i][j] += val0;
179
180
	}
      }
181
182
183
184
185
186
    }
  }

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

    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)
208
  {}
209
210
211
212
213
214
215

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

    if (firstCall) {
216
217
218
219
220
221
222
223
224
225
#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;
      }
226
227
    }

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

237
238
239
240
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
      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
266
267
268
269
270
271
  void FastQuadZOA::calculateElementMatrix(const ElInfo *rowElInfo,
					   const ElInfo *colElInfo,
					   const ElInfo *smallElInfo,
					   const ElInfo *largeElInfo,
					   ElementMatrix *mat) 
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
272
273
274
275
    FUNCNAME("FastQuadZOA::calculateElementMatrix()");

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

Thomas Witkowski's avatar
Thomas Witkowski committed
276
277
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();
278
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
281
282
283
284
285
286
287
288
289
290
291
292

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

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

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

    }    
  }

322
323
  void FastQuadZOA::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
324
325
326
    int myRank = omp_get_thread_num();
    int nPoints = quadrature->getNumPoints();

327
    if (firstCall) {
328
329
330
331
332
333
334
335
336
337
#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;
      }
338
339
    }

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

377
    std::vector<double> c(1, 0.0);
378
379
380
381
    int myRank = omp_get_thread_num();
    int size = static_cast<int>(terms[myRank].size());

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

385
    c[0] *= elInfo->getDet();
386
387
388

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

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

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

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

429
    c[0] *= elInfo->getDet();
430
431

    for (int i = 0; i < nRow; i++)
432
      (*vec)[i] += c[0] * q0->getValue(i);
433
434
435
  }

}