Lagrange.cc 179 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
#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 {

  ::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];
  ::std::list<Lagrange*> Lagrange::allBasFcts;


  Lagrange::Lagrange(int dim, int degree)
    : BasisFunction(::std::string("Lagrange"), dim, degree)
  {
    // 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)
  {
    ::std::list<Lagrange*>::iterator it;

    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) {
      if (degree == 4) {
762
763
764
765
766
767
768
769
770
771
772
773
	return &sortedCenterDeg4;
      }
    } 

    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;
774
775
    static int bound_vec_size = 0;
    BoundaryType *result;
776

777
778
779
780
781
782
783
784
785
786
    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;
      }
787

788
      result = bound_vec;
789
790
791
792
793
    }

    el_info->testFlag(Mesh::FILL_BOUND);

    DimVec<int> parts(dim, NO_INIT);
794
    for (int i = 0; i < dim + 1; i++)
795
      parts[i] = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
796
797
798
799
800
801

    int index = 0;

    // boundaries
    int offset = 0;
    BoundaryType boundaryType;
802
    for (int i = dim - 1; i > 0; i--)
803
804
      offset += parts[i];

805
806
    for (int i = 0; i < dim; i++) {
      for (int j = offset; j < offset + parts[i]; j++) {
807
	boundaryType = el_info->getBoundary(j);
808
	for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
809
810
811
812
813
814
815
	  result[index++] = boundaryType;
	}
      }
      offset -= parts[i+1];
    }

    // interior nodes in the center
816
    for (int i = 0; i < (*nDOF)[CENTER]; i++) {
817
818
819
      result[index++] = INTERIOR;
    }

820
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
821
822
823
824
825
826
827
828
829
830

    return result;
  }


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

    static double *inter = NULL;
834
835
836
837
    static int inter_size = 0;

    double *rvec =  NULL;

838
    if (vec) {
839
840
841
842
843
844
845
846
      rvec = vec;
    } else {
      if(inter) FREE_MEMORY(inter, double, nBasFcts);
      inter = GET_MEMORY(double, nBasFcts);
      inter_size = nBasFcts;
      rvec = inter;
    } 

847
    WorldVector<double> x;
848
849

    el_info->testFlag(Mesh::FILL_COORDS);
850
851
852
853
854
855
 
    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));
856
      }
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
	
      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);
      }
    }  
875
876
877
878
879
880
881
882
883
884
    
    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)
  {
885
886
    FUNCNAME("*Lagrange::interpol_d()");

887
    static WorldVector<double> *inter;
888
    WorldVector<double> *rvec = 
889
      vec ? vec : (inter?inter:inter= NEW WorldVector<double>[getNumber()]);
890
    WorldVector<double> x;
891
892
893
894
895

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

896
897
898
899
900
901
    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));
902
      }
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
      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);
      }
    }  
919
920
921
922
923
924
925
926
927
928
929
930
931
932
    
    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;
933
    int nrDOFs, n0, node0, num = 0;
934
935
936
    GeoIndex posIndex;

    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
937
    
938
    if (indices) {
939
940
      result = indices;
    } else {
941
942
      if (localVec) 
	FREE_MEMORY(localVec, DegreeOfFreedom, localVecSize);
943
944
945
946
      localVec = GET_MEMORY(DegreeOfFreedom, nBasFcts);
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
947
    
948
    for (int pos = 0, j = 0; pos <= dim; pos++) {
949
950
951
952
953
954
      posIndex = INDEX_OF_DIM(pos, dim);
      n0 = admin->getNumberOfPreDOFs(posIndex);
      node0 = admin->getMesh()->getNode(posIndex);
      num = Global::getGeo(posIndex, dim);
      nrDOFs = admin->getNumberOfDOFs(posIndex);

955
956
      for (int i = 0; i < num; node0++, i++) {
	if (nrDOFs) {
957
958
	  indi = orderOfPositionIndices(el, posIndex, i);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
968
969
970
971
972
973
974
975
976
977
978
  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]));
  }

979
980
981
982
983
984
985
986
987
988
989
  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)
  {
990
991
992
993
994
995
996
    FUNCNAME("Lagrange::l2ScpFctBas()");

    TraverseStack stack;
    double *wdetf_qp = NULL;
    double val, det;
    WorldVector<double> x;
    int iq, j;
997
998
999
    const DegreeOfFreedom  *dof;
    const Parametric *parametric;

1000
1001
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
    TEST_EXIT_DBG(fh->getFESpace())("no fe_space in DOF_REAL_VEC %s\n",
1002
				fh->getName().c_str());
1003
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)("wrong basis fcts for fh\n");
1004
1005

    Mesh* mesh = fh->getFESpace()->getMesh();
1006
    int n_phi = nBasFcts;
1007

1008
1009
1010
    if (!q) {
      q = Quadrature::provideQuadrature(dim, 2 * degree - 2);
    }
1011

1012
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
1013

1014
1015
1016
1017
1018
    if ((parametric = mesh->getParametric())) {
      ERROR_EXIT("not yet implemented\n");
    } else {
      wdetf_qp = GET_MEMORY(double, q->getNumPoints());
    }
1019

1020
    ElInfo *el_info = stack.traverseFirst(mesh, -1, Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1021
1022
1023

    int numPoints = q->getNumPoints();
     
1024
1025
1026
    while (el_info) {
      dof = getLocalIndices(el_info->getElement(), (fh->getFESpace()->getAdmin()), NULL);      
      det = el_info->getDet();
1027

1028
1029
1030
1031
1032
1033
1034
1035
1036
      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;
1037
1038
      }

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

1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
    FREE_MEMORY(wdetf_qp, double, q->getNumPoints());

    return;
  }


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

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

1054
1055
1056
1057
1058
1059
1060
1061
    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 */
1062
1063
1064
1065
1066
1067
1068
1069

    (*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)
  {
1070
    FUNCNAME("Lagrange::refineInter1_1d()");
1071

1072
1073
    if (n < 1) 
      return;
1074

1075
1076
1077
1078
1079
1080
1081
    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]);
1082
1083
1084
1085
  }

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

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

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

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

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

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

1132
1133
    if (n < 1) 
      return;
1134
1135

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

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

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

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

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

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186
1187
1188
      -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]]);
    }
1189
1190
1191
1192
1193
1194
  }

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

1197
1198
    if (n < 1) 
      return;
1199

1200
1201
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1202

1203
    DegreeOfFreedom pdof[10];
1204
1205
    basFct->getLocalIndices(el, admin, pdof);

1206
1207
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1208
1209
1210
1211
1212

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

1213
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1214
1215
1216

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

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

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
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
1265
1266
1267
1268
    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]]));
	
1269
      }
1270
    }
1271
1272
1273
1274
1275
1276
1277
    return;
  }

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

1280
1281
    if (n < 1) 
      return;
1282

1283
1284
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1285

1286
    DegreeOfFreedom *pdof = GET_MEMORY(DegreeOfFreedom, basFct->getNumber());
1287
1288
1289
1290
1291
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1292
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1293
1294

    (*drv)[cdof[2]] = 
1295
1296
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1297
    (*drv)[cdof[3]] = 
1298
1299
1300
1301
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
1302
    (*drv)[cdof[6]] =  
1303
1304
1305
1306
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1307
    (*drv)[cdof[9]] = 
1308
1309
1310
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1311
1312
1313
1314
1315
1316
1317
1318
1319
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1320
1321
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1322
    (*drv)[cdof[9]] =  
1323
1324
1325
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342
1343
1344

    if (n <= 1) {
      FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber());
      return;
    }
    /****************************************************************************/
    /*  adjust the value on the neihgbour                                       */
    /****************************************************************************/

    el = list->getElement(1);
    basFct->getLocalIndices(el, admin, pdof);

    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/
    cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
  
    (*drv)[cdof[5]] =  (*drv)[pdof[9]];
    (*drv)[cdof[6]] = 
1345
1346
1347
1348
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1349
    (*drv)[cdof[9]] = 
1350
1351
1352
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.12500 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1353
1354
1355
1356
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

1357
1358
1359
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1360
1361

    (*drv)[dof9] =  
1362
1363
1364
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) +  0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75  *(*drv)[pdof[9]]);
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374

    FREE_MEMORY(pdof, DegreeOfFreedom, basFct->getNumber());
    return;
  }

  void Lagrange::refineInter3_2d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, 
				 BasisFunction* basFct)
  {
1375
    FUNCNAME("Lagrange::refineInter3_2d()");
1376

1377
1378
1379
1380
1381
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1382

1383
    DegreeOfFreedom pdof[10];
1384
1385
1386
1387
1388
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1389
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1390
1391

    (*drv)[cdof[2]] =
1392
1393
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1394
    (*drv)[cdof[3]] = 
1395
1396
1397
1398
      (0.3125 * ((*drv)[pdof[0]] - (*drv)[pdof[8]]) + 0.0625 * (*drv)[pdof[1]]
       + 0.9375 * (*drv)[pdof[7]]);
    (*drv)[cdof[4]] = (*drv)[pdof[7]];
    (*drv)[cdof[5]] = (*drv)[pdof[9]];
1399
    (*drv)[cdof[6]] =  
1400
1401
1402
1403
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]])
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1404
    (*drv)[cdof[9]] = 
1405
1406
1407
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.125 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1408
1409
1410
1411
1412
1413
1414

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

    cdof = basFct->getLocalIndices(el->getChild(1), admin, NULL);

1415
    (*drv)[cdof[5]] = (*drv)[pdof[8]];
1416
    (*drv)[cdof[6]] =  
1417
1418
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1419
    (*drv)[cdof[9]] =  
1420
1421
1422
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1423

1424
1425
    if (n <= 1)
      return;
1426
1427
1428
1429
1430
1431
1432
1433
1434
1435
1436
1437
1438
1439
1440

    /****************************************************************************/
    /*  adjust the value on the neihgbour                                       */
    /****************************************************************************/

    el = list->getElement(1);
    basFct->getLocalIndices(el, admin, pdof);

    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/
    cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
  
    (*drv)[cdof[5]] =  (*drv)[pdof[9]];
    (*drv)[cdof[6]] =  
1441
1442
1443
1444
      (0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       - 0.25 * ((*drv)[pdof[3]] + (*drv)[pdof[6]])
       + 0.5 * ((*drv)[pdof[4]] + (*drv)[pdof[5]] + (*drv)[pdof[9]])
       - 0.0625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1445
    (*drv)[cdof[9]] = 
1446
1447
1448
      (0.0625 * (-(*drv)[pdof[0]] + (*drv)[pdof[1]]) - 0.12500 * (*drv)[pdof[3]]
       + 0.375 * (*drv)[pdof[6]] + 0.1875 * ((*drv)[pdof[7]] - (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1449
1450
1451
1452
    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/

1453
1454
1455
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1456
1457

    (*drv)[dof9] =  
1458
1459
1460
      (0.0625 * ((*drv)[pdof[0]] - (*drv)[pdof[1]]) +  0.375 * (*drv)[pdof[3]]
       - 0.125 * (*drv)[pdof[6]] + 0.1875 * (-(*drv)[pdof[7]] + (*drv)[pdof[8]])
       + 0.75 * (*drv)[pdof[9]]);
1461
1462
1463
1464
1465
1466
1467
1468
1469
1470
1471
1472
1473
1474
1475
1476
1477
1478
1479
1480
1481
1482
1483
1484
1485
1486
1487
1488
1489
1490
1491
1492
1493
1494
1495
1496
1497
1498
1499
1500
1501
1502
1503
1504
1505
1506
1507
1508
1509
1510
1511
1512
1513
1514
1515
1516
1517
1518
1519
1520
1521
1522
1523
1524
1525
1526
1527
1528
1529
1530
1531
1532
1533
1534
1535
1536
1537
1538
1539
1540
1541
1542
1543
1544
1545
1546
1547
1548
1549
1550
1551
1552
1553
1554
1555
1556
1557
1558
1559
1560
1561
1562
1563
1564
1565
1566
1567
1568
1569
1570
1571
1572
1573
1574
1575
1576
1577
1578
1579
1580
1581
1582
1583
1584
1585

    return;
  }

  void Lagrange::refineInter3_3d(DOFIndexed<double> *drv, 
				 RCNeighbourList* list, 
				 int n, BasisFunction* basFct)
  {
    FUNCNAME("Lagrange::refineInter3_3d");
    Element *el;
    const DegreeOfFreedom *cd;
    DegreeOfFreedom pd[20], cdi;
    int       i, typ, lr_set, node0, n0;
    const DOFAdmin *admin;

    if (n < 1) return;
    el = list->getElement(0);
    typ = list->getType(0);

    admin = drv->getFESpace()->getAdmin();

    basFct->getLocalIndices(el, admin, pd);

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

    cd = basFct->getLocalIndices(el->getChild(0), admin, NULL);

    (*drv)[cd[3]] = 
      (0.0625*(-(*drv)[pd[0]] - (*drv)[pd[1]]) 
       + 0.5625*((*drv)[pd[4]] + (*drv)[pd[5]]));
    (*drv)[cd[8]] = 
      (0.3125*((*drv)[pd[0]] - (*drv)[pd[5]]) + 0.0625*(*drv)[pd[1]]
       + 0.937500*(*drv)[pd[4]]);
    (*drv)[cd[9]] = (*drv)[pd[4]];
    (*drv)[cd[12]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.25*(-(*drv)[pd[6]] - (*drv)[pd[10]])
       + 0.5*((*drv)[pd[7]] + (*drv)[pd[11]] + (*drv)[pd[19]]));
    (*drv)[cd[13]] = (*drv)[pd[19]];
    (*drv)[cd[14]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.25*(-(*drv)[pd[8]] - (*drv)[pd[12]]) 
       + 0.5*((*drv)[pd[9]] + (*drv)[pd[13]] + (*drv)[pd[18]]));
    (*drv)[cd[15]] = ((*drv)[pd[18]]);
    (*drv)[cd[16]] = 
      (0.0625*((*drv)[pd[0]] + (*drv)[pd[1]] - (*drv)[pd[4]] - (*drv)[pd[5]])
       + 0.125*(-(*drv)[pd[6]]-(*drv)[pd[8]]-(*drv)[pd[10]]-(*drv)[pd[12]])
       + 0.5*((*drv)[pd[16]] + (*drv)[pd[17]]) 
       + 0.25*((*drv)[pd[18]] + (*drv)[pd[19]]));
    (*drv)[cd[17]] = 
      (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) 
       + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[8]]
       - 0.125*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
    (*drv)[cd[18]] = 
      (0.0625*(-(*drv)[pd[0]] + (*drv)[pd[1]]) 
       + 0.1875*((*drv)[pd[4]] - (*drv)[pd[5]]) + 0.375*(*drv)[pd[6]]
       - 0.125*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

    cd = basFct->getLocalIndices(el->getChild(1), admin, NULL);

    if (typ == 0)
      {
	/****************************************************************************/
	/*  parent of el_type 0                                                     */
	/****************************************************************************/

	(*drv)[cd[8]] = 
	  (0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
	   + 0.9375*(*drv)[pd[5]]);
	(*drv)[cd[9]] = (*drv)[pd[5]];
	(*drv)[cd[17]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]])
	   - 0.125*(*drv)[pd[6]] + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
	(*drv)[cd[18]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]] 
	   + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
      }
    else
      {
	/****************************************************************************/
	/*  parent of el_type 0                                                     */
	/****************************************************************************/

	(*drv)[cd[8]] = 
	  (0.0625*(*drv)[pd[0]] + 0.3125*((*drv)[pd[1]] - (*drv)[pd[4]])
	   + 0.9375*(*drv)[pd[5]]);
	(*drv)[cd[9]] = (*drv)[pd[5]];
	(*drv)[cd[17]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) 
	   + 0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[8]]
	   + 0.375*(*drv)[pd[12]] + 0.75*(*drv)[pd[18]]);
	(*drv)[cd[18]] = 
	  (0.0625*((*drv)[pd[0]] - (*drv)[pd[1]]) + 
	   0.1875*(-(*drv)[pd[4]] + (*drv)[pd[5]]) - 0.125*(*drv)[pd[6]]
	   + 0.375*(*drv)[pd[10]] + 0.75*(*drv)[pd[19]]);
      }

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/

    node0 = drv->getFESpace()->getMesh()->getNode(FACE);
    n0 = admin->getNumberOfPreDOFs(FACE);
  
    for (i = 1; i < n; i++)
      {
	el = list->getElement(i);
	typ = list->getType(i);
	basFct->getLocalIndices(el, admin, pd);

	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;

1586
	TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n");
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601
1602
1603
1604
1605
1606
1607
1608
1609
1610
1611
1612
1613
1614
1615
1616
1617
1618
1619
1620
1621
1622
1623
1624
1625
1626
1627
1628
1629