Assembler.cc 12.1 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()");

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

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

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

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

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

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

76

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

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

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

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

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

      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat_so(m, rowFESpace->getBasisFcts()->getDegree());
114

115
116
      ElementMatrix tmpMat(nRow, nRow);
      tmpMat = m * mat;
117
      mat = tmpMat;      
118
119
120
    }

    if (firstOrderAssemblerGrdPsi) {
121
      ERROR_EXIT("Not yet implemented for first order assembler!\n");
122
      firstOrderAssemblerGrdPsi->calculateElementMatrix(smallElInfo, mat);
123
124
125
    }

    if (firstOrderAssemblerGrdPhi) {
126
      ERROR_EXIT("Not yet implemented for first order assembler!\n");
127
      firstOrderAssemblerGrdPhi->calculateElementMatrix(smallElInfo, mat);
128
    }
129

130
131
    if (zeroOrderAssembler) {
      zeroOrderAssembler->calculateElementMatrix(smallElInfo, mat);
132

133
134
      ElementMatrix m(nRow, nRow);
      smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
135
136
137
138
139
140
141
142

#if 0
      std::cout << "ASM MAT:" << std::endl;
      std::cout << mat << std::endl;
      std::cout << "MULT MAT:" << std::endl;
      std::cout << m << std::endl;
#endif

143
      ElementMatrix tmpMat(nRow, nRow);
144
145
      //tmpMat = m * mat;
      tmpMat = mat * m;
146
147
      mat = tmpMat;
    }
148

149
150
    if (rememberElMat)
      userMat += factor * elementMatrix;
151
152
  }

153

154
  void Assembler::calculateElementVector(const ElInfo *elInfo, 
155
					 ElementVector& userVec,
156
157
158
159
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

160
    if (remember && factor != 1.0)
161
162
163
164
      rememberElVec = true;

    Element *el = elInfo->getElement();

165
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
166
      initElement(elInfo);
167
    
Thomas Witkowski's avatar
Thomas Witkowski committed
168
    if (el != lastVecEl || !operat->isOptimized()) {
169
170
171
      if (rememberElVec)
	set_to_zero(elementVector);
	
172
173
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
174
      if (rememberElVec) {
175
	userVec += factor * elementVector;
176
177
178
	return;
      }
    }
179
180

    ElementVector& vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
181
    if (operat->uhOld && remember) {
182
      matVecAssemble(elInfo, vec);
183
      if (rememberElVec)
184
	userVec += factor * elementVector;
185

186
187
      return;
    } 
188
189

    if (firstOrderAssemblerGrdPsi)
190
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
191
    if (zeroOrderAssembler)
192
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
193
      
194
195
    if (rememberElVec)
      userVec += factor * elementVector;
196
197
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
198
199
200
201
  void Assembler::calculateElementVector(const ElInfo *mainElInfo, 
					 const ElInfo *auxElInfo,
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
202
					 ElementVector& userVec, 
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
205
206
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementVector()");

207
    if (remember && factor != 1.0)
Thomas Witkowski's avatar
Thomas Witkowski committed
208
209
210
211
      rememberElVec = true;

    Element *el = mainElInfo->getElement();

212
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized())
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
215
      initElement(auxElInfo);
    
    if (el != lastVecEl || !operat->isOptimized()) {
216
217
218
      if (rememberElVec)
	set_to_zero(elementVector);

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

    if (operat->uhOld && remember) {

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

236
      if (rememberElVec)
237
238
	userVec += factor * elementVector;

Thomas Witkowski's avatar
Thomas Witkowski committed
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
      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);
//     }      

  }

255
  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector& vec)
256
257
258
259
  {
    FUNCNAME("Assembler::matVecAssemble()");

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
275
276
277
278
279
    delete [] uhOldLoc;
  }

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

300
301
    mtl::dense2D<double> m(nRow, nRow);
    smallElInfo->getSubElemCoordsMat(m, rowFESpace->getBasisFcts()->getDegree());
Thomas Witkowski's avatar
Thomas Witkowski committed
302
303
304
305

    for (int i = 0; i < nBasFcts; i++) {
      uhOldLoc2[i] = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
306
	uhOldLoc2[i] += m[j][i] * uhOldLoc[i];
Thomas Witkowski's avatar
Thomas Witkowski committed
307
308
309
310
311
      }      
    }

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
382
383
384
385
386
  void Assembler::finishAssembling()
  {
    lastVecEl = NULL;
    lastMatEl = NULL;
  }
Thomas Witkowski's avatar
Thomas Witkowski committed
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
432
433
434
435

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

436
}