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
	nBasFcts = getNumberOfDofs(i, degree);
196
	(*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
  int Lagrange::getNumberOfDofs(int dim, int degree)
678
679
  {
    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
    int verticesOfPosition = dimOfPosition + 1;
721

722
723
724
    for (int i = 0; i < verticesOfPosition; i++)
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);   
Thomas Witkowski's avatar
Thomas Witkowski committed
725
    
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
 
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	rvec[0] = 0.0;
836
	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
    
    return(const_cast<const double *>( rvec));
  }

860

861
  const WorldVector<double>* 
862
863
  Lagrange::interpol(const ElInfo *el_info, 
		     int no, 
864
865
866
867
		     const int *b_no,
		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
		     WorldVector<double> *vec)
  {
868
869
    FUNCNAME("*Lagrange::interpol_d()");

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

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

879
880
881
882
883
884
    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));
885
      }
886
887
888
889
      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 {
890
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
891
892
893
894
895
896
897
	  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++) {
898
	el_info->coordToWorld(*(*bary)[i], x);
899
900
901
	rvec[i] = (*f)(x);
      }
    }  
902
903
904
905
906
907
908
909
910
911
912
913
    
    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;

914
    const DegreeOfFreedom **dof = el->getDof();
Thomas Witkowski's avatar
Thomas Witkowski committed
915

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

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

933
934
935
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
936

937
    for (int pos = 0, j = 0; pos <= dim; pos++) {
938
      posIndex = INDEX_OF_DIM(pos, dim);
939
      nrDOFs = admin->getNumberOfDofs(posIndex);
940

Thomas Witkowski's avatar
Thomas Witkowski committed
941
      if (nrDOFs) {
942
	n0 = admin->getNumberOfPreDofs(posIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
943
944
945
946
	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
947
	  const int *indi = orderOfPositionIndices(el, posIndex, i);
948

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

955
956
957
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
958
959
960
961
962
963
964
  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);

965
    const DegreeOfFreedom **dof = el->getDof();
Thomas Witkowski's avatar
Thomas Witkowski committed
966
967
968
969
970
    GeoIndex posIndex;
    int nrDOFs, n0, node0, num = 0;

    for (int pos = 0, j = 0; pos <= dim; pos++) {
      posIndex = INDEX_OF_DIM(pos, dim);
971
      nrDOFs = admin->getNumberOfDofs(posIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
972
973

      if (nrDOFs) {
974
	n0 = admin->getNumberOfPreDofs(posIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
975
976
977
978
979
980
981
982
983
984
985
986
987
	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]]);
	}
      }
    }
  }

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

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

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

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

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

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

1026
1027
1028
      for (int iq = 0; iq < nPoints; iq++) {
	el_info->coordToWorld(quad->getLambda(iq), x);
	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
1029
1030
      }
      
1031
1032
1033
1034
      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
1035

1036
	(*fh)[dof[j]] += val;
1037
1038
      }

1039
1040
1041
      el_info = stack.traverseNext(el_info);
    }

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


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

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

1052
1053
1054
    if (n < 1) 
      return;

1055
    int n0 = drv->getFeSpace()->getAdmin()->getNumberOfPreDofs(CENTER);
1056
    Element *el = list->getElement(0);
1057
    int node = drv->getFeSpace()->getMesh()->getNode(CENTER);
1058
    // Parent center
1059
    DegreeOfFreedom dof0 = el->getDof(node, n0);           
1060
    // Newest vertex is center 
Thomas Witkowski's avatar