Assembler.cc 10.5 KB
Newer Older
1
2
3
4
5
6
7
#include "Assembler.h"
#include "Operator.h"
#include "Element.h"
#include "QPsiPhi.h"
#include "ElementMatrix.h"
#include "ElementVector.h"
#include "DOFVector.h"
8
#include "OpenMP.h"
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#include <vector>
#include <algorithm>

namespace AMDiS {

  Assembler::Assembler(Operator  *op,
		       const FiniteElemSpace *rowFESpace_,
		       const FiniteElemSpace *colFESpace_) 
    : operat(op),
      rowFESpace(rowFESpace_),
      colFESpace(colFESpace_ ? colFESpace_ : rowFESpace_),
      nRow(rowFESpace->getBasisFcts()->getNumber()),
      nCol(colFESpace->getBasisFcts()->getNumber()),
      remember(true),
      rememberElMat(false),
      rememberElVec(false),
      elementMatrix(NULL),
      elementVector(NULL),
      lastMatEl(NULL),
      lastVecEl(NULL),
      lastTraverseId(-1)
  
Thomas Witkowski's avatar
Thomas Witkowski committed
31
  {}
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47

  void Assembler::calculateElementMatrix(const ElInfo *elInfo, 
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = elInfo->getElement();

Thomas Witkowski's avatar
Thomas Witkowski committed
48
    
49
50
51
    checkForNewTraverse();

    checkQuadratures();
Thomas Witkowski's avatar
Thomas Witkowski committed
52
    
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(elInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

    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
80
    if (rememberElMat && userMat) {
81
      axpy(factor, *elementMatrix, *userMat);
82
    }      
83
84
  }

85
86
  void Assembler::calculateElementMatrix(const ElInfo *rowElInfo,
					 const ElInfo *colElInfo,
87
88
					 const ElInfo *smallElInfo,
					 const ElInfo *largeElInfo,
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
					 ElementMatrix *userMat,
					 double factor)
  {
    FUNCNAME("Assembler::calculateElementMatrix()");

    if (remember && ((factor != 1.0) || (operat->uhOld))) {
      rememberElMat = true;
    }
  
    if (rememberElMat && !elementMatrix)
      elementMatrix = NEW ElementMatrix(nRow, nCol);

    Element *el = rowElInfo->getElement();

    
    //    checkForNewTraverse();
    lastVecEl = lastMatEl = NULL;

    checkQuadratures();
    
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
      initElement(rowElInfo);
    }

    if (el != lastMatEl || !operat->isOptimized()) {
      if (rememberElMat) {
	elementMatrix->set(0.0);
      }
      lastMatEl = el;
    } else {
      if (rememberElMat) {
	axpy(factor, *elementMatrix, *userMat);
	return;
      }
    }
  
    ElementMatrix *mat = rememberElMat ? elementMatrix : userMat;

127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
    if (secondOrderAssembler) {
      secondOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
						   smallElInfo, largeElInfo, mat);
    }

    if (firstOrderAssemblerGrdPsi) {
      firstOrderAssemblerGrdPsi->calculateElementMatrix(rowElInfo, colElInfo, 
							smallElInfo, largeElInfo, mat);
    }

    if (firstOrderAssemblerGrdPhi) {
      firstOrderAssemblerGrdPhi->calculateElementMatrix(rowElInfo, colElInfo, 
							smallElInfo, largeElInfo, mat);
    }

142
    if (zeroOrderAssembler) {
143
144
      zeroOrderAssembler->calculateElementMatrix(rowElInfo, colElInfo, 
						 smallElInfo, largeElInfo, mat);
145
    }
146
 
147
148
149
150
151
    if (rememberElMat && userMat) {
      axpy(factor, *elementMatrix, *userMat);
    }      
  }

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

Thomas Witkowski's avatar
Thomas Witkowski committed
158
    if (remember && factor != 1.0) {
159
160
161
      rememberElVec = true;
    }

162
    if (rememberElVec && !elementVector) {
163
      elementVector = NEW ElementVector(nRow);
164
    }
165
166
167
168
169

    Element *el = elInfo->getElement();
    checkForNewTraverse();
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
170
    if ((el != lastMatEl && el != lastVecEl) || !operat->isOptimized()) {
171
172
      initElement(elInfo);
    }
173
    
Thomas Witkowski's avatar
Thomas Witkowski committed
174
175
    if (el != lastVecEl || !operat->isOptimized()) {
      if (rememberElVec) {
176
177
178
179
	elementVector->set(0.0);
      }
      lastVecEl = el;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
180
      if (rememberElVec) {
181
182
183
184
185
	axpy(factor, *elementVector, *userVec);
	return;
      }
    }
    ElementVector *vec = rememberElVec ? elementVector : userVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
186
    if (operat->uhOld && remember) {
187
      matVecAssemble(elInfo, vec);
Thomas Witkowski's avatar
Thomas Witkowski committed
188
      if (rememberElVec) {
189
190
191
192
	axpy(factor, *elementVector, *userVec);
      }
      return;
    } 
193
    if (firstOrderAssemblerGrdPsi) {
194
      firstOrderAssemblerGrdPsi->calculateElementVector(elInfo, vec);
195
196
    }
    if (zeroOrderAssembler) {
197
      zeroOrderAssembler->calculateElementVector(elInfo, vec);
198
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
199
    if (rememberElVec) {
200
      axpy(factor, *elementVector, *userVec);
201
    }      
202
203
204
205
206
207
208
209
  }

  void Assembler::matVecAssemble(const ElInfo *elInfo, ElementVector *vec)
  {
    FUNCNAME("Assembler::matVecAssemble()");

    Element *el = elInfo->getElement(); 
    const BasisFunction *basFcts = rowFESpace->getBasisFcts();
210
211
    int nBasFcts = basFcts->getNumber();
    double *uhOldLoc = new double[nBasFcts];
212

213
214
    operat->uhOld->getLocalVector(el, uhOldLoc);
    
Thomas Witkowski's avatar
Thomas Witkowski committed
215
    if (el != lastMatEl) {
216
217
      calculateElementMatrix(elInfo, NULL);
    }
218
219
220
221
222
    
    for (int i = 0; i < nBasFcts; i++) {
      double val = 0.0;
      for (int j = 0; j < nBasFcts; j++) {
	val += (*elementMatrix)[i][j] * uhOldLoc[j];
223
224
225
      }
      (*vec)[i] += val;
    }
226
227
228
    

    delete [] uhOldLoc;       
229
230
231
232
233
234
  }

  void Assembler::initElement(const ElInfo *elInfo, Quadrature *quad)
  {
    checkQuadratures();

Thomas Witkowski's avatar
Thomas Witkowski committed
235
    if (secondOrderAssembler) 
236
      secondOrderAssembler->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
237
    if (firstOrderAssemblerGrdPsi)
238
      firstOrderAssemblerGrdPsi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
239
    if (firstOrderAssemblerGrdPhi)
240
      firstOrderAssemblerGrdPhi->initElement(elInfo, quad);
Thomas Witkowski's avatar
Thomas Witkowski committed
241
    if (zeroOrderAssembler)
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
      zeroOrderAssembler->initElement(elInfo, quad);
  }

  OptimizedAssembler::OptimizedAssembler(Operator  *op,
					 Quadrature *quad2,
					 Quadrature *quad1GrdPsi,
					 Quadrature *quad1GrdPhi,
					 Quadrature *quad0,
					 const FiniteElemSpace *rowFESpace_,
					 const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    bool opt = (rowFESpace_ == colFESpace_);

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
257
    secondOrderAssembler = 
258
      SecondOrderAssembler::getSubAssembler(op, this, quad2, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
259
    firstOrderAssemblerGrdPsi = 
260
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
261
    firstOrderAssemblerGrdPhi = 
262
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, opt);
Thomas Witkowski's avatar
Thomas Witkowski committed
263
    zeroOrderAssembler = 
264
265
266
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, opt);
  }

267
  StandardAssembler::StandardAssembler(Operator *op,
268
269
270
271
272
273
274
275
276
277
278
				       Quadrature *quad2,
				       Quadrature *quad1GrdPsi,
				       Quadrature *quad1GrdPhi,
				       Quadrature *quad0,
				       const FiniteElemSpace *rowFESpace_,
				       const FiniteElemSpace *colFESpace_) 
    : Assembler(op, rowFESpace_, colFESpace_)
  {
    remember = false;

    // create sub assemblers
Thomas Witkowski's avatar
Thomas Witkowski committed
279
    secondOrderAssembler = 
280
      SecondOrderAssembler::getSubAssembler(op, this, quad2, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    firstOrderAssemblerGrdPsi = 
282
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPsi, GRD_PSI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
283
    firstOrderAssemblerGrdPhi = 
284
      FirstOrderAssembler::getSubAssembler(op, this, quad1GrdPhi, GRD_PHI, false);
Thomas Witkowski's avatar
Thomas Witkowski committed
285
    zeroOrderAssembler = 
286
287
288
289
      ZeroOrderAssembler::getSubAssembler(op, this, quad0, false);
  }

  ElementMatrix *Assembler::initElementMatrix(ElementMatrix *elMat, 
290
291
					      const ElInfo *rowElInfo,
					      const ElInfo *colElInfo)
292
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
    if (!elMat) {
294
295
296
297
      elMat = NEW ElementMatrix(nRow, nCol);
    }

    elMat->set(0.0);
Thomas Witkowski's avatar
Thomas Witkowski committed
298
         
299
    rowFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
Thomas Witkowski's avatar
Thomas Witkowski committed
300
						   rowFESpace->getAdmin(),
Thomas Witkowski's avatar
Thomas Witkowski committed
301
						   &(elMat->rowIndices));
302

Thomas Witkowski's avatar
Thomas Witkowski committed
303
304
305
    if (rowFESpace == colFESpace) {
      elMat->colIndices = elMat->rowIndices;
    } else {
306
307
308
309
310
311
312
313
314
315
      if (colElInfo) {
	colFESpace->getBasisFcts()->getLocalIndicesVec(colElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      } else {
	// If there is no colElInfo pointer, use rowElInfo the get the indices.
	colFESpace->getBasisFcts()->getLocalIndicesVec(rowElInfo->getElement(),
						       colFESpace->getAdmin(),
						       &(elMat->colIndices));
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
316
    }
317
    
318
319
320
321
322
323
    return elMat;
  }

  ElementVector *Assembler::initElementVector(ElementVector *elVec, 
					      const ElInfo *elInfo)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
324
    if (!elVec) {
325
326
327
328
329
330
331
332
      elVec = NEW ElementVector(nRow);
    }

    elVec->set(0.0);
  
    DOFAdmin *rowAdmin = rowFESpace->getAdmin();
    Element *element = elInfo->getElement();

333
    rowFESpace->getBasisFcts()->getLocalIndicesVec(element, rowAdmin, &(elVec->dofIndices));
334
335
336
337
338

    return elVec;
  }

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

}