Assembler.cc 13.5 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
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
45
46
  
    checkQuadratures();
47
48

    Element *el = elInfo->getElement();
Thomas Witkowski's avatar
Thomas Witkowski committed
49

50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
    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
65
 
66
67
68
69
70
71
72
73
74
75
76
    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
77
    if (rememberElMat && userMat) {     
78
      axpy(factor, *elementMatrix, *userMat);
79
    }      
Thomas Witkowski's avatar
Thomas Witkowski committed
80

81
82
  }

83
84
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
85
86
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
87
88
89
90
91
92
93
94
95
96
97
98
					 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
99
    Element *el = smallElInfo->getElement();
100
101
102
103
104
105
    
    lastVecEl = lastMatEl = NULL;

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

    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
123
124
125
126
127
128
129
130
131
132
133
    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, 
134
						 smallElInfo, largeElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
135
  
136
137
138
139
140
    if (rememberElMat && userMat) {
      axpy(factor, *elementMatrix, *userMat);
    }      
  }

141
142
143
144
145
146
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
					 ElementVector *userVec,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
147
    if (remember && factor != 1.0) {
148
149
150
      rememberElVec = true;
    }

151
    if (rememberElVec && !elementVector) {
152
      elementVector = NEW ElementVector(nRow);
153
    }
154
155
156
157

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

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
159
160
      initElement(elInfo);
    }
161
    
Thomas Witkowski's avatar
Thomas Witkowski committed
162
163
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
164
165
166
167
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
      if (rememberElVec) {
169
170
171
172
173
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
174
    if (operat->uhOld && remember) {
175
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
176
      if (rememberElVec) {
177
178
179
180
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
181
    if (firstOrderAssemblerGrdPsi) {
182
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
183
184
    }
    if (zeroOrderAssembler) {
185
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
186
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
187
    if (rememberElVec) {
188
      axpy(factor, *elementVector, *userVec);
189
    }      
190
191
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
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
  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();
    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);
//     }      

  }

256
257
258
259
260
261
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

265
266
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
267
    if (el != lastMatEl) {
268
269
      calculateElementMatrix(elInfo, NULL);
    }
270
271
272
273
274
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
275
276
277
      }
      (*vec)[i] += val;
    }
278
279
    

Thomas Witkowski's avatar
Thomas Witkowski committed
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
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
    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;
330
331
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
332
333
334
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
335
336
337
  {
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
338
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
339
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
340
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
341
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
342
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
343
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
344
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
345
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
346
347
348
349
350
351
352
353
354
355
356
357
358
359
  }

  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
360
    secondOrderAssembler = 
361
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    firstOrderAssemblerGrdPsi = 
363
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
364
    firstOrderAssemblerGrdPhi = 
365
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
366
    zeroOrderAssembler = 
367
368
369
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
  }

370
  StandardAssembler::StandardAssembler(Operator *op,
371
372
373
374
375
376
377
378
379
380
381
				       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
382
    secondOrderAssembler = 
383
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
384
    firstOrderAssemblerGrdPsi = 
385
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
386
    firstOrderAssemblerGrdPhi = 
387
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
388
    zeroOrderAssembler = 
389
390
391
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

392
393
394
  void Assembler::initElementMatrix(ElementMatrix *elMat, 
				    const ElInfo *rowElInfo,
				    const ElInfo *colElInfo)
395
  {
396
    TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
397
398

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
399
         
400
    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
401
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
402
						   &(elMat->rowIndices));
403

Thomas Witkowski's avatar
Thomas Witkowski committed
404
405
406
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
407
408
409
410
411
412
413
414
415
416
      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
417
    }
418
419
  }

420
  void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
421
  {
422
    TEST_EXIT_DBG(elVec)("No ElementVector!\n");
423
424

    elVec->set(0.0);
425
426
427
428
 
    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						   rowFESpace->getAdmin(), 
						   &(elVec->dofIndices));
429
430
431
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
  { 
    if (secondOrderAssembler) {
434
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
435
436
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
437
438
439
440
441
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
442
    if (firstOrderAssemblerGrdPsi) {
443
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
444
445
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
446
447
448
449
450
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
451
    if (firstOrderAssemblerGrdPhi) {
452
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
453
454
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
455
456
457
458
459
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
460
    if (zeroOrderAssembler) {
461
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
462
463
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
464
465
466
467
468
469
470
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
471
472
473
474
475
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
476
}