Assembler.cc 45.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
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
    double *c = GET_MEMORY(double, 1);
684

685
    if (firstCall) {
686
687
688
689
690
691
692
693
694
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    c[0] = 0.0;
695
696
697
698
    int myRank = omp_get_thread_num();

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

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

    if (symmetric) {
704
705
706
707
      for (int i = 0; i < nRow; i++) {
	(*mat)[i][i] += c[0] * q00->getValue(i,i);
	for (int j = i + 1; j < nCol; j++) {
	  double val = c[0] * 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
716
717
718
719
720
721
722
	  (*mat)[i][j] += c[0]*q00->getValue(i,j);
    }

    FREE_MEMORY(c, double, 1);
  }

  void Pre0::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
723
    double *c = GET_MEMORY(double, 1);
724

725
    if (firstCall) {
726
727
728
729
730
731
732
733
734
735
      q00 = Q00PsiPhi::provideQ00PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q0 = Q0Psi::provideQ0Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

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

736
    int myRank = omp_get_thread_num();
737
    c[0] = 0.0;
738
    for (termIt = terms[myRank].begin(); termIt != terms[myRank].end(); ++termIt) {
739
740
741
742
743
      (static_cast<ZeroOrderTerm*>( *termIt))->getC(elInfo, 1, c);
    }

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

744
    for (int i = 0; i < nRow; i++)
745
746
747
748
749
750
751
752
753
754
755
756
      (*vec)[i] += c[0] * q0->getValue(i);

    FREE_MEMORY(c, double, 1);
  }

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


  void Stand10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
757
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
758
759
760
761
762
763
764
    double *phival = GET_MEMORY(double, nCol);

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

    int numPoints = quadrature->getNumPoints();

765
    VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
766
    for (int iq = 0; iq < numPoints; iq++) {
767
768
      Lb[iq].set(0.0);
    }
769
770
771
772

    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);
773
774
    }
  
775
    for (int iq = 0; iq < numPoints; iq++) {
776
777
      Lb[iq] *= elInfo->getDet();

778
      for (int i = 0; i < nCol; i++) {
779
780
781
	phival[i] = (*(phi->getPhi(i)))(quadrature->getLambda(iq));
      }

782
      for (int i = 0; i < nRow; i++) {
783
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
784
785
	for (int j = 0; j < nCol; j++) {
	  (*mat)[i][j] += quadrature->getWeight(iq) * (Lb[iq] * grdPsi) * phival[j];
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
	}
      }
    }
    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;

804
    if (firstCall) {
805
806
807
808
809
810
811
812
813
814
      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);
815
    for (int iq = 0; iq < numPoints; iq++) {
816
817
      Lb[iq].set(0.0);
    }
818

819
820
821
    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);
822
823
    }
  
824
    for (int iq = 0; iq < numPoints; iq++) {
825
826
827
828
829
      Lb[iq] *= elInfo->getDet();

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

830
831
832
      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];
833
834
835
836
837
838
839
840
841
842
      }
    }
  }


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

843

844
845
  void Pre10::calculateElementMatrix(const ElInfo *elInfo, ElementMatrix *mat)
  {
846
    VectorOfFixVecs<DimVec<double> > Lb(dim, 1, NO_INIT);
847
848
849
    const int *k;
    const double *values;

850
    if (firstCall) {
851
852
853
854
855
856
857
858
859
      q10 = Q10PsiPhi::provideQ10PsiPhi(owner->getRowFESpace()->getBasisFcts(), 
					owner->getColFESpace()->getBasisFcts(), 
					quadrature);
      q1 = Q1Psi::provideQ1Psi(owner->getRowFESpace()->getBasisFcts(),
			       quadrature);
      firstCall = false;
    }

    const int **nEntries = q10->getNumberEntries();
860
    int myRank = omp_get_thread_num();
861
    Lb[0].set(0.0);
862
863
864

    for (int i = 0; i < static_cast<int>( terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, 1, Lb);
865
866
867
868
    }

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

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


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

887

888
889
890
891
892
893
894
895
  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();
896
897
    VectorOfFixVecs<DimVec<double> > Lb(dim, numPoints, NO_INIT);
    int myRank = omp_get_thread_num();
898

899
    for (int iq = 0; iq < numPoints; iq++) {
900
901
      Lb[iq].set(0.0);
    }
902
903
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
904
905
    }
  
906
    for (int iq = 0; iq < numPoints; iq++) {
907
908
      Lb[iq] *= elInfo->getDet();

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

913
      for (int i = 0; i < nRow; i++) {
914
	double psival = (*(psi->getPhi(i)))(quadrature->getLambda(iq));
915
916
	for (int j = 0; j < nCol; j++)
	  (*mat)[i][j] += quadrature->getWeight(iq) * ((Lb[iq] * psival) * grdPhi[j]);
917
      }
918
    } 
919
920
921
922
  }

  void Stand10::calculateElementVector(const ElInfo *elInfo, ElementVector *vec)
  {
923
    DimVec<double> grdPsi(dim, DEFAULT_VALUE, 0.0);
924
925
926
    const BasisFunction *psi = owner->getRowFESpace()->getBasisFcts();
    int numPoints = quadrature->getNumPoints();
    VectorOfFixVecs<DimVec<double> > Lb(dim,numPoints,NO_INIT);
927
    int myRank = omp_get_thread_num();
928

929
    for (int iq = 0; iq < numPoints; iq++) {
930
931
      Lb[iq].set(0.0);
    }
932
933
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
934
935
    }
  
936
    for (int iq = 0; iq < numPoints; iq++) {
937
938
      Lb[iq] *= elInfo->getDet();

939
      for (int i = 0; i < nRow; i++) {
940
	(*(psi->getGrdPhi(i)))(quadrature->getLambda(iq), grdPsi);
941
942
943
944
945
946
947
948
949
950
951
952
953
954
	(*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;

955
    if (firstCall) {
956
957
958
959
960
961
962
963
964
      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);
965
    int myRank = omp_get_thread_num();
966

967
    for (int iq = 0; iq < numPoints; iq++) {
968
969
      Lb[iq].set(0.0);
    }
970
971
    for (int i = 0; i < static_cast<int>(terms[myRank].size()); i++) {
      (static_cast<FirstOrderTerm*>((terms[myRank][i])))->getLb(elInfo, numPoints, Lb);
972
973
    }
  
974
    for (int iq = 0; iq < numPoints; iq++) {
975
976
      Lb[iq] *= elInfo->getDet();

977
      const double *psi = psiFast->getPhi(iq);
978
979
      grdPhi = phiFast->getGradient(iq);

980
981
982
      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];
983
984
985
986
987
988
989
990
      }
    }
  }

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

991
    if (firstCall) {
992
993
994
995
996
997
998
999
1000
      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);