Assembler.cc 12.3 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
      tmpMat(nRow, nCol),
27
28
29
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
30
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
31

32

Thomas Witkowski's avatar
Thomas Witkowski committed
33
  Assembler::~Assembler()
34
  {}
35

36

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

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

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

48
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
49
50
51
      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
80
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
81
82
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
83
					 ElementMatrix& userMat,
84
85
86
87
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

88
    if (remember && (factor != 1.0 || operat->uhOld))
89
90
      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

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

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

111
    if (secondOrderAssembler) {
112
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
113

114
      ElementMatrix &m = 
115
	smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
116
117
118
119
120
121
122
      
      if (smallElInfo == colElInfo)
	tmpMat = m * mat;	
      else
  	tmpMat = mat * trans(m);
      
      mat = tmpMat;
123
124
125
    }

    if (firstOrderAssemblerGrdPsi) {
126
127
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);

128
      if (largeElInfo == rowElInfo) {
129
130
	ElementMatrix &m = 
	  smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());
131

132
133
	tmpMat = m * mat;
      } else {
134
135
136
137
138
	ElementMatrix &m = 
	  smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
	
	tmpMat = mat * trans(m);
      }
139
140
	
      mat = tmpMat;
141
142
143
    }

    if (firstOrderAssemblerGrdPhi) {
144
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
145

146
      if (largeElInfo == colElInfo) {
147
148
149
150
151
152
153
154
155
156
157
158
	ElementMatrix &m = 
	  smallElInfo->getSubElemGradCoordsMat(rowFESpace->getBasisFcts()->getDegree());

	tmpMat = mat * trans(m);
      } else {
	ElementMatrix &m = 
	  smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
	
	tmpMat = m * mat;	
      }

      mat = tmpMat;
159
    }
160

161
162
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
Thomas Witkowski's avatar
Thomas Witkowski committed
163
      
164
165
      ElementMatrix &m = 
	smallElInfo->getSubElemCoordsMat(rowFESpace->getBasisFcts()->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
166
      
Thomas Witkowski's avatar
Thomas Witkowski committed
167
      if (smallElInfo == colElInfo)
168
169
	tmpMat = m * mat;	
      else
Thomas Witkowski's avatar
Thomas Witkowski committed
170
171
172
  	tmpMat = mat * trans(m);
      
      mat = tmpMat;
173
    }
174

175
176
    if (rememberElMat && &userMat != &elementMatrix)
      userMat += factor * elementMatrix;   
177
178
  }

179

180
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
181
					 ElementVector& userVec,
182
183
184
185
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

186
    if (remember && factor != 1.0)
187
188
189
190
      rememberElVec = true;

    Element *el = elInfo->getElement();

191
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
192
      initElement(elInfo);
193
    
Thomas Witkowski's avatar
Thomas Witkowski committed
194
    if (el != lastVecEl || !operat->isOptimized()) {
195
196
197
      if (rememberElVec)
	set_to_zero(elementVector);
	
198
199
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
200
      if (rememberElVec) {
201
	userVec += factor * elementVector;
202
203
204
	return;
      }
    }
205
206

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
207
    if (operat->uhOld && remember) {
208
      matVecAssemble(elInfo, vec);
209
      if (rememberElVec)
210
	userVec += factor * elementVector;      
211

212
213
      return;
    } 
214
215

    if (firstOrderAssemblerGrdPsi)
216
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
217
    if (zeroOrderAssembler)
218
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
219
      
220
    if (rememberElVec)
221
      userVec += factor * elementVector;    
222
223
  }

224

Thomas Witkowski's avatar
Thomas Witkowski committed
225
226
227
228
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
229
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
230
231
232
233
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

234
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
235
236
237
238
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

239
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
240
241
      initElement(smallElInfo, largeElInfo);
   
Thomas Witkowski's avatar
Thomas Witkowski committed
242
    if (el != lastVecEl || !operat->isOptimized()) {
243
244
245
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
246
247
248
      lastVecEl = el;
    } else {
      if (rememberElVec) {
249
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
250
251
252
	return;
      }
    }
253
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255

    if (operat->uhOld && remember) {
256
      if (smallElInfo->getLevel() == largeElInfo->getLevel())
Thomas Witkowski's avatar
Thomas Witkowski committed
257
	matVecAssemble(auxElInfo, vec);
258
259
      else
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);      
Thomas Witkowski's avatar
Thomas Witkowski committed
260

261
      if (rememberElVec)
262
	userVec += factor * elementVector;      
263

Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
266
267
268
      return;
    } 

    ERROR_EXIT("Not yet implemented!\n");

269
270
271
272
273
274
275
276
#if 0
    if (firstOrderAssemblerGrdPsi)
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);    
    if (zeroOrderAssembler)
      zeroOrderAssembler->calculateElementVector(elInfo, vec);    
    if (rememberElVec)
      axpy(factor, *elementVector, *userVec);
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
277
278
  }

279

280
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
281
282
283
  {
    FUNCNAME("Assembler::matVecAssemble()");

284
    Element *el = elInfo->getElement(); 
285
    double *uhOldLoc = new double[nRow];
286

287
    operat->uhOld->getLocalVector(el, uhOldLoc);
288
    
289
    if (el != lastMatEl) {
290
      set_to_zero(elementMatrix);
291
      calculateElementMatrix(elInfo, elementMatrix);
292
293
    }

294
    for (int i = 0; i < nRow; i++) {
295
      double val = 0.0;
296
      for (int j = 0; j < nRow; j++)
297
	val += elementMatrix[i][j] * uhOldLoc[j];
298
      
299
      vec[i] += val;
300
    }   
301

Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
    delete [] uhOldLoc;
  }

305

Thomas Witkowski's avatar
Thomas Witkowski committed
306
307
  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
308
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
309
310
311
312
313
314
315
316
317
318
319
320
321
322
  {
    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();
323
    std::vector<double> uhOldLoc(nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
324

325
    operat->uhOld->getLocalVector(auxEl, &(uhOldLoc[0]));
Thomas Witkowski's avatar
Thomas Witkowski committed
326
327

    if (mainEl != lastMatEl) {
328
      set_to_zero(elementMatrix);
329
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
330
 			     elementMatrix);    
Thomas Witkowski's avatar
Thomas Witkowski committed
331
    }
332

Thomas Witkowski's avatar
Thomas Witkowski committed
333
334
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
335
      for (int j = 0; j < nBasFcts; j++)
336
 	val += elementMatrix[i][j] * uhOldLoc[j];
337
      vec[i] += val;
338
    }   
339
340
  }

341

Thomas Witkowski's avatar
Thomas Witkowski committed
342
343
344
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
345
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
346
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
347
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
348
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
349
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
350
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
351
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
352
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
353
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
354
355
  }

356

357
  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
358
359
  { 
    if (secondOrderAssembler) {
360
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
361
362
      if (!secondOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
363
364
365
366
367
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
368
    if (firstOrderAssemblerGrdPsi) {
369
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
370
371
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
372
373
374
375
376
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
377
    if (firstOrderAssemblerGrdPhi) {
378
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
379
380
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
381
382
383
384
385
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
386
    if (zeroOrderAssembler) {
387
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
388
389
      if (!zeroOrderAssembler->getQuadrature()) {
	int dim = rowFESpace->getMesh()->getDim();
390
391
392
393
394
395
396
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

397

Thomas Witkowski's avatar
Thomas Witkowski committed
398
399
400
401
402
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
403

404

Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
410
411
412
					 const FiniteElemSpace *rowFeSpace,
					 const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
413
  {
414
    bool opt = (rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts());
Thomas Witkowski's avatar
Thomas Witkowski committed
415
416
417
418
419
420
421
422
423
424
425
426
427
428

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

429

Thomas Witkowski's avatar
Thomas Witkowski committed
430
431
432
433
434
  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
435
436
437
				       const FiniteElemSpace *rowFeSpace,
				       const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
  {
    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();
  }

454
}