Assembler.cc 45.4 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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
#include <vector>
#include <algorithm>

namespace AMDiS {

  ::std::vector<SubAssembler*> ZeroOrderAssembler::optimizedSubAssemblers;
  ::std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPhi;
  ::std::vector<SubAssembler*> FirstOrderAssembler::optimizedSubAssemblersGrdPsi;
  ::std::vector<SubAssembler*> SecondOrderAssembler::optimizedSubAssemblers;
  
  ::std::vector<SubAssembler*> ZeroOrderAssembler::standardSubAssemblers;
  ::std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPhi;
  ::std::vector<SubAssembler*> FirstOrderAssembler::standardSubAssemblersGrdPsi;
  ::std::vector<SubAssembler*> SecondOrderAssembler::standardSubAssemblers;

  SubAssembler::SubAssembler(Operator *op,
			     Assembler *assembler,
			     Quadrature *quadrat,
			     int order, 
			     bool optimized,
			     FirstOrderType type) 
    : nRow(0),
      nCol(0),
      coordsAtQPs(NULL),
      quadrature(quadrat),
      psiFast(NULL),
      phiFast(NULL),
      owner(assembler),
      symmetric(true),
      opt(optimized),
      firstCall(true)
  {
    const BasisFunction *psi = assembler->rowFESpace->getBasisFcts();
    const BasisFunction *phi = assembler->colFESpace->getBasisFcts();

    nRow = psi->getNumber();
    nCol = phi->getNumber();

47
48
49
    int maxThreads = omp_get_max_threads();
    terms.resize(maxThreads);

50
    switch (order) {
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
    case 0:
      terms = op->zeroOrder;
      break;
    case 1:
      if(type == GRD_PHI)
	terms = op->firstOrderGrdPhi;
      else 
	terms = op->firstOrderGrdPsi;
      break;
    case 2:
      terms = op->secondOrder;
      break;
    }

    // check if all terms are symmetric
    symmetric = true;
67
68
    for (int i = 0; i < static_cast<int>(terms[0].size()); i++) {
      if (!(terms[0][i])->isSymmetric()) {
69
70
71
72
73
74
75
76
77
78
79
80
	symmetric = false;
	break;
      }
    }  

    dim = assembler->rowFESpace->getMesh()->getDim();
  }

  FastQuadrature *SubAssembler::updateFastQuadrature(FastQuadrature *quadFast,
						     const BasisFunction *psi,
						     Flag updateFlag)
  {
81
82
83
84
    if (!quadFast) {
      quadFast = FastQuadrature::provideFastQuadrature(psi,
						       *quadrature,
						       updateFlag);
85
    } else {
86
      if (!quadFast->initialized(updateFlag))
87
88
89
90
91
92
93
94
95
96
97
98
99
100
	quadFast->init(updateFlag);
    }

    return quadFast;
  }

  void SubAssembler::initElement(const ElInfo* elInfo,
				 Quadrature *quad)
  {
    // set corrdsAtQPs invalid
    coordsValid = false;

    // set values at QPs invalid
    ::std::map<const DOFVectorBase<double>*, ValuesAtQPs*>::iterator it1;
101
    for (it1 = valuesAtQPs.begin(); it1 != valuesAtQPs.end(); ++it1) {
102
103
104
105
106
      ((*it1).second)->valid = false;
    }

    // set gradients at QPs invalid
    ::std::map<const DOFVectorBase<double>*, GradientsAtQPs*>::iterator it2;
107
    for (it2 = gradientsAtQPs.begin(); it2 != gradientsAtQPs.end(); ++it2) {
108
109
110
      ((*it2).second)->valid = false;
    }

111
    int myRank = omp_get_thread_num();
112
113
    // calls initElement of each term
    ::std::vector<OperatorTerm*>::iterator it;
114
    for (it = terms[myRank].begin(); it != terms[myRank].end(); ++it) {
115
116
117
118
119
120
121
122
123
124
125
126
      (*it)->initElement(elInfo, this, quad);
    }
  }

  WorldVector<double>* SubAssembler::getCoordsAtQPs(const ElInfo* elInfo,
						    Quadrature *quad)
  {
    Quadrature *localQuad = quad ? quad : quadrature;
  
    const int numPoints = localQuad->getNumPoints();

    // already calculated for this element ?
127
    if (coordsValid) {
128
129
130
131
      return coordsAtQPs;
    }
   
    // not yet calcualted !
132
133
    if (coordsAtQPs) 
      DELETE [] coordsAtQPs;
134
135
136
137
138
139
    // allocate memory
    coordsAtQPs = NEW WorldVector<double>[numPoints];

    // set new values
    WorldVector<double>* k = NULL;
    int l;
140
    for (k = &(coordsAtQPs[0]), l = 0; k < &(coordsAtQPs[numPoints]); ++k, ++l) {
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
      elInfo->coordToWorld(localQuad->getLambda(l), k);
    }

    // mark values as valid
    coordsValid = true;

    return coordsAtQPs;
  }

  double* SubAssembler::getVectorAtQPs(DOFVectorBase<double>* dv, 
				       const ElInfo* elInfo,
				       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getVectorAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

158
    TEST_EXIT_DBG(vec)("no dof vector!\n");
159

160
    if (valuesAtQPs[vec] && valuesAtQPs[vec]->valid) 
161
162
163
164
      return valuesAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

165
    if (!valuesAtQPs[vec]) {
166
167
168
169
170
171
172
173
174
175
      valuesAtQPs[vec] = new ValuesAtQPs;
    }
    valuesAtQPs[vec]->values.resize(localQuad->getNumPoints());

    double *values = valuesAtQPs[vec]->values.getValArray();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

176
    if (opt && !quad && sameFESpaces) {
177
178
      const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
      const BasisFunction *phi = owner->colFESpace->getBasisFcts();
179
      if (vec->getFESpace()->getBasisFcts() == psi) {
180
181
182
183
184
185
186
187
188
	psiFast = updateFastQuadrature(psiFast, psi, INIT_PHI);
      } else if(vec->getFESpace()->getBasisFcts() == phi) {
	phiFast = updateFastQuadrature(phiFast, phi, INIT_PHI);
      }
    }

    // calculate new values
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();

189
190
    if (opt && !quad && sameFESpaces) {
      if (psiFast->getBasisFunctions() == basFcts) {
191
192
193
194
195
196
197
198
199
	vec->getVecAtQPs(elInfo, NULL, psiFast, values);
      } else if(phiFast->getBasisFunctions() == basFcts) {
	vec->getVecAtQPs(elInfo, NULL, phiFast, values);
      } else {
	vec->getVecAtQPs(elInfo, localQuad, NULL, values);
      }
    } else {
      vec->getVecAtQPs(elInfo, localQuad, NULL, values);
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
200
  
201
202
203
204
205
206
207
208
209
210
211
212
213
    valuesAtQPs[vec]->valid = true;

    return values;
  }

  WorldVector<double>* SubAssembler::getGradientsAtQPs(DOFVectorBase<double>* dv, 
						       const ElInfo* elInfo,
						       Quadrature *quad)
  {
    FUNCNAME("SubAssembler::getGradientsAtQPs()");

    const DOFVectorBase<double>* vec = dv ? dv : owner->operat->getUhOld();

214
    TEST_EXIT_DBG(vec)("no dof vector!\n");
215

216
    if (gradientsAtQPs[vec] && gradientsAtQPs[vec]->valid) 
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
      return gradientsAtQPs[vec]->values.getValArray();

    Quadrature *localQuad = quad ? quad : quadrature;

    if(!gradientsAtQPs[vec]) {
      gradientsAtQPs[vec] = new GradientsAtQPs;
    } 
    gradientsAtQPs[vec]->values.resize(localQuad->getNumPoints());

    WorldVector<double> *values = gradientsAtQPs[vec]->values.getValArray();

    const BasisFunction *psi = owner->rowFESpace->getBasisFcts();
    const BasisFunction *phi = owner->colFESpace->getBasisFcts();

    bool sameFESpaces = 
      (vec->getFESpace() == owner->rowFESpace) || 
      (vec->getFESpace() == owner->colFESpace);

235
236
    if (opt && !quad && sameFESpaces) {
      if (vec->getFESpace()->getBasisFcts() == psi) {
237
238
239
240
241
242
243
244
245
	psiFast = updateFastQuadrature(psiFast, psi, INIT_GRD_PHI);
      } else if(vec->getFESpace()->getBasisFcts() == phi) {
	phiFast = updateFastQuadrature(phiFast, phi, INIT_GRD_PHI);
      }
    }
  
    // calculate new values
    const BasisFunction *basFcts = vec->getFESpace()->getBasisFcts();

246
    if (opt && !quad && sameFESpaces) {
247
248
249
250
251
252
253
254
      if(psiFast->getBasisFunctions() == basFcts) {
	vec->getGrdAtQPs(elInfo, NULL, psiFast, values);
      } else {
	vec->getGrdAtQPs(elInfo, NULL, phiFast, values);
      }
    } else {
      vec->getGrdAtQPs(elInfo, localQuad, NULL, values);
    }
255
   
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
    gradientsAtQPs[vec]->valid = true;

    return values;
  }

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

  FirstOrderAssembler::FirstOrderAssembler(Operator *op,
					   Assembler *assembler,
					   Quadrature *quad,
					   bool optimized,
					   FirstOrderType type)
    : SubAssembler(op, assembler, quad, 1, optimized, type)
  {}

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

  ZeroOrderAssembler* 
  ZeroOrderAssembler::getSubAssembler(Operator* op,
				      Assembler *assembler,
				      Quadrature *quad,
				      bool optimized)
  {
    // check if a assembler is needed at all
290
    if (op->zeroOrder.size() == 0) {
291
292
293
294
295
296
297
298
299
300
      return NULL;
    }

    ZeroOrderAssembler *newAssembler;

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

301
302
    int myRank = omp_get_thread_num();
    ::std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];
303
304
305
306

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

    // check if a new assembler is needed
307
308
    if (quad) {
      for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
309
310
311
312
	::std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());

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

313
314
315
316
317
	if ((opTerms == assTerms) && 
	    ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	  return dynamic_cast<ZeroOrderAssembler*>((*subAssemblers)[i]);
	}
318
319
      }
    }
320
 
321
322
    // check if all terms are pw_const
    bool pwConst = true;
323
324
    for (int i = 0; i < static_cast<int>( op->zeroOrder[myRank].size()); i++) {
      if (!op->zeroOrder[myRank][i]->isPWConst()) {
325
326
327
328
329
330
	pwConst = false;
	break;
      }
    }  

    // create new assembler
331
    if (!optimized) {
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
      newAssembler = NEW Stand0(op, assembler, quad);
    } else {
      if(pwConst) {
	newAssembler = NEW Pre0(op, assembler, quad);
      } else {
	newAssembler = NEW Quad0(op, assembler, quad);
      }
    }

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

  FirstOrderAssembler* 
  FirstOrderAssembler::getSubAssembler(Operator* op,
				       Assembler *assembler,
				       Quadrature *quad,
				       FirstOrderType type,
				       bool optimized)
  {
    ::std::vector<SubAssembler*> *subAssemblers =
	optimized ?
	(type == GRD_PSI ? 
	 &optimizedSubAssemblersGrdPsi : 
	 &optimizedSubAssemblersGrdPhi) :
    (type == GRD_PSI ? 
     &standardSubAssemblersGrdPsi :
     &standardSubAssemblersGrdPhi);

361
    int myRank = omp_get_thread_num();
362
    ::std::vector<OperatorTerm*> opTerms 
363
364
365
	= (type == GRD_PSI) ? 
	    op->firstOrderGrdPsi[myRank] : 
            op->firstOrderGrdPhi[myRank];
366
367

    // check if a assembler is needed at all
368
    if (opTerms.size() == 0) {
369
370
371
372
373
374
375
376
      return NULL;
    }

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

    FirstOrderAssembler *newAssembler;

    // check if a new assembler is needed
377
    for (int i = 0; i < static_cast<int>( subAssemblers->size()); i++) {
378
379
380

      ::std::vector<OperatorTerm*> assTerms = *((*subAssemblers)[i]->getTerms());
    
381
      sort(assTerms.begin(), assTerms.end());   
382

383
384
      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
385

386
387
	return dynamic_cast<FirstOrderAssembler*>((*subAssemblers)[i]);
      }
388
389
390
391
    }

    // check if all terms are pw_const
    bool pwConst = true;
392
393
    for (int i = 0; i < static_cast<int>( opTerms.size()); i++) {
      if (!(opTerms[i])->isPWConst()) {
394
395
396
397
398
399
	pwConst = false;
	break;
      }
    }  

    // create new assembler
400
    if (!optimized) {
401
402
403
404
405
      newAssembler = 
	(type == GRD_PSI) ?
	dynamic_cast<FirstOrderAssembler*>(NEW Stand10(op, assembler, quad)) :
	dynamic_cast<FirstOrderAssembler*>(NEW Stand01(op, assembler, quad));    
    } else {
406
      if (pwConst) {
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>(NEW Pre01(op, assembler, quad));
      } else {
	newAssembler = 
	  (type == GRD_PSI) ?
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad10(op, assembler, quad)) :
	  dynamic_cast<FirstOrderAssembler*>( NEW Quad01(op, assembler, quad));
      }
    }

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

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

431
    // check if a assembler is needed at all
432
    if (op->secondOrder[myRank].size() == 0) {
433
434
435
436
437
438
439
440
441
442
      return NULL;
    }

    SecondOrderAssembler *newAssembler;

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

443
    ::std::vector<OperatorTerm*> opTerms  = op->zeroOrder[myRank];
444
445
446
447

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

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

453
454
455
456
457
      if ((opTerms == assTerms) && 
	  ((*subAssemblers)[i]->getQuadrature() == quad)) {
	
	return dynamic_cast<SecondOrderAssembler*>((*subAssemblers)[i]);
      }
458
459
460
461
    }

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

    // create new assembler
470
    if (!optimized) {
471
472
      newAssembler = NEW Stand2(op, assembler, quad);
    } else {
473
      if (pwConst) {
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
	newAssembler = NEW Pre2(op, assembler, quad);
      } else {
	newAssembler = NEW Quad2(op, assembler, quad);
      }
    }

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

  Stand0::Stand0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, false)
  {
  }

  void Stand0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    double val;

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

    double  psival;
    double *phival = GET_MEMORY(double, nCol);
    int numPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, numPoints);
Thomas Witkowski's avatar
Thomas Witkowski committed
500

501
    for (int iq = 0; iq < numPoints; iq++) {
502
503
504
      c[iq] = 0.0;
    }

505
    int myRank = omp_get_thread_num();
506
    ::std::vector<OperatorTerm*>::iterator termIt;
507
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
508
509
510
511
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
    }
      
    if (symmetric) {
512
      for (int iq = 0; iq < numPoints; iq++) {
513
514
515
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
516
	for (int i = 0; i < nCol; i++) {
517
518
519
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

520
	for (int i = 0; i < nRow; i++) {
521
	  psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
Thomas Witkowski's avatar
Thomas Witkowski committed
522
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psival * phival[i];
523
	  for (int j = i + 1; j < nCol; j++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
524
	    val = quadrature->getWeight(iq) * c[iq] * psival * phival[j];
525
526
527
528
529
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
Thomas Witkowski's avatar
Thomas Witkowski committed
530
    } else {      //  non symmetric assembling 
531
      for (int iq = 0; iq < numPoints; iq++) {
532
533
534
	c[iq] *= elInfo->getDet();

	// calculate phi at QPs only once!
535
	for (int i = 0; i < nCol; i++) {
536
537
538
	  phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
	}

539
	for (int i = 0; i < nRow; i++) {
540
	  psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
541
	  for (int j = 0; j < nCol; j++) {
542
543
544
545
546
	    (*mat)[i][j] += quadrature->getWeight(iq)*c[iq]*psival*phival[j];
	  }
	}
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
547

548
    FREE_MEMORY(phival, double, nCol);
Thomas Witkowski's avatar
Thomas Witkowski committed
549
    FREE_MEMORY(c, double, numPoints);
550
551
552
553
554
555
556
  }

  void Stand0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    int numPoints = quadrature->getNumPoints();

    double *c = GET_MEMORY(double, numPoints);
557

558
    for (int iq = 0; iq < numPoints; iq++) {
559
560
561
      c[iq] = 0.0;
    }

562
    int myRank = omp_get_thread_num();
563
    ::std::vector<OperatorTerm*>::iterator termIt;
564
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
565
566
567
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
    }

568
    for (int iq = 0; iq < numPoints; iq++) {
569
570
      c[iq] *= elInfo->getDet();

571
572
      for (int i = 0; i < nRow; i++) {
	double psi = (*(owner->getRowFESpace()->getBasisFcts()->getPhi(i)))
573
	  (quadrature->getLambda(iq));
574
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi;
575
576
      }
    }
577
    
578
579
580
581
582
583
584
585
586
587
588
589
    FREE_MEMORY(c, double, numPoints);
  }

  Quad0::Quad0(Operator *op, Assembler *assembler, Quadrature *quad)
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

  void Quad0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    const double *psi, *phi;

590
    if (firstCall) {
591
592
593
594
595
596
597
598
599
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int numPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, numPoints);
600

601
    for (int iq = 0; iq < numPoints; iq++) {
602
603
604
      c[iq] = 0.0;
    }

605
    int myRank = omp_get_thread_num();
606
    ::std::vector<OperatorTerm*>::iterator termIt;
607
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
608
609
610
611
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
    }

    if (symmetric) {
612
      for (int iq = 0; iq < numPoints; iq++) {
613
614
615
616
	c[iq] *= elInfo->getDet();

	psi = psiFast->getPhi(iq);
	phi = phiFast->getPhi(iq);
617
618
619
620
	for (int i = 0; i < nRow; i++) {
	  (*mat)[i][i] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[i];
	  for (int j = i + 1; j < nCol; j++) {
	    double val = quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
621
622
623
624
625
626
	    (*mat)[i][j] += val;
	    (*mat)[j][i] += val;
	  }
	}
      }
    } else {      /*  non symmetric assembling   */
627
      for (int iq = 0; iq < numPoints; iq++) {
628
629
630
631
	c[iq] *= elInfo->getDet();

	psi = psiFast->getPhi(iq);
	phi = phiFast->getPhi(iq);
632
633
634
	for (int i = 0; i < nRow; i++) {
	  for (int j = 0; j < nCol; j++) {
	    (*mat)[i][j] += quadrature->getWeight(iq) * c[iq] * psi[i] * phi[j];
635
636
637
638
639
640
641
642
643
	  }
	}
      }
    }
    FREE_MEMORY(c, double, numPoints);
  }

  void Quad0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
644
    if (firstCall) {
645
646
647
648
649
650
651
652
653
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int numPoints = quadrature->getNumPoints();
    double *c = GET_MEMORY(double, numPoints);
654

655
    for (int iq = 0; iq < numPoints; iq++) {
656
657
658
      c[iq] = 0.0;
    }

659
    int myRank = omp_get_thread_num();
660
    ::std::vector<OperatorTerm*>::iterator termIt;
661
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
662
663
664
      (static_cast<ZeroOrderTerm*>((*termIt)))->getC(elInfo, numPoints, c);
    }

665
    for (int iq = 0; iq < numPoints; iq++) {
666
667
      c[iq] *= elInfo->getDet();

668
669
670
      const double *psi = psiFast->getPhi(iq);
      for (int i = 0; i < nRow; i++) {
	(*vec)[i] += quadrature->getWeight(iq) * c[iq] * psi[i];
671
672
673
674
675
676
677
678
679
680
681
682
      }
    }
    FREE_MEMORY(c, double, numPoints);
  }

  Pre0::Pre0(Operator *op, Assembler *assembler, Quadrature *quad) 
    : ZeroOrderAssembler(op, assembler, quad, true)
  {
  }

  void Pre0::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
683
    if (firstCall) {
684
685
686
687
688
689
690
691
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

692
693
    //    c[0] = 0.0;
    double c = 0.0;
694
    int myRank = omp_get_thread_num();
695
    int size = static_cast<int>(terms[myRank].size());
696

697
698
    for (int i = 0; i < size; i++) {
      (static_cast<ZeroOrderTerm*>((terms[myRank][i])))->getC(elInfo, 1, &c);
699
700
    }

701
    c *= elInfo->getDet();
702
703

    if (symmetric) {
704
      for (int i = 0; i < nRow; i++) {
705
	(*mat)[i][i] += c * q00->getValue(i,i);
706
	for (int j = i + 1; j < nCol; j++) {
707
	  double val = c * q00->getValue(i, j);
708
709
710
711
712
	  (*mat)[i][j] += val;
	  (*mat)[j][i] += val;
	}
      }
    } else {
713
714
      for (int i = 0; i < nRow; i++)
	for (int j = 0; j < nCol; j++)
715
	  (*mat)[i][j] += c * q00->getValue(i, j);
716
717
718
719
720
    }
  }

  void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
721
    if (firstCall) {
722
723
724
725
726
727
728
729
730
731
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    ::std::vector<OperatorTerm*>::iterator termIt;

732
    int myRank = omp_get_thread_num();
733
    double c = 0.0;
734
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
735
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, &c);
736
737
    }

738
    c *= elInfo->getDet();
739

740
    for (int i = 0; i < nRow; i++)
741
      (*vec)[i] += c * q0->getValue(i);
742
743
744
745
746
747
748
749
750
  }

  Stand10::Stand10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PSI)
  {}


  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
751
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
752
753
754
755
756
757
758
    double *phival = GET_MEMORY(double, nCol);

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

    int numPoints = quadrature->getNumPoints();

759
    VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
760
    for (int iq = 0; iq < numPoints; iq++) {
761
762
      Lb[iq].set(0.0);
    }
763
764
765
766

    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
767
768
    }
  
769
    for (int iq = 0; iq < numPoints; iq++) {
770
771
      Lb[iq] *= elInfo->getDet();

772
      for (int i = 0; i < nCol; i++) {
773
774
775
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
      }

776
      for (int i = 0; i < nRow; i++) {
777
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
778
779
	for (int j = 0; j < nCol; j++) {
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
	}
      }
    }
    FREE_MEMORY(phival, double, nCol);
  }


  Quad10::Quad10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }


  void Quad10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;
    const double *phi;

798
    if (firstCall) {
799
800
801
802
803
804
805
806
807
808
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int numPoints = quadrature->getNumPoints();

    VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
809
    for (int iq = 0; iq < numPoints; iq++) {
810
811
      Lb[iq].set(0.0);
    }
812

813
814
815
    int myRank = omp_get_thread_num();
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
816
817
    }
  
818
    for (int iq = 0; iq < numPoints; iq++) {
819
820
821
822
823
      Lb[iq] *= elInfo->getDet();

      phi    = phiFast->getPhi(iq);
      grdPsi = psiFast->getGradient(iq);

824
825
826
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPsi)[i]) * phi[j];
827
828
829
830
831
832
833
834
835
836
      }
    }
  }


  Pre10::Pre10(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PSI)
  {
  }

837

838
839
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
840
    VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
841
842
843
    const int *k;
    const double *values;

844
    if (firstCall) {
845
846
847
848
849
850
851
852
853
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q10->getNumberEntries();
854
    int myRank = omp_get_thread_num();
855
    Lb[0].set(0.0);
856
857
858

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
859
860
861
862
    }

    Lb[0] *= elInfo->getDet();

863
864
865
    for (int i = 0; i < nRow; i++) {
      for (int j = 0; j < nCol; j++) {
	k = q10->getKVec(i, j);
866
	values = q10->getValVec(i, j);
867
	double val = 0.0;
Thomas Witkowski's avatar
Thomas Witkowski committed
868
869
870
	for (int m = 0; m < nEntries[i][j]; m++) {
	  val += values[m] * Lb[0][k[m]];
	}
871
872
873
874
875
876
877
878
879
880
	(*mat)[i][j] += val;
      }
    }
  }


  Stand01::Stand01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, false, GRD_PHI)
  {}

881

882
883
884
885
886
887
888
889
  void Stand01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > grdPhi(dim, nCol, NO_INIT);

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

    int numPoints = quadrature->getNumPoints();
890
891
    VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
    int myRank = omp_get_thread_num();
892

893
    for (int iq = 0; iq < numPoints; iq++) {
894
895
      Lb[iq].set(0.0);
    }
896
897
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
898
899
    }
  
900
    for (int iq = 0; iq < numPoints; iq++) {
901
902
      Lb[iq] *= elInfo->getDet();

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

907
      for (int i = 0; i < nRow; i++) {
908
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
909
910
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
911
      }
912
    } 
913
914
915
916
  }

  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
917
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
918
919
920
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    int numPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
921
    int myRank = omp_get_thread_num();
922

923
    for (int iq = 0; iq < numPoints; iq++) {
924
925
      Lb[iq].set(0.0);
    }
926
927
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
928
929
    }
  
930
    for (int iq = 0; iq < numPoints; iq++) {
931
932
      Lb[iq] *= elInfo->getDet();

933
      for (int i = 0; i < nRow; i++) {
934
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
935
936
937
938
939
940
941
942
943
944
945
946
947
948
	(*vec)[i] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi);
      }
    }
  }

  Quad01::Quad01(Operator *op, Assembler *assembler, Quadrature *quad) 
    : FirstOrderAssembler(op, assembler, quad, true, GRD_PHI)
  {
  }

  void Quad01::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
    VectorOfFixVecs<DimVec<double> > *grdPhi;

949
    if (firstCall) {
950
951
952
953
954
955
956
957
958
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_GRD_PHI);
      firstCall = false;
    }

    int numPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
959
    int myRank = omp_get_thread_num();
960

961
    for (int iq = 0; iq < numPoints; iq++) {
962
963
      Lb[iq].set(0.0);
    }
964
965
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
966
967
    }
  
968
    for (int iq = 0; iq < numPoints; iq++) {
969
970
      Lb[iq] *= elInfo->getDet();

971
      const double *psi = psiFast->getPhi(iq);
972
973
      grdPhi = phiFast->getGradient(iq);

974
975
976
      for (int i = 0; i < nRow; i++) {
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * (*grdPhi)[j]) * psi[i];
977
978
979
980
981
982
983
984
      }
    }
  }

  void Quad10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
    VectorOfFixVecs<DimVec<double> > *grdPsi;

985
    if (firstCall) {
986
987
988
989
990
991
992
993
994
      const BasisFunction *basFcts = owner->getRowFESpace()->getBasisFcts();
      psiFast = updateFastQuadrature(psiFast, basFcts, INIT_GRD_PHI);
      basFcts = owner->getColFESpace()->getBasisFcts();
      phiFast = updateFastQuadrature(phiFast, basFcts, INIT_PHI);
      firstCall = false;
    }

    int numPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
995
    int myRank = omp_get_thread_num();
996

997
    for (int iq = 0; iq < numPoints; iq++) {
998
999
      Lb[iq].set(0.0);
    }
1000
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
For faster browsing, not all history is shown. View entire blame