Assembler.cc 13.2 KB
Newer Older
1
2
3
4
#include <vector>
#include <algorithm>
#include <boost/numeric/mtl/mtl.hpp>

5
6
7
8
9
10
11
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
12
#include "OpenMP.h"
13
14
15

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
16
17
18
  Assembler::Assembler(Operator *op,
		       const FiniteElemSpace *row,
		       const FiniteElemSpace *col) 
19
    : operat(op),
Thomas Witkowski's avatar
Thomas Witkowski committed
20
21
      rowFESpace(row),
      colFESpace(col ? col : row),
22
23
24
25
26
27
28
29
30
31
32
      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
33
34
35
36
37
38
39
40
41
42
43
44
45
  {
    elementMatrix = NEW ElementMatrix(nRow, nCol);
    elementVector = NEW ElementVector(nRow);
  }

  Assembler::~Assembler()
  {
    if (elementMatrix)
      DELETE elementMatrix;
      
    if (elementVector)
      DELETE elementVector;
  }
46
47
48
49
50
51
52

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

Thomas Witkowski's avatar
Thomas Witkowski committed
53
    if (remember && (factor != 1.0 || operat->uhOld)) {     
54
55
      rememberElMat = true;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
56

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

59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
    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
74
 
75
76
77
78
79
80
81
82
83
84
85
    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
86
    if (rememberElMat && userMat) {     
87
      axpy(factor, *elementMatrix, *userMat);
88
    }      
Thomas Witkowski's avatar
Thomas Witkowski committed
89

90
91
  }

92
93
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
94
95
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
96
97
98
99
100
101
102
103
104
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
105
    Element *el = smallElInfo->getElement();   
106
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
107
   
108
    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
    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, 
137
						 smallElInfo, largeElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
138
  
139
140
141
142
143
    if (rememberElMat && userMat) {
      axpy(factor, *elementMatrix, *userMat);
    }      
  }

144
145
146
147
148
149
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
					 ElementVector *userVec,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

Thomas Witkowski's avatar
Thomas Witkowski committed
150
    if (remember && factor != 1.0) {
151
152
153
154
155
      rememberElVec = true;
    }

    Element *el = elInfo->getElement();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
190
191
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
  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;
    }

    Element *el = mainElInfo->getElement();

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

  }

249
250
251
252
253
254
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

258
259
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
260
    if (el != lastMatEl) {
261
262
      calculateElementMatrix(elInfo, NULL);
    }
263
264
265
266
267
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
268
269
270
      }
      (*vec)[i] += val;
    }
271
272
    

Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
    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);

298
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
299
300
301
302

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
303
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
      }      
    }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
325
326
327
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
328
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
329
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
330
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
332
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
333
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
334
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
335
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
336
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
337
338
  }

339
340
341
  void Assembler::initElementMatrix(ElementMatrix *elMat, 
				    const ElInfo *rowElInfo,
				    const ElInfo *colElInfo)
342
  {
343
    TEST_EXIT_DBG(elMat)("No ElementMatrix!\n");
344
345

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
346
         
347
    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
348
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
349
						   &(elMat->rowIndices));
350

Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
353
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
354
355
356
357
358
359
360
361
362
363
      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
364
    }
365
366
  }

367
  void Assembler::initElementVector(ElementVector *elVec, const ElInfo *elInfo)
368
  {
369
    TEST_EXIT_DBG(elVec)("No ElementVector!\n");
370
371

    elVec->set(0.0);
372
373
374
375
 
    rowFESpace->getBasisFcts()->getLocalIndicesVec(elInfo->getElement(), 
						   rowFESpace->getAdmin(), 
						   &(elVec->dofIndices));
376
377
378
  }

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
  { 
    if (secondOrderAssembler) {
381
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
384
385
386
387
388
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
389
    if (firstOrderAssemblerGrdPsi) {
390
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
391
392
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
393
394
395
396
397
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
398
    if (firstOrderAssemblerGrdPhi) {
399
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
400
401
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
402
403
404
405
406
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
407
    if (zeroOrderAssembler) {
408
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
409
410
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
411
412
413
414
415
416
417
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
418
419
420
421
422
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471

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

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);

    checkQuadratures();
  }

  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *row,
				       const FiniteElemSpace *col) 
    : Assembler(op, row, col)
  {
    remember = false;

    // create sub assemblers
    secondOrderAssembler = 
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
    firstOrderAssemblerGrdPsi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
    firstOrderAssemblerGrdPhi = 
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
    zeroOrderAssembler = 
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);

    checkQuadratures();
  }

472
}