Assembler.cc 12.7 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),
17
18
19
20
      rowFeSpace(row),
      colFeSpace(col ? col : row),
      nRow(rowFeSpace->getBasisFcts()->getNumber()),
      nCol(colFeSpace->getBasisFcts()->getNumber()),
21
22
23
      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
					 bool rowColFeSpaceEqual,
84
					 ElementMatrix& userMat,
85
86
87
88
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

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

    if (el != lastMatEl || !operat->isOptimized()) {
103
104
105
      if (rememberElMat)
	set_to_zero(elementMatrix);

106
107
108
      lastMatEl = el;
    } else {
      if (rememberElMat) {
109
	userMat += factor * elementMatrix;
110
111
112
	return;
      }
    }
113
 
114
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
115

116
    if (secondOrderAssembler) {
117
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
118

119
      ElementMatrix &m = 
120
	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
121
      
122
123
124
125
126
127
128
129
      if (!rowColFeSpaceEqual) {
	if (smallElInfo == colElInfo)
	  tmpMat = m * mat;	
	else
	  tmpMat = mat * trans(m);
	
	mat = tmpMat;
      }
130
131
132
    }

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

135
136
137
138
139
140
141
142
143
144
145
146
      if (!rowColFeSpaceEqual) {
	if (largeElInfo == rowElInfo) {
	  ElementMatrix &m = 
	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	  
	  tmpMat = m * mat;
	} else {
	  ElementMatrix &m = 
	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	  
	  tmpMat = mat * trans(m);
	}
147
	
148
	mat = tmpMat;
149
      }
150
151
152
    }

    if (firstOrderAssemblerGrdPhi) {
153
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
154

155
156
157
158
159
160
161
162
163
164
165
166
      if (!rowColFeSpaceEqual) {
	if (largeElInfo == colElInfo) {
	  ElementMatrix &m = 
	    smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	  
	  tmpMat = mat * trans(m);
	} else {
	  ElementMatrix &m = 
	    smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	  
	  tmpMat = m * mat;	
	}
167
	
168
	mat = tmpMat;
169
      }
170
    }
171

172
173
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
174

175
176
177
178
179
180
181
182
183
184
185
      if (!rowColFeSpaceEqual) {
	ElementMatrix &m = 
	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	
	if (smallElInfo == colElInfo)
	  tmpMat = m * mat;
	else 
	  tmpMat = mat * trans(m);
	
	mat = tmpMat;
      }
186
    }
187

188
    if (rememberElMat && &userMat != &elementMatrix)
189
      userMat += factor * elementMatrix;       
190
191
  }

192

193
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
194
					 ElementVector& userVec,
195
196
197
198
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

199
    if (remember && factor != 1.0)
200
201
202
203
      rememberElVec = true;

    Element *el = elInfo->getElement();

204
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
205
      initElement(elInfo);
206
    
Thomas Witkowski's avatar
Thomas Witkowski committed
207
    if (el != lastVecEl || !operat->isOptimized()) {
208
209
210
      if (rememberElVec)
	set_to_zero(elementVector);
	
211
212
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
213
      if (rememberElVec) {
214
	userVec += factor * elementVector;
215
216
217
	return;
      }
    }
218
219

    ElementVector& vec = rememberElVec ? elementVector : userVec;
220

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    if (operat->uhOld && remember) {
222
      matVecAssemble(elInfo, vec);
223
      if (rememberElVec)
224
	userVec += factor * elementVector;      
225

226
227
      return;
    } 
228
229

    if (firstOrderAssemblerGrdPsi)
230
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
231
    if (zeroOrderAssembler)
232
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
233
      
234
    if (rememberElVec)
235
      userVec += factor * elementVector;    
236
237
  }

238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241
242
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
243
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
246
247
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

248
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
251
252
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

253
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
254
255
      initElement(smallElInfo, largeElInfo);
   
Thomas Witkowski's avatar
Thomas Witkowski committed
256
    if (el != lastVecEl || !operat->isOptimized()) {
257
258
259
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
260
261
262
      lastVecEl = el;
    } else {
      if (rememberElVec) {
263
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
266
	return;
      }
    }
267
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
268
269

    if (operat->uhOld && remember) {
270
      if (smallElInfo->getLevel() == largeElInfo->getLevel())
Thomas Witkowski's avatar
Thomas Witkowski committed
271
	matVecAssemble(auxElInfo, vec);
272
273
      else
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);      
Thomas Witkowski's avatar
Thomas Witkowski committed
274

275
      if (rememberElVec)
276
	userVec += factor * elementVector;      
277

Thomas Witkowski's avatar
Thomas Witkowski committed
278
279
280
      return;
    } 

281
282
283
284
285
286
287
288
289
290
    if (firstOrderAssemblerGrdPsi) {
      ERROR_EXIT("Not yet implemented!\n");
    }

    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementVector(smallElInfo, vec);
      
      if (smallElInfo != mainElInfo) {
	ElementVector tmpVec(vec);	
	ElementMatrix &m = 
291
	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
292
293
294
295
296

	tmpVec = m * vec;	
	vec = tmpVec;
      }      
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
297

298
    if (rememberElVec)
299
      userVec += factor * elementVector;    
Thomas Witkowski's avatar
Thomas Witkowski committed
300
301
  }

302

303
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
304
305
306
  {
    FUNCNAME("Assembler::matVecAssemble()");

307
    Element *el = elInfo->getElement(); 
308
    ElementVector uhOldLoc(operat->uhOld->getFeSpace() == rowFeSpace ? nRow : nCol);
309

310
    operat->uhOld->getLocalVector(el, uhOldLoc);
311
    
312
    if (el != lastMatEl) {
313
      set_to_zero(elementMatrix);
314
      calculateElementMatrix(elInfo, elementMatrix);
315
316
    }

317
    for (int i = 0; i < nRow; i++) {
318
      double val = 0.0;
319
      for (int j = 0; j < nCol; j++)
320
	val += elementMatrix[i][j] * uhOldLoc[j];
321
      
322
      vec[i] += val;
323
    }   
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
  }

326

Thomas Witkowski's avatar
Thomas Witkowski committed
327
328
  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
329
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
330
331
332
  {
    FUNCNAME("Assembler::matVecAssemble()");

333
    TEST_EXIT(rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts())
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
      ("Works only for equal basis functions for different components!\n");

336
    TEST_EXIT(operat->uhOld->getFeSpace()->getMesh() == auxElInfo->getMesh())
Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
339
340
341
      ("Da stimmt was nicht!\n");

    Element *mainEl = mainElInfo->getElement(); 
    Element *auxEl = auxElInfo->getElement();

342
    const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    int nBasFcts = basFcts->getNumber();
344
    ElementVector uhOldLoc(nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
345

346
    operat->uhOld->getLocalVector(auxEl, uhOldLoc);
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348

    if (mainEl != lastMatEl) {
349
      set_to_zero(elementMatrix);
350
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
351
			     false, elementMatrix);    
Thomas Witkowski's avatar
Thomas Witkowski committed
352
    }
353

Thomas Witkowski's avatar
Thomas Witkowski committed
354
355
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
356
      for (int j = 0; j < nBasFcts; j++)
357
 	val += elementMatrix[i][j] * uhOldLoc[j];
358
      vec[i] += val;
359
    }   
360
361
  }

362

Thomas Witkowski's avatar
Thomas Witkowski committed
363
364
365
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
366
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
367
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
368
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
369
    if (firstOrderAssemblerGrdPsi)
Thomas Witkowski's avatar
Thomas Witkowski committed
370
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
371
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
372
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
373
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
374
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
375
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
      if (!secondOrderAssembler->getQuadrature()) {
383
	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
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
392
	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
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
401
	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
      if (!zeroOrderAssembler->getQuadrature()) {
410
	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);
      }
    }
  }

418

Thomas Witkowski's avatar
Thomas Witkowski committed
419
420
421
422
423
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
424

425

Thomas Witkowski's avatar
Thomas Witkowski committed
426
427
428
429
430
  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
431
432
433
					 const FiniteElemSpace *rowFeSpace,
					 const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
434
  {
435
    bool opt = (rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts());
Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
440
441
442
443
444
445
446
447
448
449

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

450

Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
453
454
455
  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
456
457
458
				       const FiniteElemSpace *rowFeSpace,
				       const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
  {
    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();
  }

475
}