Assembler.cc 14 KB
Newer Older
1
2
3
4
5
6
7
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
8
#include "OpenMP.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <vector>
#include <algorithm>

namespace AMDiS {

  Assembler::Assembler(Operator  *op,
		       const FiniteElemSpace *rowFESpace_,
		       const FiniteElemSpace *colFESpace_) 
    : operat(op),
      rowFESpace(rowFESpace_),
      colFESpace(colFESpace_ ? colFESpace_ : rowFESpace_),
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
      elementMatrix(NULL),
      elementVector(NULL),
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
Thomas Witkowski's avatar
Thomas Witkowski committed
31
  {}
32
33
34
35
36
37
38

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

Thomas Witkowski's avatar
Thomas Witkowski committed
39
    if (remember && ((factor != 1.0) || (operat->uhOld))) {      
40
41
      rememberElMat = true;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
42

43
44
45
46
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
47
   
48
49
    checkForNewTraverse();
    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
50

51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(elInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
66
 
67
68
69
70
71
72
73
74
75
76
77
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

    if (secondOrderAssembler)
      secondOrderAssembler->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementMatrix(elInfo, mat);
    if (firstOrderAssemblerGrdPhi)
      firstOrderAssemblerGrdPhi->calculateElementMatrix(elInfo, mat);
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementMatrix(elInfo, mat);

Thomas Witkowski's avatar
Thomas Witkowski committed
78
    if (rememberElMat && userMat) {     
79
      axpy(factor, *elementMatrix, *userMat);
80
    }      
Thomas Witkowski's avatar
Thomas Witkowski committed
81

82
83
  }

84
85
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
86
87
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
88
89
90
91
92
93
94
95
96
97
98
99
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

Thomas Witkowski's avatar
Thomas Witkowski committed
100
    Element *el = smallElInfo->getElement();
101
102
103
104
105
106
107
108

    
    //    checkForNewTraverse();
    lastVecEl = lastMatEl = NULL;

    checkQuadratures();
    
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
      initElement(smallElInfo, largeElInfo);
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

Thomas Witkowski's avatar
Thomas Witkowski committed
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
    if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
      if (secondOrderAssembler)
	secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
      if (firstOrderAssemblerGrdPsi)
	firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
      if (firstOrderAssemblerGrdPhi)
	firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
      if (zeroOrderAssembler)
	zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
    } else {
      if (secondOrderAssembler)
	secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
						     smallElInfo, largeElInfo, mat);   
      if (firstOrderAssemblerGrdPsi)
	firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo, 
							  smallElInfo, largeElInfo, mat);    
      if (firstOrderAssemblerGrdPhi)
	firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo, 
							  smallElInfo, largeElInfo, mat);    

      if (zeroOrderAssembler) 
	zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
148
						 smallElInfo, largeElInfo, mat);
149
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
150

151
152
153
154
155
    if (rememberElMat && userMat) {
      axpy(factor, *elementMatrix, *userMat);
    }      
  }

156
157
158
159
160
161
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
					 ElementVector *userVec,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
162
    if (remember && factor != 1.0) {
163
164
165
      rememberElVec = true;
    }

166
    if (rememberElVec && !elementVector) {
167
      elementVector = NEW ElementVector(nRow);
168
    }
169
170
171
172
173

    Element *el = elInfo->getElement();
    checkForNewTraverse();
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
174
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
175
176
      initElement(elInfo);
    }
177
    
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
180
181
182
183
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
184
      if (rememberElVec) {
185
186
187
188
189
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
190
    if (operat->uhOld && remember) {
191
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
192
      if (rememberElVec) {
193
194
195
196
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
197
    if (firstOrderAssemblerGrdPsi) {
198
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
199
200
    }
    if (zeroOrderAssembler) {
201
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
202
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
203
    if (rememberElVec) {
204
      axpy(factor, *elementVector, *userVec);
205
    }      
206
207
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
208
209
210
211
212
213
214
215
216
217
218
219
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
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
					 ElementVector *userVec, 
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

    if (remember && factor != 1.0) {
      rememberElVec = true;
    }

    if (rememberElVec && !elementVector) {
      elementVector = NEW ElementVector(nRow);
    }

    Element *el = mainElInfo->getElement();
    //    checkForNewTraverse();
    checkQuadratures();

    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(auxElInfo);
    }
    
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
      if (rememberElVec) {
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;

    if (operat->uhOld && remember) {

      if (smallElInfo->getLevel() == largeElInfo->getLevel()) {
	matVecAssemble(auxElInfo, vec);
      } else {
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);
      }

      if (rememberElVec) {
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 

    ERROR_EXIT("Not yet implemented!\n");
//     if (firstOrderAssemblerGrdPsi) {
//       firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
//     }
//     if (zeroOrderAssembler) {
//       zeroOrderAssembler->calculateElementVector(elInfo, vec);
//     }
//     if (rememberElVec) {
//       axpy(factor, *elementVector, *userVec);
//     }      

  }

273
274
275
276
277
278
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
279
280
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
281

282
283
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
284
    if (el != lastMatEl) {
285
286
      calculateElementMatrix(elInfo, NULL);
    }
287
288
289
290
291
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
292
293
294
      }
      (*vec)[i] += val;
    }
295
296
    

Thomas Witkowski's avatar
Thomas Witkowski committed
297
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
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
    delete [] uhOldLoc;
  }

  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
				 ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

    TEST_EXIT(rowFESpace->getBasisFcts() == colFESpace->getBasisFcts())
      ("Works only for equal basis functions for different components!\n");

    TEST_EXIT(operat->uhOld->getFESpace()->getMesh() == auxElInfo->getMesh())
      ("Da stimmt was nicht!\n");

    Element *mainEl = mainElInfo->getElement(); 
    Element *auxEl = auxElInfo->getElement();

    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
    double *uhOldLoc2 = new double[nBasFcts];

    operat->uhOld->getLocalVector(auxEl, uhOldLoc);

    DimMat<double> *m = smallElInfo->getSubElemCoordsMat();

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	uhOldLoc2[i] += (*m)[j][i] * uhOldLoc[i];
      }      
    }

    
    if (mainEl != lastMatEl) {
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, NULL);
    }
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc2[j];
      }
      (*vec)[i] += val;
    }
    

    delete [] uhOldLoc;
    delete [] uhOldLoc2;
347
348
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
349
350
351
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
352
353
354
  {
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
355
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
356
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
357
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
358
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
359
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
360
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
361
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
362
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
363
364
365
366
367
368
369
370
371
372
373
374
375
376
  }

  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
					 const FiniteElemSpace *rowFESpace_,
					 const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    bool opt = (rowFESpace_ == colFESpace_);

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
377
    secondOrderAssembler = 
378
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
379
    firstOrderAssemblerGrdPsi = 
380
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    firstOrderAssemblerGrdPhi = 
382
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
383
    zeroOrderAssembler = 
384
385
386
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
  }

387
  StandardAssembler::StandardAssembler(Operator *op,
388
389
390
391
392
393
394
395
396
397
398
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *rowFESpace_,
				       const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    remember = false;

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
399
    secondOrderAssembler = 
400
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
401
    firstOrderAssemblerGrdPsi = 
402
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
403
    firstOrderAssemblerGrdPhi = 
404
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
405
    zeroOrderAssembler = 
406
407
408
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

409
410
411
  void Assembler::initElementMatrix(ElementMatrix *elMat, 
				    const ElInfo *rowElInfo,
				    const ElInfo *colElInfo)
412
  {
413
    TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
414
415

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
416
         
417
    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
418
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
419
						   &(elMat->rowIndices));
420

Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
423
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
424
425
426
427
428
429
430
431
432
433
      if (colElInfo) {
	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
434
    }
435
436
  }

437
  void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
438
  {
439
    TEST_EXIT_DBG(elVec)("No ElementVector!\n");
440
441

    elVec->set(0.0);
442
443
444
445
 
    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						   rowFESpace->getAdmin(), 
						   &(elVec->dofIndices));
446
447
448
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
  { 
    if (secondOrderAssembler) {
451
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
452
453
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
454
455
456
457
458
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
459
    if (firstOrderAssemblerGrdPsi) {
460
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
463
464
465
466
467
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
468
    if (firstOrderAssemblerGrdPhi) {
469
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
470
471
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
472
473
474
475
476
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
477
    if (zeroOrderAssembler) {
478
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
479
480
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
481
482
483
484
485
486
487
488
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

}