Assembler.cc 13 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
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "DOFVector.h"
21
#include "OpenMP.h"
22
23
24

namespace AMDiS {

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

44

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

48

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

55
    if (remember && (factor != 1.0 || operat->uhOld))
56
      rememberElMat = true;
Thomas Witkowski's avatar
Thomas Witkowski committed
57

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

60
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
61
62
63
      initElement(elInfo);

    if (el != lastMatEl || !operat->isOptimized()) {
64
65
66
      if (rememberElMat)
	set_to_zero(elementMatrix);

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

    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
86
    if (rememberElMat && &userMat != &elementMatrix)
87
      userMat += factor * elementMatrix;
88
89
  }

90

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

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

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

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

128
    if (secondOrderAssembler) {
129
      secondOrderAssembler->calculateElementMatrix(smallElInfo, mat);
130

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

    if (firstOrderAssemblerGrdPsi) {
145
146
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);

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

    if (firstOrderAssemblerGrdPhi) {
165
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
166

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

184
185
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
186

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

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

204

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

211
    if (remember && factor != 1.0)
212
213
214
215
      rememberElVec = true;

    Element *el = elInfo->getElement();

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

    ElementVector& vec = rememberElVec ? elementVector : userVec;
232

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

238
239
      return;
    } 
240
241

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

250

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

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

    Element *el = mainElInfo->getElement();

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

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

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

287
      if (rememberElVec)
288
	userVec += factor * elementVector;      
289

Thomas Witkowski's avatar
Thomas Witkowski committed
290
291
292
      return;
    } 

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

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

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

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

314

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

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

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

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

338

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

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

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

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

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

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

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

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

374

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

389

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

430

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

437

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

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

462

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

487
}