Assembler.cc 12 KB
Newer Older
1
2
3
#include <vector>
#include <algorithm>
#include <boost/numeric/mtl/mtl.hpp>
4
5
6
7
8
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "DOFVector.h"
9
#include "OpenMP.h"
10
11
12

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
13
14
15
  Assembler::Assembler(Operator *op,
		       const FiniteElemSpace *row,
		       const FiniteElemSpace *col) 
16
    : operat(op),
Thomas Witkowski's avatar
Thomas Witkowski committed
17
18
      rowFESpace(row),
      colFESpace(col ? col : row),
19
20
21
22
23
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
24
25
      elementMatrix(nRow, nCol),
      elementVector(nRow),
26
27
28
29
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
30
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
31
32

  Assembler::~Assembler()
33
  {}
34
35

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
36
					 ElementMatrix& userMat,
37
38
39
40
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

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

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

47
48
49
50
51
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(elInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
52
53
54
      if (rememberElMat)
	set_to_zero(elementMatrix);

55
56
57
      lastMatEl = el;
    } else {
      if (rememberElMat) {
58
	userMat += factor * elementMatrix;
59
60
61
	return;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
62
 
63
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
64
65
66
67
68
69
70
71
72
73

    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
74
    if (rememberElMat && &userMat != &elementMatrix)
75
      userMat += factor * elementMatrix;
76
77
  }

78
79
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
80
81
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
82
					 ElementMatrix& userMat,
83
84
85
86
87
88
89
90
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
Thomas Witkowski's avatar
Thomas Witkowski committed
91
    Element *el = smallElInfo->getElement();   
92
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
93
   
94
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
Thomas Witkowski's avatar
Thomas Witkowski committed
95
      initElement(smallElInfo, largeElInfo);
96
97
98
    }

    if (el != lastMatEl || !operat->isOptimized()) {
99
100
101
      if (rememberElMat)
	set_to_zero(elementMatrix);

102
103
104
      lastMatEl = el;
    } else {
      if (rememberElMat) {
105
	userMat += factor * elementMatrix;
106
107
108
109
	return;
      }
    }
  
110
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
111

112
    if (secondOrderAssembler) {
113
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
114
115
116
117
118
119
120
121
122
123

      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
      mat = tmpMat;
    }

    if (firstOrderAssemblerGrdPsi) {
      ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
124
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
125
126
127
128
    }

    if (firstOrderAssemblerGrdPhi) {
      ERROR_EXIT("Not yet implemented for not zero order assembler!\n");
129
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
130
    }
131

132
133
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
134

135
136
137
138
139
140
      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
      mat = tmpMat;
    }
141

142
143
    if (rememberElMat)
      userMat += factor * elementMatrix;
144
145
  }

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

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

    Element *el = elInfo->getElement();

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
    if (el != lastVecEl || !operat->isOptimized()) {
163
164
165
      if (rememberElVec)
	set_to_zero(elementVector);
	
166
167
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
168
      if (rememberElVec) {
169
	userVec += factor * elementVector;
170
171
172
	return;
      }
    }
173
174

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
175
    if (operat->uhOld && remember) {
176
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
177
      if (rememberElVec) {
178
	userVec += factor * elementVector;
179
180
181
      }
      return;
    } 
182
183

    if (firstOrderAssemblerGrdPsi)
184
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
185
    if (zeroOrderAssembler)
186
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
187
188
    if (rememberElVec)
      userVec += factor * elementVector;
189
190
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
193
194
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
195
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
					 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()) {
211
212
213
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
214
215
216
      lastVecEl = el;
    } else {
      if (rememberElVec) {
217
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
220
	return;
      }
    }
221
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
222
223
224
225
226
227
228
229
230
231

    if (operat->uhOld && remember) {

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

      if (rememberElVec) {
232
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
233
      }
234

Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
      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);
//     }      

  }

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

    Element *el = elInfo->getElement(); 
256
    double *uhOldLoc = new double[nRow];
257

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

Thomas Witkowski's avatar
Thomas Witkowski committed
271
272
273
274
275
    delete [] uhOldLoc;
  }

  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
276
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
  {
    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);

296
297
    mtl::dense2D<double> m(nRow, nRow);
    smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
298
299
300
301

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
302
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
306
307
      }      
    }

    
    if (mainEl != lastMatEl) {
308
309
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
			     elementMatrix);
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
312
313
314
    }
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
315
	val += elementMatrix[i][j] * uhOldLoc2[j];
Thomas Witkowski's avatar
Thomas Witkowski committed
316
      }
317
      vec[i] += val;
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
320
321
    }
    
    delete [] uhOldLoc;
    delete [] uhOldLoc2;
322
323
  }

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

  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
339
340
  { 
    if (secondOrderAssembler) {
341
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
344
345
346
347
348
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
349
    if (firstOrderAssemblerGrdPsi) {
350
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
351
352
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
353
354
355
356
357
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
358
    if (firstOrderAssemblerGrdPhi) {
359
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
360
361
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
362
363
364
365
366
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
367
    if (zeroOrderAssembler) {
368
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
369
370
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
371
372
373
374
375
376
377
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
378
379
380
381
382
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431

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

432
}