SecondOrderAssembler.cc 11.9 KB
Newer Older
1
#include <vector>
2
#include <boost/numeric/mtl/mtl.hpp>
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
#include "Assembler.h"
#include "SecondOrderAssembler.h"
#include "Operator.h"
#include "QPsiPhi.h"
#include "FiniteElemSpace.h"
#include "Quadrature.h"
#include "DOFVector.h"
#include "OpenMP.h"

namespace AMDiS {

  std::vector<SubAssembler*> SecondOrderAssembler::optimizedSubAssemblers; 
  std::vector<SubAssembler*> SecondOrderAssembler::standardSubAssemblers;

  SecondOrderAssembler::SecondOrderAssembler(Operator *op,
					     Assembler *assembler,
					     Quadrature *quad,
					     bool optimized)
    : SubAssembler(op, assembler, quad, 2, optimized)
  {}

  SecondOrderAssembler* 
  SecondOrderAssembler::getSubAssembler(Operator* op,
					Assembler *assembler,
					Quadrature *quad,
					bool optimized)
  {
    int myRank = omp_get_thread_num();

    // check if a assembler is needed at all
33
    if (op->secondOrder[myRank].size() == 0)
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
      return NULL;

    SecondOrderAssembler *newAssembler;

    std::vector<SubAssembler*> *subAssemblers =
	optimized ?
	&optimizedSubAssemblers :
    &standardSubAssemblers;

    std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];

    sort(opTerms.begin(), opTerms.end());

    // check if a new assembler is needed
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
      std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
      sort(assTerms.begin(), assTerms.end());

      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
      }
    }

    // check if all terms are pw_const
    bool pwConst = true;
    for (int i = 0; i < static_cast<int>( op->secondOrder[myRank].size()); i++) {
      if (!op->secondOrder[myRank][i]->isPWConst()) {
	pwConst = false;
	break;
      }
    }  

    // create new assembler
70
    if (!optimized) {
71
      newAssembler = NEW Stand2(op, assembler, quad);
72
73
74
75
76
77
78
    } else {
	if (pwConst) {
	    newAssembler = NEW Pre2(op, assembler, quad);
	} else {
	    newAssembler = NEW Quad2(op, assembler, quad);
	}
    }
79
80
81
82
83
84
85
86
87
88
89

    subAssemblers->push_back(newAssembler);
    return newAssembler;
  }

  Pre2::Pre2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
    q11 = Q11PsiPhi::provideQ11PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
				      owner->getColFESpace()->getBasisFcts(), 
				      quadrature);
90
91
    tmpLALt.resize(omp_get_overall_max_threads());
    for (int i = 0; i < omp_get_overall_max_threads(); i++) {
92
93
94
95
96
97
98
99
      tmpLALt[i] = NEW DimMat<double>*;
      *(tmpLALt[i]) = NEW DimMat<double>(dim, NO_INIT);
    }
  }

  Pre2::~Pre2()
  {
    for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
100
101
      delete *(tmpLALt[i]);
      delete tmpLALt[i];
102
103
104
    }
  }

105
  void Pre2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
106
107
108
109
110
111
  {
    const int **nEntries;
    const int *k, *l;
    const double *values;

    int myRank = omp_get_thread_num();
Thomas Witkowski's avatar
Thomas Witkowski committed
112
    DimMat<double> **LALt = tmpLALt[0];
113
114
115
    DimMat<double> &tmpMat = *LALt[0];
    tmpMat.set(0.0);

116
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
117
118
119
120
121
122
123
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, 1, LALt);
    }

    tmpMat *= elInfo->getDet();
    nEntries = q11->getNumberEntries();

    if (symmetric) {
124
125
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

126
127
128
129
130
131
132
133
      for (int i = 0; i < nRow; i++) {
	k = q11->getKVec(i, i);
	l = q11->getLVec(i, i);
	values = q11->getValVec(i, i);
	double val = 0.0;
	for (int m = 0; m < nEntries[i][i]; m++) {
	  val += values[m] * tmpMat[k[m]][l[m]];
	}
134
	mat[i][i] += val;
135
136
137
138
139
140
141
142
143

	for (int j = i + 1; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
	  values = q11->getValVec(i, j);
	  val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * tmpMat[k[m]][l[m]];
	  }
144
145
	  mat[i][j] += val;
	  mat[j][i] += val;
146
147
148
149
150
151
152
153
154
155
156
157
	}
      }
    } else {  /*  A not symmetric or psi != phi        */
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++) {
	  k = q11->getKVec(i, j);
	  l = q11->getLVec(i, j);
	  values = q11->getValVec(i, j);
	  double val = 0.0;
	  for (int m = 0; m < nEntries[i][j]; m++) {
	    val += values[m] * tmpMat[k[m]][l[m]];
	  }
158
	  mat[i][j] += val;
159
160
161
162
163
164
165
166
	}
      }
    }
  }

  Quad2::Quad2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, true)
  {
167
    tmpLALt.resize(omp_get_overall_max_threads());
168
169
170
171
172
173
174
175
  }

  Quad2::~Quad2()
  {
    if (!firstCall) {
      int nPoints = quadrature->getNumPoints();
      for (int i = 0; i < static_cast<int>(tmpLALt.size()); i++) {
	for (int j = 0; j < nPoints; j++) {
176
	  delete tmpLALt[i][j];
177
	}
178
	delete [] tmpLALt[i];
179
180
181
182
      }
    }
  }

183
  void Quad2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
  {
    int nPoints = quadrature->getNumPoints();
    int myRank = omp_get_thread_num();

    if (firstCall) {
      tmpLALt[myRank] = NEW DimMat<double>*[nPoints];
      for (int j = 0; j < nPoints; j++) {
	tmpLALt[myRank][j] = NEW DimMat<double>(dim, NO_INIT);
      }
    
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

    DimMat<double> **LALt = tmpLALt[myRank];
    
    for (int i = 0; i < nPoints; i++) {
      LALt[i]->set(0.0);
    }
    
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
    }
        
    VectorOfFixVecs< DimVec<double> > *grdPsi, *grdPhi;

    if (symmetric) {
214
215
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

216
217
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();
Thomas Witkowski's avatar
Thomas Witkowski committed
218
	
219
220
221
222
	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

	for (int i = 0; i < nRow; i++) {
223
	  mat[i][i] += quadrature->getWeight(iq) * 
224
225
226
227
228
	    xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[i]);

	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * 
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
229
230
	    mat[i][j] += val;
	    mat[j][i] += val;
231
232
233
234
235
236
237
238
239
240
241
242
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

	grdPsi = psiFast->getGradient(iq);
	grdPhi = phiFast->getGradient(iq);

	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
243
	    mat[i][j] += quadrature->getWeight(iq) * 
244
245
246
247
248
249
250
251
252
253
254
	      xAy((*grdPsi)[i], (*LALt[iq]), (*grdPhi)[j]);
	  }
	}
      }
    }
  }

  Stand2::Stand2(Operator *op, Assembler *assembler, Quadrature *quad) 
    : SecondOrderAssembler(op, assembler, quad, false)
  {}

255
  void Stand2::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix& mat)
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
  {
    DimVec<double> grdPsi(dim, NO_INIT);
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();

    int nPoints = quadrature->getNumPoints();

    DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
      LALt[iq]->set(0.0);
    }

    int myRank = omp_get_thread_num();

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))->getLALt(elInfo, nPoints, LALt);
    }

    if (symmetric) {
278
279
      TEST_EXIT_DBG(nCol == nRow)("nCol != nRow, but symmetric assembling!\n");

280
281
282
283
284
285
286
287
288
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

	for (int i = 0; i < nCol; i++) {
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
	}

	for (int i = 0; i < nRow; i++) {
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
289
	  mat[i][i] += quadrature->getWeight(iq) * 
290
291
292
	    (grdPsi * ((*LALt[iq]) * grdPhi[i]));

	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
293
294
	    double val = quadrature->getWeight(iq) * 
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
295
296
	    mat[i][j] += val;
	    mat[j][i] += val;
297
298
299
300
301
302
303
304
305
306
307
308
309
310
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
      for (int iq = 0; iq < nPoints; iq++) {
	(*LALt[iq]) *= elInfo->getDet();

	for (int i = 0; i < nCol; i++) {
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
	}

	for (int i = 0; i < nRow; i++) {
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	  for (int j = 0; j < nCol; j++) {
311
	    mat[i][j] += quadrature->getWeight(iq) *
312
313
314
315
316
317
318
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }
    }

    for (int iq = 0; iq < nPoints; iq++) 
319
      delete LALt[iq];
320

321
    delete [] LALt;
322
323
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
326
327
  void Stand2::calculateElementMatrix(const ElInfo *rowElInfo, 
				      const ElInfo *colElInfo,
				      const ElInfo *smallElInfo,
				      const ElInfo *largeElInfo,
328
				      ElementMatrix& mat)
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
  {
    FUNCNAME("Stand2::calculateElementMatrix()");

    TEST_EXIT((nRow <= 3) && (nCol <= 3))("not yet!\n");

    bool p = false;
    /*
    if (smallElInfo->getLevel() == largeElInfo->getLevel() + 1) {
      p = true;
      std::cout << "--------------------------------------" << std::endl;
      std::cout << "i-th child = " << smallElInfo->getIChild() << std::endl;
      std::cout << "Mesh row = " << owner->getRowFESpace()->getMesh() << std::endl;
      std::cout << "Mesh col = " << owner->getColFESpace()->getMesh() << std::endl;
      std::cout << "Mesh main = " << rowElInfo->getMesh() << std::endl;
      std::cout << "Mesh aux = " << colElInfo->getMesh() << std::endl;
    }
    */

    DimVec<double> grdPsi(dim, NO_INIT);
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    const BasisFunction *phi = owner->getColFESpace()->getBasisFcts();
352
    const mtl::dense2D<double>& m = smallElInfo->getSubElemCoordsMat();
Thomas Witkowski's avatar
Thomas Witkowski committed
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397

    int nPoints = quadrature->getNumPoints();

    DimMat<double> **LALt = NEW DimMat<double>*[nPoints];
    for (int iq = 0; iq < nPoints; iq++) {
      LALt[iq] = NEW DimMat<double>(dim, NO_INIT);
      LALt[iq]->set(0.0);
    }

    int myRank = omp_get_thread_num();

    if (p) {
      std::cout << "Terms = " << terms[myRank].size() << std::endl;
    }

    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<SecondOrderTerm*>(terms[myRank][i]))
	->getLALt(smallElInfo, nPoints, LALt);
    }

    if (p) {
      for (int i = 0; i < nPoints; i++) {
	std::cout << "LALt for iq = " << i << std::endl;
	LALt[i]->print();
      }
    }

    DimVec<double> val1(dim, DEFAULT_VALUE, 0.0);

    if (p)
      std::cout << "Multiply LALt with det = " << smallElInfo->getDet() << std::endl;

    for (int iq = 0; iq < nPoints; iq++) {
      (*LALt[iq]) *= smallElInfo->getDet();
    }
     
    /*
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {	
	for (int iq = 0; iq < nPoints; iq++) {

	  val1.set(0.0);

	  for (int k = 0; k < nCol; k++) {
	    (*(phi->getGrdPhi(k)))(quadrature->getLambda(iq), grdPsi);
398
	    grdPsi *= m[j][k];
Thomas Witkowski's avatar
Thomas Witkowski committed
399
400
401
402
403
	    val1 += grdPsi;
	  }

	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);

404
	  mat[i][j] += quadrature->getWeight(iq) *
Thomas Witkowski's avatar
Thomas Witkowski committed
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
	      (grdPsi * ((*LALt[iq]) * val1));
	}
      }
      }*/

      for (int iq = 0; iq < nPoints; iq++) {
	for (int i = 0; i < nCol; i++) {
	  (*(phi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPhi[i]);
	  if (p) {
	    std::cout << "grdPhi[" << i << "]" << std::endl;
	    grdPhi[i].print();
	  }
	}

	for (int i = 0; i < nRow; i++) {
	  (*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
	  for (int j = 0; j < nCol; j++) {
422
	    mat[i][j] += quadrature->getWeight(iq) * 
Thomas Witkowski's avatar
Thomas Witkowski committed
423
424
425
426
427
428
	      (grdPsi * ((*LALt[iq]) * grdPhi[j]));
	  }
	}
      }

    for (int iq = 0; iq < nPoints; iq++) 
429
      delete LALt[iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
430

431
    delete [] LALt;
Thomas Witkowski's avatar
Thomas Witkowski committed
432
433
  }

434
}