Assembler.cc 12.9 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
#include <vector>
#include <algorithm>
#include <boost/numeric/mtl/mtl.hpp>
16
17
18
19
20
21
22
23
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "DOFVector.h"

namespace AMDiS {

Thomas Witkowski's avatar
Thomas Witkowski committed
24
25
26
  Assembler::Assembler(Operator *op,
		       const FiniteElemSpace *row,
		       const FiniteElemSpace *col) 
27
    : operat(op),
28
29
30
31
      rowFeSpace(row),
      colFeSpace(col ? col : row),
      nRow(rowFeSpace->getBasisFcts()->getNumber()),
      nCol(colFeSpace->getBasisFcts()->getNumber()),
32
33
34
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
35
36
      elementMatrix(nRow, nCol),
      elementVector(nRow),
37
      tmpMat(nRow, nCol),
38
39
40
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
41
  {}
Thomas Witkowski's avatar
Thomas Witkowski committed
42

43

Thomas Witkowski's avatar
Thomas Witkowski committed
44
  Assembler::~Assembler()
45
  {}
46

47

48
  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
49
					 ElementMatrix& userMat,
50
51
52
53
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

54
    if (remember && (factor != 1.0 || operat->uhOld))
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
    if (el != lastMatEl || !operat->isOptimized()) {
60
61
      initElement(elInfo);

62
63
64
      if (rememberElMat)
	set_to_zero(elementMatrix);

65
66
67
      lastMatEl = el;
    } else {
      if (rememberElMat) {
68
	userMat += factor * elementMatrix;
69
70
71
	return;
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
72
 
73
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
74
75
76
77
78
79
80
81
82
83

    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 != &elementMatrix)
85
      userMat += factor * elementMatrix;
86
87
  }

88

89
90
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
91
92
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
93
					 bool rowColFeSpaceEqual,
94
					 ElementMatrix& userMat,
95
96
97
98
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

99
    if (remember && (factor != 1.0 || operat->uhOld))
100
101
      rememberElMat = true;
  
Thomas Witkowski's avatar
Thomas Witkowski committed
102
    Element *el = smallElInfo->getElement();   
103
    lastVecEl = lastMatEl = NULL;
Thomas Witkowski's avatar
Thomas Witkowski committed
104
   
105
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
106
      if (smallElInfo == colElInfo)
107
	initElement(smallElInfo);	
108
      else
109
	initElement(smallElInfo, largeElInfo);      
110
    }      
111
112

    if (el != lastMatEl || !operat->isOptimized()) {
113
114
115
      if (rememberElMat)
	set_to_zero(elementMatrix);

116
117
118
      lastMatEl = el;
    } else {
      if (rememberElMat) {
119
	userMat += factor * elementMatrix;
120
121
122
	return;
      }
    }
123
 
124
    ElementMatrix& mat = rememberElMat ? elementMatrix : userMat;
125

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

129
      ElementMatrix &m = 
130
	smallElInfo->getSubElemGradCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
131
      
132
133
134
135
136
137
138
139
      if (!rowColFeSpaceEqual) {
	if (smallElInfo == colElInfo)
	  tmpMat = m * mat;	
	else
	  tmpMat = mat * trans(m);
	
	mat = tmpMat;
      }
140
141
142
    }

    if (firstOrderAssemblerGrdPsi) {
143
144
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);

145
146
147
148
149
150
151
152
153
154
155
156
      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);
	}
157
	
158
	mat = tmpMat;
159
      }
160
161
162
    }

    if (firstOrderAssemblerGrdPhi) {
163
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
164

165
166
167
168
169
170
171
172
173
174
175
176
      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;	
	}
177
	
178
	mat = tmpMat;
179
      }
180
    }
181

182
183
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
184

185
186
187
188
189
190
191
192
193
194
195
      if (!rowColFeSpaceEqual) {
	ElementMatrix &m = 
	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
	
	if (smallElInfo == colElInfo)
	  tmpMat = m * mat;
	else 
	  tmpMat = mat * trans(m);
	
	mat = tmpMat;
      }
196
    }
197

198
    if (rememberElMat && &userMat != &elementMatrix)
199
      userMat += factor * elementMatrix;       
200
201
  }

202

203
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
204
					 ElementVector& userVec,
205
206
207
208
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

209
    if (remember && factor != 1.0)
210
211
212
213
      rememberElVec = true;

    Element *el = elInfo->getElement();

214
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
215
      initElement(elInfo);
216
    
Thomas Witkowski's avatar
Thomas Witkowski committed
217
    if (el != lastVecEl || !operat->isOptimized()) {
218
219
220
      if (rememberElVec)
	set_to_zero(elementVector);
	
221
222
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
223
      if (rememberElVec) {
224
	userVec += factor * elementVector;
225
226
227
	return;
      }
    }
228
229

    ElementVector& vec = rememberElVec ? elementVector : userVec;
230

Thomas Witkowski's avatar
Thomas Witkowski committed
231
    if (operat->uhOld && remember) {
232
      matVecAssemble(elInfo, vec);
233
      if (rememberElVec)
234
	userVec += factor * elementVector;      
235

236
237
      return;
    } 
238
239

    if (firstOrderAssemblerGrdPsi)
240
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
241
    if (zeroOrderAssembler)
242
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
243
      
244
    if (rememberElVec)
245
      userVec += factor * elementVector;    
246
247
  }

248

Thomas Witkowski's avatar
Thomas Witkowski committed
249
250
251
252
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
253
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
254
255
256
257
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

258
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
259
260
261
262
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

263
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
264
265
      initElement(smallElInfo, largeElInfo);
   
Thomas Witkowski's avatar
Thomas Witkowski committed
266
    if (el != lastVecEl || !operat->isOptimized()) {
267
268
269
      if (rememberElVec)
	set_to_zero(elementVector);

Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
272
      lastVecEl = el;
    } else {
      if (rememberElVec) {
273
	userVec += factor * elementVector;
Thomas Witkowski's avatar
Thomas Witkowski committed
274
275
276
	return;
      }
    }
277
    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
278
279

    if (operat->uhOld && remember) {
280
      if (smallElInfo->getLevel() == largeElInfo->getLevel())
Thomas Witkowski's avatar
Thomas Witkowski committed
281
	matVecAssemble(auxElInfo, vec);
282
283
      else
	matVecAssemble(mainElInfo, auxElInfo, smallElInfo, largeElInfo, vec);      
Thomas Witkowski's avatar
Thomas Witkowski committed
284

285
      if (rememberElVec)
286
	userVec += factor * elementVector;      
287

Thomas Witkowski's avatar
Thomas Witkowski committed
288
289
290
      return;
    } 

291
292
293
294
295
296
297
298
299
300
    if (firstOrderAssemblerGrdPsi) {
      ERROR_EXIT("Not yet implemented!\n");
    }

    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementVector(smallElInfo, vec);
      
      if (smallElInfo != mainElInfo) {
	ElementVector tmpVec(vec);	
	ElementMatrix &m = 
301
	  smallElInfo->getSubElemCoordsMat(rowFeSpace->getBasisFcts()->getDegree());
302
303
304
305
306

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

308
    if (rememberElVec)
309
      userVec += factor * elementVector;    
Thomas Witkowski's avatar
Thomas Witkowski committed
310
311
  }

312

313
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
314
315
316
  {
    FUNCNAME("Assembler::matVecAssemble()");

317
    Element *el = elInfo->getElement(); 
318
    ElementVector uhOldLoc(operat->uhOld->getFeSpace() == rowFeSpace ? nRow : nCol);
319

320
    operat->uhOld->getLocalVector(el, uhOldLoc);
321
    
322
    if (el != lastMatEl) {
323
      set_to_zero(elementMatrix);
324
      calculateElementMatrix(elInfo, elementMatrix);
325
326
    }

327
    for (int i = 0; i < nRow; i++) {
328
      double val = 0.0;
329
      for (int j = 0; j < nCol; j++)
330
	val += elementMatrix[i][j] * uhOldLoc[j];
331
      
332
      vec[i] += val;
333
    }   
Thomas Witkowski's avatar
Thomas Witkowski committed
334
335
  }

336

Thomas Witkowski's avatar
Thomas Witkowski committed
337
338
  void Assembler::matVecAssemble(const ElInfo *mainElInfo, const ElInfo *auxElInfo,
				 const ElInfo *smallElInfo, const ElInfo *largeElInfo,
339
				 ElementVector& vec)
Thomas Witkowski's avatar
Thomas Witkowski committed
340
341
342
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

346
    TEST_EXIT(operat->uhOld->getFeSpace()->getMesh() == auxElInfo->getMesh())
Thomas Witkowski's avatar
Thomas Witkowski committed
347
348
349
350
351
      ("Da stimmt was nicht!\n");

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

352
    const BasisFunction *basFcts = rowFeSpace->getBasisFcts();
Thomas Witkowski's avatar
Thomas Witkowski committed
353
    int nBasFcts = basFcts->getNumber();
354
    ElementVector uhOldLoc(nBasFcts);
Thomas Witkowski's avatar
Thomas Witkowski committed
355

356
    operat->uhOld->getLocalVector(auxEl, uhOldLoc);
Thomas Witkowski's avatar
Thomas Witkowski committed
357
358

    if (mainEl != lastMatEl) {
359
      set_to_zero(elementMatrix);
360
      calculateElementMatrix(mainElInfo, auxElInfo, smallElInfo, largeElInfo, 
361
			     false, elementMatrix);    
Thomas Witkowski's avatar
Thomas Witkowski committed
362
    }
363

Thomas Witkowski's avatar
Thomas Witkowski committed
364
365
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
366
      for (int j = 0; j < nBasFcts; j++)
367
 	val += elementMatrix[i][j] * uhOldLoc[j];
368
      vec[i] += val;
369
    }   
370
371
  }

372

Thomas Witkowski's avatar
Thomas Witkowski committed
373
374
375
  void Assembler::initElement(const ElInfo *smallElInfo, 
			      const ElInfo *largeElInfo,
			      Quadrature *quad)
376
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
377
    if (secondOrderAssembler) 
Thomas Witkowski's avatar
Thomas Witkowski committed
378
      secondOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
379
    if (firstOrderAssemblerGrdPsi)
380
      firstOrderAssemblerGrdPsi->initElement(smallElInfo, largeElInfo, quad);    
Thomas Witkowski's avatar
Thomas Witkowski committed
381
    if (firstOrderAssemblerGrdPhi)
Thomas Witkowski's avatar
Thomas Witkowski committed
382
      firstOrderAssemblerGrdPhi->initElement(smallElInfo, largeElInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
383
    if (zeroOrderAssembler)
Thomas Witkowski's avatar
Thomas Witkowski committed
384
      zeroOrderAssembler->initElement(smallElInfo, largeElInfo, quad);
385
386
  }

387

388
  void Assembler::checkQuadratures()
Thomas Witkowski's avatar
Thomas Witkowski committed
389
390
  { 
    if (secondOrderAssembler) {
391
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
392
      if (!secondOrderAssembler->getQuadrature()) {
393
	int dim = rowFeSpace->getMesh()->getDim();
394
395
396
397
398
	int degree = operat->getQuadratureDegree(2);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	secondOrderAssembler->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
399
    if (firstOrderAssemblerGrdPsi) {
400
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
401
      if (!firstOrderAssemblerGrdPsi->getQuadrature()) {
402
	int dim = rowFeSpace->getMesh()->getDim();
403
404
405
406
407
	int degree = operat->getQuadratureDegree(1, GRD_PSI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPsi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
408
    if (firstOrderAssemblerGrdPhi) {
409
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
410
      if (!firstOrderAssemblerGrdPhi->getQuadrature()) {
411
	int dim = rowFeSpace->getMesh()->getDim();
412
413
414
415
416
	int degree = operat->getQuadratureDegree(1, GRD_PHI);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	firstOrderAssemblerGrdPhi->setQuadrature(quadrature);
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
417
    if (zeroOrderAssembler) {
418
      // create quadrature
Thomas Witkowski's avatar
Thomas Witkowski committed
419
      if (!zeroOrderAssembler->getQuadrature()) {
420
	int dim = rowFeSpace->getMesh()->getDim();
421
422
423
424
425
426
427
	int degree = operat->getQuadratureDegree(0);
	Quadrature *quadrature = Quadrature::provideQuadrature(dim, degree);
	zeroOrderAssembler->setQuadrature(quadrature);
      }
    }
  }

428

Thomas Witkowski's avatar
Thomas Witkowski committed
429
430
431
432
433
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
434

435

Thomas Witkowski's avatar
Thomas Witkowski committed
436
437
438
439
440
  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
441
442
443
					 const FiniteElemSpace *rowFeSpace,
					 const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
444
  {
445
    bool opt = (rowFeSpace->getBasisFcts() == colFeSpace->getBasisFcts());
Thomas Witkowski's avatar
Thomas Witkowski committed
446
447
448
449
450
451
452
453
454
455
456
457
458
459

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

460

Thomas Witkowski's avatar
Thomas Witkowski committed
461
462
463
464
465
  StandardAssembler::StandardAssembler(Operator *op,
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
466
467
468
				       const FiniteElemSpace *rowFeSpace,
				       const FiniteElemSpace *colFeSpace) 
    : Assembler(op, rowFeSpace, colFeSpace)
Thomas Witkowski's avatar
Thomas Witkowski committed
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
  {
    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();
  }

485
}