Lagrange.cc 179 KB
Newer Older
1
2
3
4
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"

#include <stdio.h>
#include <algorithm>
#include <list>

namespace AMDiS {

19
  std::vector<DimVec<double>* > Lagrange::baryDimDegree[MAX_DIM+1][MAX_DEGREE+1];
20
21
  DimVec<int>* Lagrange::ndofDimDegree[MAX_DIM+1][MAX_DEGREE+1];
  int Lagrange::nBasFctsDimDegree[MAX_DIM+1][MAX_DEGREE+1];
22
23
24
25
  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];
  std::list<Lagrange*> Lagrange::allBasFcts;
26
27
28


  Lagrange::Lagrange(int dim, int degree)
29
    : BasisFunction(std::string("Lagrange"), dim, degree)
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
  {
    // set name
    char dummy[3];
    sprintf(dummy, "%d", dim);
    name.append(dummy);
    name.append(" ");
    sprintf(dummy, "%d", degree);  
    name.append(dummy);

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();

  }

  Lagrange::~Lagrange()
  {
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
56
    std::list<Lagrange*>::iterator it;
57
58
59
60
61
62
63
64
65
66
67
68
69
70

    for(it = allBasFcts.begin(); it != allBasFcts.end(); it++) {
      if(((*it)->dim == dim) && ((*it)->degree == degree)) {
	return (*it);
      } 
    }

    Lagrange* newLagrange = NEW Lagrange(dim, degree);
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

  void Lagrange::setFunctionPointer()
  {
71
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
72
      // for all positions
73
      for (int i = 0; i < dim+1; i++) {
74
	// no vertex dofs for degree 0 ?
75
76
	if (degree == 0 && i != dim) 
	  continue;
77
	// for all vertices/edges/...
78
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
79
	  // for all dofs at this position
80
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
81
82
	    // basis function
	    phiDimDegree[dim][degree].push_back(NEW Phi(this, 
83
							INDEX_OF_DIM(i ,dim),j, k));
84
85
	    // gradients
	    grdPhiDimDegree[dim][degree].push_back(NEW GrdPhi(this, 
86
87
							      INDEX_OF_DIM(i, dim),
							      j, k));
88
89
	    // D2-Matrices
	    D2PhiDimDegree[dim][degree].push_back(NEW D2Phi(this, 
90
91
							    INDEX_OF_DIM(i, dim),
							    j, k));
92
93
94
95
	  }
	}
      }
    }
96
    phi = &phiDimDegree[dim][degree];
97
    grdPhi = &grdPhiDimDegree[dim][degree];
98
    d2Phi = &D2PhiDimDegree[dim][degree];
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177

    switch(degree) {
    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:
      switch(dim) {
      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:
      switch(dim) {
      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:
      switch(dim) {
      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()
  {
178
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
179
180
      ndofDimDegree[dim][degree] = NEW DimVec<int>(dim, DEFAULT_VALUE, 0);

181
182
183
184
185
      if (degree != 0) {
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
      } else {
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;
      }
186

187
      for (int i = 1; i < dim + 1; i++) {
188
189
	nBasFcts = getNumberOfDOFs(i, degree);
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
190
	for (int j = 0; j < i; j++) {
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
	  (*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)
  {
211
    TEST_EXIT_DBG((*vertices) == NULL)("vertices != NULL\n");
212
213
214
215
    int dimOfPosition = DIM_OF_INDEX(position, dim);

    *vertices = GET_MEMORY(int, dimOfPosition + 1);

216
    if ((degree == 4) && (dimOfPosition==1)) {
217
218
219
220
221
222
223
224
225
      (*vertices)[(nodeIndex != 2) ? 0 : 1] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex,
							      0);
      (*vertices)[(nodeIndex != 2) ? 1 : 0] = 
	Global::getReferenceElement(dim)->getVertexOfPosition(position,
							      positionIndex, 
							      1);
    } else if ((degree==4) && (dimOfPosition==2)) {
226
227
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
228
229
230
231
232
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
233
234
      for (int i = 0; i < dimOfPosition + 1; i++) {
	(*vertices)[(i+nodeIndex) % (dimOfPosition + 1)] = 
235
236
237
238
239
240
241
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

242
  Lagrange::Phi::Phi(Lagrange* owner, 
243
244
245
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
246
    : vertices(NULL)
247
248
249
250
251
252
253
254
255
256
257
  {
    FUNCNAME("Lagrange::Phi::Phi")
      // get relevant vertices
      Lagrange::setVertices(owner->getDim(), 
			    owner->getDegree(), 
			    position, 
			    positionIndex, 
			    nodeIndex, 
			    &vertices);
  
    // set function pointer
258
    switch (owner->getDegree()) {
259
    case 0:
260
      switch (position) {
261
262
263
264
265
266
267
268
      case CENTER:
	func = phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
269
      switch (position) {
270
271
272
273
274
275
276
277
      case VERTEX:
	func = phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
278
      switch (position) {
279
280
281
282
      case VERTEX:
	func = phi2v;
	break;
      case EDGE:
283
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
284
285
286
	func = phi2e;
	break;
      case CENTER:
287
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
288
289
290
291
292
293
294
	func = phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
295
      switch (position) {
296
297
298
299
300
301
302
      case VERTEX:
	func = phi3v;
	break;
      case EDGE:
	func = phi3e;
	break;
      case FACE:
303
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
304
305
306
	func = phi3f;
	break;
      case CENTER:
307
	switch (owner->getDim()) {
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
	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:
324
      switch (position) {
325
326
327
328
      case VERTEX:
	func = phi4v;
	break;
      case EDGE:
329
	if (nodeIndex == 1)
330
331
332
333
334
	  func = phi4e1;
	else 
	  func = phi4e0;
	break;
      case FACE:
335
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
336
337
338
	func = phi4f;      
	break;
      case CENTER:
339
	switch (owner->getDim()) {
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
	case 1:
	  if(nodeIndex == 1)
	    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");
    }
  }

366
367
368
  Lagrange::Phi::~Phi() { 
    DELETE [] vertices; 
  }
369
370


371
  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
372
373
374
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
375
    : vertices(NULL)
376
377
378
379
380
381
382
383
384
385
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

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

493
494
495
  Lagrange::GrdPhi::~GrdPhi() { 
    DELETE [] vertices; 
  }
496

497
  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
498
499
500
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
501
    : vertices(NULL)
502
503
504
505
506
507
508
509
510
511
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

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

620
621
622
  Lagrange::D2Phi::~D2Phi() { 
    DELETE [] vertices; 
  }
623
624
625
626
627
628
629

  void Lagrange::createCoords(int* coordInd, 
			      int numCoords, 
			      int dimIndex, 
			      int rest, 
			      DimVec<double>* vec)
  {
630
631
632
    if (vec == NULL) {
      vec = NEW DimVec<double>(dim, DEFAULT_VALUE, 0.0);
    }
633

634
    if (dimIndex == numCoords-1) {
635
636
637
638
      (*vec)[coordInd[dimIndex]] = double(rest)/degree;
      DimVec<double>* newCoords = NEW DimVec<double>(*vec);
      bary->push_back(newCoords);
    } else {
639
      for (int i = rest - 1; i >= 1; i--) {
640
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
641
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
642
643
644
645
646
647
648
      }
    }
  }

  void Lagrange::setBary()
  {
    bary  =  &baryDimDegree[dim][degree];
649
650
651

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
652
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
653
654
655
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
	  int *coordInd = GET_MEMORY(int, i + 1); // indices of relevant coords
	  for (int k = 0; k < i+1; k++) { 
656
657
658
	    coordInd[k] = Global::getReferenceElement(dim)->
	      getVertexOfPosition(INDEX_OF_DIM(i,dim), j, k);
	  }
659
660
661
	  createCoords(coordInd, i + 1, 0, degree);
	  FREE_MEMORY(coordInd, int, i + 1);
	  if(static_cast<int>(bary->size()) == nBasFcts) return;
662
663
664
665
666
667
668
669
	}
      }
    }
  }

  int Lagrange::getNumberOfDOFs(int dim, int degree)
  {
    int result = 0;
670
671
    for (int i = 0; i <= degree; i++) {
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
    }
    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
701
    if (dimOfPosition == 0) {
702
703
704
705
      return &sortedVertex;
    }

    // edge
706
    if ((dimOfPosition == 1) && (degree == 2)) {
707
708
709
710
711
712
713
      return &sortedEdgeDeg2;
    }
      
    int vertex[3];
    int** dof = const_cast<int**>( el->getDOF());
    int verticesOfPosition = dimOfPosition + 1;

714
    for (int i = 0; i < verticesOfPosition; i++) {
715
716
717
718
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);
    }

719
720
721
    if (dimOfPosition == 1) {
      if (degree == 3) {
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
722
723
724
725
726
	  return sortedEdgeDeg3[0];
	} else {
	  return sortedEdgeDeg3[1];
	}
      } else { // degree == 4
727
	if (dof[vertex[0]][0] < dof[vertex[1]][0]) {
728
729
730
731
732
733
734
735
	  return sortedEdgeDeg4[0];
	} else {
	  return sortedEdgeDeg4[1];
	}
      }
    }

    // face
736
737
    if (dimOfPosition == 2) {
      if (degree == 3) {
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
	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
760
761
    if ((dimOfPosition == 3) && (degree == 4)) {
      return &sortedCenterDeg4;
762
763
764
765
766
767
768
769
770
771
    } 

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

  const BoundaryType* Lagrange::getBound(const ElInfo* el_info, 
					 BoundaryType* bound) const
  {
    static BoundaryType *bound_vec = NULL;
772
773
    static int bound_vec_size = 0;
    BoundaryType *result;
774

775
776
777
778
779
780
781
782
783
784
    if (bound) {
      result = bound;
    } else {
      // realloc memory if neccessary
      if (bound_vec_size < nBasFcts) {
	if (bound_vec) 
	  FREE_MEMORY(bound_vec, BoundaryType, bound_vec_size);
	bound_vec = GET_MEMORY(BoundaryType, nBasFcts);
	bound_vec_size = nBasFcts;
      }
785

786
      result = bound_vec;
787
788
789
790
791
    }

    el_info->testFlag(Mesh::FILL_BOUND);

    // boundaries
792
    int index = 0;
793
794
    int offset = 0;
    BoundaryType boundaryType;
795
    for (int i = dim - 1; i > 0; i--)
796
      offset += Global::getGeo(INDEX_OF_DIM(i, dim), dim);
797

798
    for (int i = 0; i < dim; i++) {
799
800
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
801
	boundaryType = el_info->getBoundary(j);
802
803
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
	for (int k = 0; k < kto; k++) {
804
805
806
	  result[index++] = boundaryType;
	}
      }
807
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
808
809
810
    }

    // interior nodes in the center
811
    for (int i = 0; i < (*nDOF)[CENTER]; i++) {
812
813
814
      result[index++] = INTERIOR;
    }

815
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
816
817
818
819
820
821
822
823
824
825

    return result;
  }


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

    static double *inter = NULL;
829
830
831
832
    static int inter_size = 0;

    double *rvec =  NULL;

833
    if (vec) {
834
835
836
837
838
839
840
841
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

842
    WorldVector<double> x;
843
844

    el_info->testFlag(Mesh::FILL_COORDS);
845
846
847
848
849
850
 
    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));
851
      }
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
	
      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 {
	  el_info->coordToWorld(*(*bary)[b_no[i]], &x);
	  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++) {
	el_info->coordToWorld(*(*bary)[i], &x);
	rvec[i] = (*f)(x);
      }
    }  
870
871
872
873
874
875
876
877
878
879
    
    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)
  {
880
881
    FUNCNAME("*Lagrange::interpol_d()");

882
    static WorldVector<double> *inter;
883
    WorldVector<double> *rvec = 
884
      vec ? vec : (inter?inter:inter= NEW WorldVector<double>[getNumber()]);
885
    WorldVector<double> x;
886
887
888
889
890

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

891
892
893
894
895
896
    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));
897
      }
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
      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 {
	  el_info->coordToWorld(*(*bary)[b_no[i]], &x);
	  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++) {
	el_info->coordToWorld(*(*bary)[i], &x);
	rvec[i] = (*f)(x);
      }
    }  
914
915
916
917
918
919
920
921
922
923
924
925
926
927
    
    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();
    const int *indi;
928
    int nrDOFs, n0, node0, num = 0;
929
930
931
    GeoIndex posIndex;

    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
932
    
933
    if (indices) {
934
935
      result = indices;
    } else {
936
937
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
938
939
940
941
      localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
942
    
943
    for (int pos = 0, j = 0; pos <= dim; pos++) {
944
945
946
      posIndex = INDEX_OF_DIM(pos, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

Thomas Witkowski's avatar
Thomas Witkowski committed
947
948
949
950
951
952
      if (nrDOFs) {
	n0 = admin->getNumberOfPreDOFs(posIndex);
	node0 = admin->getMesh()->getNode(posIndex);
	num = Global::getGeo(posIndex, dim);

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

955
956
	  for (int k = 0; k < nrDOFs; k++)
	    result[j++] = dof[node0][n0 + indi[k]];
957
958
959
	}
      }
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
960
    
961
962
963
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
964
965
966
967
968
969
970
971
972
973
974
  void Lagrange::getLocalIndicesVec(const Element* el,
				    const DOFAdmin *admin,
				    Vector<DegreeOfFreedom> *indices) const
  {
    if (indices->getSize() < nBasFcts) {
      indices->resize(nBasFcts);
    }

    getLocalIndices(el, admin, &((*indices)[0]));
  }

975
976
977
978
979
980
981
982
983
984
985
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

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

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
993
994
995
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

996
997
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
998
				fh->getName().c_str());
999
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
1000
1001

    Mesh* mesh = fh->getFESpace()->getMesh();
1002
    int n_phi = nBasFcts;
1003

1004
1005
1006
    if (!q) {
      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
    }
1007

1008
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
1009

1010
1011
1012
1013
1014
    if ((parametric = mesh->getParametric())) {
      ERROR_EXIT("not yet implemented\n");
    } else {
      wdetf_qp = GET_MEMORY(double, q->getNumPoints());
    }
1015

1016
    ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1017
1018
1019

    int numPoints = q->getNumPoints();
     
1020
1021
1022
    while (el_info) {
      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
      det = el_info->getDet();
1023

1024
1025
1026
1027
1028
1029
1030
1031
1032
      for (iq = 0; iq < numPoints; iq++) {
	el_info->coordToWorld(q->getLambda(iq), &x);
	wdetf_qp[iq] = q->getWeight(iq)*det*((*f)(x));
      }
      
      for (j = 0; j < n_phi; j++) {
	for (val = iq = 0; iq < numPoints; iq++)
	  val += quad_fast->getPhi(iq,j)*wdetf_qp[iq];
	(*fh)[dof[j]] += val;
1033
1034
      }

1035
1036
1037
      el_info = stack.traverseNext(el_info);
    }

1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
    FREE_MEMORY(wdetf_qp, double, q->getNumPoints());

    return;
  }


  // ===== 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
1056
1057
    if (n < 1) 
      return;

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

    (*drv)[dof_new] = (*drv)[dof0];
    dof_new = el->getChild(1)->getDOF(node, n0);  /*      newest vertex is center */
    (*drv)[dof_new] = (*drv)[dof0];
  }

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

1068
1069
    if (n < 1) 
      return;
1070

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

  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, BasisFunction* basFct)
  {
1082
    FUNCNAME("Lagrange::refineInter2_1d()");
1083

1084
1085
1086
1087
1088
1089
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1090
  
1091
1092
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1093
1094
1095
1096
1097

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

1098
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+1, n0);  /*      newest vertex is DIM */
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
    (*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] = 
1110
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1111
1112
1113
1114
1115
1116
1117

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

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

1128
1129
    if (n < 1) 
      return;
1130
1131

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1132
1133
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1134
  
1135
1136
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1137
1138
1139
1140
1141

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

1142
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
    (*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] = 
1154
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1155
1156
1157
1158
1159
1160
1161

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

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1162
1163
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1164
1165
1166
1167
1168
1169
1170

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]]  + 0.75 * (*drv)[pdof[5]];

    if (n > 1) {
      /****************************************************************************/
      /*  adjust the value at the midpoint of the common edge of neigh's children */
      /****************************************************************************/
      el = list->getElement(1);
      pdof = basFct->getLocalIndices(el, admin, NULL);
      
      cdof = el->getChild(0)->getDOF(node+1, n0); 
      (*drv)[cdof] = 
	-0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
	+ 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
    }
1185
1186
1187
1188
1189
1190
  }

  void Lagrange::refineInter2_3d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1191
    FUNCNAME("Lagrange::refineInter2_3d()");
1192

1193
1194
    if (n < 1) 
      return;
1195

1196
1197
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1198

1199
    DegreeOfFreedom pdof[10];
1200
1201
    basFct->getLocalIndices(el, admin, pdof);

1202
1203
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1204
1205
1206
1207
1208

    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/

1209
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1210
1211
1212

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1213
1214
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1215
    (*drv)[cdof[8]] = 
1216
1217
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1218
    (*drv)[cdof[9]] = 
1219
1220
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1221
1222
1223
1224
1225

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1226
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1227
    (*drv)[cdofi] = 
1228
1229
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1230
1231
1232
1233
1234

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1235
1236
1237
1238
1239
1240
1241
1242
1243
1244
1245
1246
1247
1248
1249
1250
1251
1252
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
    for (int i = 1; i < n; i++) {
      el = list->getElement(i);
      basFct->getLocalIndices(el, admin, pdof);
      
      int lr_set = 0;
      if (list->getNeighbourElement(i, 0) && list->getNeighbourNr(i, 0) < i)
	lr_set = 1;
      
      if (list->getNeighbourElement(i, 1) && list->getNeighbourNr(i, 1) < i)
	lr_set += 2;
      
      TEST_EXIT(lr_set)("no values set on both neighbours\n");
      
      /****************************************************************************/
      /*  values on child[0]                                                      */
      /****************************************************************************/
      
      switch (lr_set) {
      case 1:
	cdofi = el->getChild(0)->getDOF(node0+4, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
	break;
      case 2:
	cdofi = el->getChild(0)->getDOF(node0+5, n0);      
	(*drv)[cdofi] = 
	  (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
	   + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
	
1265
      }
1266
    }
1267
1268
1269
1270
1271
1272
1273
    return;
  }

  void Lagrange::refineInter3_1d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
1274
    FUNCNAME("Lagrange::refineInter3_1d()");
1275

1276
1277
    if (n < 1) 
      return;
1278

1279