Lagrange.cc 169 KB
Newer Older
1
2
3
4
#include <stdio.h>
#include <algorithm>
#include <list>
#include "boost/lexical_cast.hpp"
5
6
7
8
9
10
11
12
13
14
15
16
17
18
#include "Mesh.h"
#include "Element.h"
#include "Lagrange.h"
#include "DOFAdmin.h"
#include "RCNeighbourList.h"
#include "DOFVector.h"
#include "Traverse.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "Parametric.h"

namespace AMDiS {

19
20
  using boost::lexical_cast;

21
22
23
24
25
26
  std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  int Lagrange::nBasFctsDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<BasFctType*> Lagrange::phiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<GrdBasFctType*> Lagrange::grdPhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
  std::vector<D2BasFctType*> Lagrange::D2PhiDimDegree[MAX_DIM + 1][MAX_DEGREE + 1];
27
  std::list<Lagrange*> Lagrange::allBasFcts;
28
29
30


  Lagrange::Lagrange(int dim, int degree)
31
    : BasisFunction(std::string("Lagrange"), dim, degree)
32
33
  {
    // set name
34
    name += lexical_cast<std::string>(dim) + " " + lexical_cast<std::string>(degree);
35
36
37
38
39
40
41
42
43
44
45
46
47

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();
  }

  Lagrange::~Lagrange()
  {
48
    for (int i = 0; i < static_cast<int>(bary->size()); i++)
49
      if ((*bary)[i]) {
50
	delete (*bary)[i];
51
52
	(*bary)[i] = NULL;
      }
53
54
55
56
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
57
    std::list<Lagrange*>::iterator it;
Thomas Witkowski's avatar
Thomas Witkowski committed
58
59
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++)
      if (((*it)->dim == dim) && ((*it)->degree == degree))
60
61
	return (*it);

62
    Lagrange* newLagrange = new Lagrange(dim, degree);
63
64
65
66
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

67
68
  void Lagrange::clear()
  {    
69
70
    for (std::list<Lagrange*>::iterator it = allBasFcts.begin(); 
	 it != allBasFcts.end(); it++)
71
      if (*it) {
72
	delete *it;
73
74
	*it = NULL;
      }     
75
76
  }

77
78
  void Lagrange::setFunctionPointer()
  {
79
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
80
      // for all positions
81
      for (int i = 0; i < dim+1; i++) {
82
	// no vertex dofs for degree 0 ?
83
84
	if (degree == 0 && i != dim) 
	  continue;
85
	// for all vertices/edges/...
86
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
87
	  // for all dofs at this position
88
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
89
	    // basis function
90
	    phiDimDegree[dim][degree].push_back(new Phi(this, 
91
							INDEX_OF_DIM(i ,dim),j, k));
92
	    // gradients
93
	    grdPhiDimDegree[dim][degree].push_back(new GrdPhi(this, 
94
95
							      INDEX_OF_DIM(i, dim),
							      j, k));
96
	    // D2-Matrices
97
	    D2PhiDimDegree[dim][degree].push_back(new D2Phi(this, 
98
99
							    INDEX_OF_DIM(i, dim),
							    j, k));
100
101
102
103
	  }
	}
      }
    }
104
    phi = &phiDimDegree[dim][degree];
105
    grdPhi = &grdPhiDimDegree[dim][degree];
106
    d2Phi = &D2PhiDimDegree[dim][degree];
107

Thomas Witkowski's avatar
Thomas Witkowski committed
108
    switch (degree) {
109
110
111
112
113
114
115
116
117
118
119
    case 0:
      refineInter_fct = refineInter0;
      coarseRestr_fct = coarseRestr0;
      coarseInter_fct = coarseInter0;  // not yet implemented
      break;
    case 1:
      refineInter_fct = refineInter1;
      coarseRestr_fct = coarseRestr1;
      coarseInter_fct = NULL;  // not yet implemented
      break;
    case 2:
Thomas Witkowski's avatar
Thomas Witkowski committed
120
      switch (dim) {
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
      case 1:
	refineInter_fct = refineInter2_1d;
	coarseRestr_fct = NULL; // not yet implemented
	coarseInter_fct = coarseInter2_1d;
	break;
      case 2: 
	refineInter_fct = refineInter2_2d;
	coarseRestr_fct = coarseRestr2_2d;
	coarseInter_fct = coarseInter2_2d;
	break;
      case 3:
	refineInter_fct = refineInter2_3d;
	coarseRestr_fct = coarseRestr2_3d;
	coarseInter_fct = coarseInter2_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 3:
Thomas Witkowski's avatar
Thomas Witkowski committed
140
      switch (dim) {
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
      case 1:
	refineInter_fct = refineInter3_1d;
	coarseRestr_fct = coarseRestr3_1d;
	coarseInter_fct = coarseInter3_1d;
	break;
      case 2: 
	refineInter_fct = refineInter3_2d;
	coarseRestr_fct = coarseRestr3_2d;
	coarseInter_fct = coarseInter3_2d;
	break;
      case 3:
	refineInter_fct = refineInter3_3d;
	coarseRestr_fct = coarseRestr3_3d;
	coarseInter_fct = coarseInter3_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    case 4:
Thomas Witkowski's avatar
Thomas Witkowski committed
160
      switch (dim) {
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
      case 1: 
	refineInter_fct = refineInter4_1d;
	coarseRestr_fct = coarseRestr4_1d;
	coarseInter_fct = coarseInter4_1d;
	break;
      case 2:
	refineInter_fct = refineInter4_2d;
	coarseRestr_fct = coarseRestr4_2d;
	coarseInter_fct = coarseInter4_2d;
	break;
      case 3:
	refineInter_fct = refineInter4_3d;
	coarseRestr_fct = coarseRestr4_3d;
	coarseInter_fct = coarseInter4_3d;
	break;
      default: ERROR_EXIT("invalid dim\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

  void Lagrange::setNDOF()
  {
186
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
187
      ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);
188

Thomas Witkowski's avatar
Thomas Witkowski committed
189
      if (degree != 0)
190
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
191
192
      else
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;      
193

194
      for (int i = 1; i < dim + 1; i++) {
195
196
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
197
	for (int j = 0; j < i; j++) {
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
	  (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] -= 
	    Global::getGeo(INDEX_OF_DIM(j, dim), i) * 
	    (*ndofDimDegree[dim][degree])[INDEX_OF_DIM(j, dim)];
	}
      }
      nBasFctsDimDegree[dim][degree] = nBasFcts;
    }
    nBasFcts = nBasFctsDimDegree[dim][degree];
    nDOF = ndofDimDegree[dim][degree];
  } 

  DimVec<double> *Lagrange::getCoords(int i) const
  {
    return (*bary)[i];
  }

  void Lagrange::setVertices(int dim, int degree, 
			     GeoIndex position, int positionIndex, int nodeIndex, 
			     int** vertices)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
218
219
    FUNCNAME("Lagrange::setVertices()");

Thomas Witkowski's avatar
Thomas Witkowski committed
220
    TEST_EXIT_DBG(*vertices == NULL)("vertices != NULL\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
221

222
223
    int dimOfPosition = DIM_OF_INDEX(position, dim);

Thomas Witkowski's avatar
Thomas Witkowski committed
224
    *vertices = new int[dimOfPosition + 1];
225

Thomas Witkowski's avatar
Thomas Witkowski committed
226
    if (degree == 4 && dimOfPosition == 1) {
227
228
229
230
231
232
233
234
      (*vertices)[(nodeIndex != 2) ? 0 : 1] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex,
							      0);
      (*vertices)[(nodeIndex != 2) ? 1 : 0] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex, 
							      1);
Thomas Witkowski's avatar
Thomas Witkowski committed
235
    } else if (degree==4 && dimOfPosition == 2) {
236
      for (int i = 0; i < dimOfPosition + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
237
	(*vertices)[(i + dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
238
239
240
241
242
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
243
      for (int i = 0; i < dimOfPosition + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
244
	(*vertices)[(i + nodeIndex) % (dimOfPosition + 1)] = 
245
246
247
248
249
250
251
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

252
  Lagrange::Phi::Phi(Lagrange* owner, 
253
254
255
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
256
    : vertices(NULL)
257
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
258
259
260
261
262
263
264
265
266
    FUNCNAME("Lagrange::Phi::Phi()");

    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);
267
268
  
    // set function pointer
269
    switch (owner->getDegree()) {
270
    case 0:
271
      switch (position) {
272
273
274
275
276
277
278
279
      case CENTER:
	func = phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
280
      switch (position) {
281
282
283
284
285
286
287
288
      case VERTEX:
	func = phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
289
      switch (position) {
290
291
292
293
      case VERTEX:
	func = phi2v;
	break;
      case EDGE:
294
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
295
296
297
	func = phi2e;
	break;
      case CENTER:
298
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
299
300
301
302
303
304
305
	func = phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
306
      switch (position) {
307
308
309
310
311
312
313
      case VERTEX:
	func = phi3v;
	break;
      case EDGE:
	func = phi3e;
	break;
      case FACE:
314
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
315
316
317
	func = phi3f;
	break;
      case CENTER:
318
	switch (owner->getDim()) {
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
	case 1:
	  func = phi3e;
	  break;
	case 2:
	  func = phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
335
      switch (position) {
336
337
338
339
      case VERTEX:
	func = phi4v;
	break;
      case EDGE:
340
	if (nodeIndex == 1)
341
342
343
344
345
	  func = phi4e1;
	else 
	  func = phi4e0;
	break;
      case FACE:
346
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
347
348
349
	func = phi4f;      
	break;
      case CENTER:
350
	switch (owner->getDim()) {
351
	case 1:
Thomas Witkowski's avatar
Thomas Witkowski committed
352
	  if (nodeIndex == 1)
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
	    func = phi4e1;
	  else 
	    func = phi4e0;
	  break;
	case 2:
	  func = phi4f;
	  break;
	case 3:
	  func = phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

377
  Lagrange::Phi::~Phi() { 
378
    delete [] vertices; 
379
  }
380
381


382
  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
383
384
385
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
386
    : vertices(NULL)
387
388
389
390
391
392
393
394
395
396
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
397
    switch (owner->getDegree()) {
398
    case 0:
399
      switch (position) {
400
401
402
403
404
405
406
407
      case CENTER:
	func = grdPhi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
408
      switch (position) {
409
410
411
412
413
414
415
416
      case VERTEX:
	func = grdPhi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
417
      switch (position) {
418
419
420
421
422
423
424
      case VERTEX:
	func = grdPhi2v;
	break;
      case EDGE:
	func = grdPhi2e;
	break;
      case CENTER:
425
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
426
427
428
429
430
431
432
	func = grdPhi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
433
      switch (position) {
434
435
436
437
438
439
440
      case VERTEX:
	func = grdPhi3v;
	break;
      case EDGE:
	func = grdPhi3e;
	break;
      case FACE:
441
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
442
443
444
	func = grdPhi3f;
	break;
      case CENTER:
445
	switch (owner->getDim()) {
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
	case 1:
	  func = grdPhi3e;
	  break;
	case 2:
	  func = grdPhi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
462
      switch (position) {
463
464
465
466
      case VERTEX:
	func = grdPhi4v;
	break;
      case EDGE:
467
	if (nodeIndex == 1)
468
469
470
471
472
	  func = grdPhi4e1;
	else 
	  func = grdPhi4e0;
	break;
      case FACE:
473
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
474
475
476
	func = grdPhi4f;      
	break;
      case CENTER:
477
	switch (owner->getDim()) {
478
	case 1:
479
	  if (nodeIndex == 1)
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
	    func = grdPhi4e1;
	  else 
	    func = grdPhi4e0;
	  break;
	case 2:
	  func = grdPhi4f;
	  break;
	case 3:
	  func = grdPhi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
504
505
  Lagrange::GrdPhi::~GrdPhi() 
  { 
506
    delete [] vertices; 
507
  }
508

509
  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
510
511
512
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
513
    : vertices(NULL)
514
515
516
517
518
519
520
521
522
523
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
524
    switch (owner->getDegree()) {
525
    case 0:
526
      switch (position) {
527
528
529
530
531
532
533
534
      case CENTER:
	func = D2Phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
535
      switch (position) {
536
537
538
539
540
541
542
543
      case VERTEX:
	func = D2Phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
544
      switch (position) {
545
546
547
548
      case VERTEX:
	func = D2Phi2v;
	break;
      case EDGE:
549
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
550
551
552
	func = D2Phi2e;
	break;
      case CENTER:
553
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
554
555
556
557
558
559
560
	func = D2Phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
561
      switch (position) {
562
563
564
565
566
567
568
      case VERTEX:
	func = D2Phi3v;
	break;
      case EDGE:
	func = D2Phi3e;
	break;
      case FACE:
569
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
570
571
572
	func = D2Phi3f;
	break;
      case CENTER:
573
	switch (owner->getDim()) {
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
	case 1:
	  func = D2Phi3e;
	  break;
	case 2:
	  func = D2Phi3f;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 4:
590
      switch (position) {
591
592
593
594
      case VERTEX:
	func = D2Phi4v;
	break;
      case EDGE:
Thomas Witkowski's avatar
Thomas Witkowski committed
595
	if (nodeIndex == 1)
596
597
598
599
600
	  func = D2Phi4e1;
	else 
	  func = D2Phi4e0;
	break;
      case FACE:
601
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
602
603
604
	func = D2Phi4f;      
	break;
      case CENTER:
605
	switch (owner->getDim()) {
606
	case 1:
607
	  if (nodeIndex == 1)
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
	    func = D2Phi4e1;
	  else 
	    func = D2Phi4e0;
	  break;
	case 2:
	  func = D2Phi4f;
	  break;
	case 3:
	  func = D2Phi4c;
	  break;
	default:
	  ERROR_EXIT("invalid dim\n");
	  break;
	}
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    default:
      ERROR_EXIT("invalid degree\n");
    }
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
632
633
  Lagrange::D2Phi::~D2Phi() 
  { 
634
    delete [] vertices; 
635
  }
636

Thomas Witkowski's avatar
Thomas Witkowski committed
637
  void Lagrange::createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
638
639
			      DimVec<double>* vec)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
640
    if (vec == NULL)
641
      vec = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
642

643
644
    if (dimIndex == numCoords - 1) {
      (*vec)[coordInd[dimIndex]] = double(rest) / degree;
645
      DimVec<double>* newCoords = new DimVec<double>(*vec);
646
647
      bary->push_back(newCoords);
    } else {
648
      for (int i = rest - 1; i >= 1; i--) {
649
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
650
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
651
652
653
654
655
656
      }
    }
  }

  void Lagrange::setBary()
  {
657
    bary = &baryDimDegree[dim][degree];
658
659
660

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
661
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
662
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
Thomas Witkowski's avatar
Thomas Witkowski committed
663
	  int *coordInd = new int[i + 1];      // indices of relevant coords
Thomas Witkowski's avatar
Thomas Witkowski committed
664
	  for (int k = 0; k < i + 1; k++)
665
	    coordInd[k] = Global::getReferenceElement(dim)->
Thomas Witkowski's avatar
Thomas Witkowski committed
666
	      getVertexOfPosition(INDEX_OF_DIM(i, dim), j, k);
Thomas Witkowski's avatar
Thomas Witkowski committed
667
	  
668
	  createCoords(coordInd, i + 1, 0, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
669
670
671
	  delete [] coordInd;
	  if (static_cast<int>(bary->size()) == nBasFcts) 
	    return;
672
673
674
675
676
677
678
679
	}
      }
    }
  }

  int Lagrange::getNumberOfDOFs(int dim, int degree)
  {
    int result = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
680
    for (int i = 0; i <= degree; i++)
681
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
Thomas Witkowski's avatar
Thomas Witkowski committed
682

683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
    return result;
  }


  const DegreeOfFreedom *Lagrange::getDOFIndices(const Element *el, 
						 const DOFAdmin& admin,
						 DegreeOfFreedom *idof) const
  {
    WARNING("depricated: use getLocalIndices() instead\n");
    return getLocalIndices(el, &admin, idof);
  }

  int* Lagrange::orderOfPositionIndices(const Element* el, 
					GeoIndex position, 
					int positionIndex) const
  {
    static int sortedVertex = 0;
    static int sortedEdgeDeg2 = 0;
    static int sortedEdgeDeg3[2][2] = {{0, 1}, {1, 0}};
    static int sortedEdgeDeg4[2][3] = {{0, 1, 2}, {2, 1, 0}};
    static int sortedFaceDeg3 = 0;
    static int sortedFaceDeg4[7][3] = {{0,0,0}, {0,2,1}, {1,0,2}, {0,1,2},
				       {2,1,0}, {1,2,0}, {2,0,1}};
    static int sortedCenterDeg4 = 0;

    int dimOfPosition = DIM_OF_INDEX(position, dim);

    // vertex
Thomas Witkowski's avatar
Thomas Witkowski committed
711
712
    if (dimOfPosition == 0)
      return &sortedVertex;   
713
714

    // edge
715
    if (dimOfPosition == 1 && degree == 2)
Thomas Witkowski's avatar
Thomas Witkowski committed
716
      return &sortedEdgeDeg2;    
717
718
      
    int vertex[3];
719
    int** dof = const_cast<int**>(el->getDOF());
720
721
    int verticesOfPosition = dimOfPosition + 1;

Thomas Witkowski's avatar
Thomas Witkowski committed
722
    for (int i = 0; i < verticesOfPosition; i++)
723
      vertex[i] = Global::getReferenceElement(dim)->
Thomas Witkowski's avatar
Thomas Witkowski committed
724
725
	getVertexOfPosition(position, positionIndex, i);   
    
726
727
    if (dimOfPosition == 1) {
      if (degree == 3) {
Thomas Witkowski's avatar
Thomas Witkowski committed
728
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
729
	  return sortedEdgeDeg3[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
730
731
	else
	  return sortedEdgeDeg3[1];	
732
      } else { // degree == 4
Thomas Witkowski's avatar
Thomas Witkowski committed
733
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
734
	  return sortedEdgeDeg4[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
735
736
	else
	  return sortedEdgeDeg4[1];	
737
738
739
740
      }
    }

    // face
741
742
    if (dimOfPosition == 2) {
      if (degree == 3) {
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
	return &sortedFaceDeg3;
      } else {  // degree == 4!
	int no = 0;
	const Element *refElem = Global::getReferenceElement(dim);

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0]) 
	  no++;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 1)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0])
	  no += 2;

	if (dof[refElem->getVertexOfPosition(position, positionIndex, 2)][0] <
	    dof[refElem->getVertexOfPosition(position, positionIndex, 0)][0])
	  no += 4;

	return sortedFaceDeg4[no];
      }
    }

    // center
765
    if (dimOfPosition == 3 && degree == 4)
Thomas Witkowski's avatar
Thomas Witkowski committed
766
      return &sortedCenterDeg4;    
767
768
769
770
771

    ERROR_EXIT("should not be reached\n");
    return NULL;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
772
  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
773
774
775
776
  {
    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
777
    int index = 0;
778
779
    int offset = 0;
    BoundaryType boundaryType;
780
    for (int i = dim - 1; i > 0; i--)
781
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
782

783
    for (int i = 0; i < dim; i++) {
784
785
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
786
	boundaryType = el_info->getBoundary(j);
787
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
Thomas Witkowski's avatar
Thomas Witkowski committed
788
789
	for (int k = 0; k < kto; k++)
	  bound[index++] = boundaryType;	
790
      }
791
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
792
793
794
    }

    // interior nodes in the center
Thomas Witkowski's avatar
Thomas Witkowski committed
795
796
    for (int i = 0; i < (*nDOF)[CENTER]; i++)
      bound[index++] = INTERIOR;    
797

798
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
799
800
801
802
803
804
805
806
  }


  const double* Lagrange::interpol(const ElInfo *el_info, 
				   int no, const int *b_no, 
				   AbstractFunction<double, WorldVector<double> > *f, 
				   double *vec)
  {
807
808
    FUNCNAME("Lagrange::interpol()");

809
810
    static double* localVec = NULL;
    static int localVecSize = 0;
811

812
    double *rvec = NULL;
813

814
    if (vec) {
815
816
      rvec = vec;
    } else {
817
818
819
820
821
822
823
824
825
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new double[nBasFcts]; 
      }
      if (!localVec)
	localVec = new double[nBasFcts]; 

      localVecSize = nBasFcts;
      rvec = localVec;
826
827
    } 

828
    WorldVector<double> x;
829
830

    el_info->testFlag(Mesh::FILL_COORDS);
831
832
833
834
835
836
 
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	rvec[0] = 0.0;
	return(const_cast<const double *>( rvec));
837
      }
838
839
840
841
842
	
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
843
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
844
845
846
847
848
849
850
851
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      int vertices = Global::getGeo(VERTEX, dim);
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
852
	el_info->coordToWorld(*(*bary)[i], x);
853
854
855
	rvec[i] = (*f)(x);
      }
    }  
856
857
858
859
860
861
862
863
864
865
    
    return(const_cast<const double *>( rvec));
  }

  const WorldVector<double>* 
  Lagrange::interpol(const ElInfo *el_info, int no, 
		     const int *b_no,
		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
		     WorldVector<double> *vec)
  {
866
867
    FUNCNAME("*Lagrange::interpol_d()");

868
    static WorldVector<double> *inter;
869
    WorldVector<double> *rvec = 
870
      vec ? vec : (inter ? inter : inter = new WorldVector<double>[getNumber()]);
871
    WorldVector<double> x;
872
873
874
875
876

    el_info->testFlag(Mesh::FILL_COORDS);
  
    int vertices = Global::getGeo(VERTEX, dim);

877
878
879
880
881
882
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	for (int i = 0; i < dow; i++)
	  rvec[0][i] = 0.0;
	return(const_cast<const WorldVector<double> *>( rvec));
883
      }
884
885
886
887
      for (int i = 0; i < no; i++) {
	if (b_no[i] < Global::getGeo(VERTEX, dim)) {
	  rvec[i] = (*f)(el_info->getCoord(b_no[i]));
	} else {
888
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
889
890
891
892
893
894
895
	  rvec[i] = (*f)(x);
	}
      }
    } else {
      for (int i = 0; i < vertices; i++)
	rvec[i] = (*f)(el_info->getCoord(i));
      for (int i = vertices; i < nBasFcts; i++) {
896
	el_info->coordToWorld(*(*bary)[i], x);
897
898
899
	rvec[i] = (*f)(x);
      }
    }  
900
901
902
903
904
905
906
907
908
909
910
911
912
    
    return(const_cast<const WorldVector<double> *>( rvec));
  }


  const DegreeOfFreedom *Lagrange::getLocalIndices(const Element* el,
						   const DOFAdmin *admin,
						   DegreeOfFreedom *indices) const
  {
    static DegreeOfFreedom *localVec = NULL;
    static int localVecSize = 0;

    const DegreeOfFreedom **dof = el->getDOF();
Thomas Witkowski's avatar
Thomas Witkowski committed
913

914
    int nrDOFs, n0, node0, num = 0;
915
916
    GeoIndex posIndex;
    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
917

918
    if (indices) {
919
920
      result = indices;
    } else {
921
922
923
#ifdef _OPENMP
      ERROR_EXIT("Using static variable while using OpenMP parallelization!\n");
#endif
Thomas Witkowski's avatar
Thomas Witkowski committed
924
      if (localVec && nBasFcts > localVecSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
925
	delete [] localVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
926
927
928
929
930
	localVec = new DegreeOfFreedom[nBasFcts];
      }
      if (!localVec)
	localVec = new DegreeOfFreedom[nBasFcts];

931
932
933
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
934

935
    for (int pos = 0, j = 0; pos <= dim; pos++) {
936
937
938
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

Thomas Witkowski's avatar
Thomas Witkowski committed
939
940
941
942
943
944
      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
945
	  const int *indi = orderOfPositionIndices(el, posIndex, i);
946

Thomas Witkowski's avatar
Thomas Witkowski committed
947
	  for (int k = 0; k < nrDOFs; k++)
Thomas Witkowski's avatar
Thomas Witkowski committed
948
	    result[j++] = dof[node0][n0 + indi[k]];  
949
950
951
	}
      }
    }
952

953
954
955
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
  void Lagrange::getLocalDofPtrVec(const Element *el, 
				   const DOFAdmin *admin,
				   std::vector<const DegreeOfFreedom*>& vec) const
  {
    if (static_cast<int>(vec.size()) < nBasFcts)
      vec.resize(nBasFcts);

    const DegreeOfFreedom **dof = el->getDOF();
    GeoIndex posIndex;
    int nrDOFs, n0, node0, num = 0;

    for (int pos = 0, j = 0; pos <= dim; pos++) {
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

	for (int i = 0; i < num; node0++, i++) {
	  const int *indi = orderOfPositionIndices(el, posIndex, i);

	  for (int k = 0; k < nrDOFs; k++)
	    vec[j++] = &(dof[node0][n0 + indi[k]]);
	}
      }
    }
  }

986
987
988
989
990
991
992
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

993
  void Lagrange::l2ScpFctBas(Quadrature *quad,
994
995
996
			     AbstractFunction<double, WorldVector<double> >* f,
			     DOFVector<double>* fh)
  {
997
998
    FUNCNAME("Lagrange::l2ScpFctBas()");

999
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1000
1001
1002
1003
    TEST_EXIT_DBG(fh->getFESpace())
      ("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
      ("wrong basis fcts for fh\n");
1004
1005
    TEST_EXIT_DBG(!fh->getFESpace()->getMesh()->getParametric())
      ("Not yet implemented!");
1006

1007
1008
    if (!quad)
      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);
1009

Thomas Witkowski's avatar
Thomas Witkowski committed
1010
    const FastQuadrature *quad_fast = 
1011
1012
1013
1014
1015
      FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
    double *wdetf_qp = new double[quad->getNumPoints()];
    int nPoints = quad->getNumPoints();
    DOFAdmin *admin = fh->getFESpace()->getAdmin();
    WorldVector<double> x;
1016

1017
1018
    TraverseStack stack;
    ElInfo *el_info = stack.traverseFirst(fh->getFESpace()->getMesh(), -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1019
					  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1020
    while (el_info) {
1021
1022
      const DegreeOfFreedom  *dof = getLocalIndices(el_info->getElement(), admin, NULL);
      double det = el_info->getDet();
1023

1024
1025
1026
      for (int iq = 0; iq < nPoints; iq++) {
	el_info->coordToWorld(quad->getLambda(iq), x);
	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
1027
1028
      }
      
1029
1030
1031
1032
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int iq = 0; iq < nPoints; iq++)
	  val += quad_fast->getPhi(iq, j) * wdetf_qp[iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
1033

1034
	(*fh)[dof[j]] += val;
1035
1036
      }

1037
1038
1039
      el_info = stack.traverseNext(el_info);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1040
    delete [] wdetf_qp;
1041
1042
1043
1044
1045
1046
1047
  }


  // ===== refineInter functions =====

  void  Lagrange::refineInter0(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
1048
    FUNCNAME("Lagrange::refineInter0()");
1049

1050
1051
1052
1053
1054
1055
    if (n < 1) 
      return;

    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
    Element *el = list->getElement(0);
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);
1056
1057
1058
1059
    // Parent center
    DegreeOfFreedom dof0 = el->getDOF(node, n0);           
    // Newest vertex is center 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);  
1060
1061

    (*drv)[dof_new] = (*drv)[dof0];
1062
1063
    // Newest vertex is center 
    dof_new = el->getChild(1)->getDOF(node, n0);  
1064
1065
1066
    (*drv)[dof_new] = (*drv)[dof0];
  }

1067
1068
  void Lagrange::refineInter1(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
			      BasisFunction* basFct)
1069
  {
1070
    FUNCNAME("Lagrange::refineInter1_1d()");
1071

1072
1073
    if (n < 1) 
      return;
1074

1075
1076
1077
    int dim = drv->getFESpace()->getMesh()->getDim();
    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    Element *el = list->getElement(0);
1078
1079
1080
1081
1082
1083
    // 1st endpoint of refinement edge
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    // 2nd endpoint of refinement edge          
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
    // newest vertex is DIM 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0); 
1084
    (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
1085
1086
  }

1087
1088
  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
				 BasisFunction* basFct)
1089
  {
1090
    FUNCNAME("Lagrange::refineInter2_1d()");
1091

1092
1093
1094
1095
1096
1097
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1098
  
1099
1100
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1101
1102
1103
1104
1105

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

1106
1107
    // newest vertex is DIM
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node + 1, n0);
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
    (*drv)[cdof] = (*drv)[pdof[2]];

    node = drv->getFESpace()->getMesh()->getNode(CENTER);        
    n0 = admin->getNumberOfPreDOFs(CENTER);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
1119
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1120
1121
1122
1123
1124
1125
1126

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1127
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1128
1129
1130
1131
1132
1133
  }

  void Lagrange::refineInter2_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1134
    FUNCNAME("Lagrange::refineInter2_2d()");
1135

1136
1137
    if (n < 1) 
      return;
1138
1139

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1140
1141
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1142
  
1143
1144
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1145
1146
1147
1148
1149

    /****************************************************************************/
    /*  newest vertex of child[0] and child[1]                                  */
    /****************************************************************************/

1150
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
    (*drv)[cdof] = (*drv)[pdof[5]];

    node = drv->getFESpace()->getMesh()->getNode(EDGE);        
    n0 = admin->getNumberOfPreDOFs(EDGE);

    /****************************************************************************/
    /*  midpoint of edge on child[0] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(0)->getDOF(node, n0); 
    (*drv)[cdof] = 
1162
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1163
1164
1165
1166
1167
1168
1169

    /****************************************************************************/
    /* node in the common edge of child[0] and child[1]                         */
    /****************************************************************************/

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1170
1171
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1172
1173
1174
1175
1176
1177
1178

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/<