Assembler.cc 13.1 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
#include <vector>
#include <algorithm>

namespace AMDiS {

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

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

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

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

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

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

88
89
  }

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

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
103
    Element *el = smallElInfo->getElement();   
104
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
105
   
106
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
107
      initElement(smallElInfo, largeElInfo);
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
    }

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

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

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

    Element *el = elInfo->getElement();

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

Thomas Witkowski's avatar
Thomas Witkowski committed
188
189
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
  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);
//     }      

  }

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
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
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
    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;
321
322
  }

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

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
416
417
418
419
420
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
421
422
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

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

470
}