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);

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

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

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

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

    return result;
  }


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

    static double *inter = NULL;
831
832
833
834
    static int inter_size = 0;

    double *rvec =  NULL;

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

844
    WorldVector<double> x;
845
846

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

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

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

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

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

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

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

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

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

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

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

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

    Mesh* mesh = fh->getFESpace()->getMesh();
1004
    int n_phi = nBasFcts;
1005

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

1010
    const FastQuadrature *quad_fast = FastQuadrature::provideFastQuadrature(this, *q, INIT_PHI);
1011

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

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

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

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

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

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

    return;
  }


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

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

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

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

1070
1071
    if (n < 1) 
      return;
1072

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

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

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

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

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

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

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

1130
1131
    if (n < 1) 
      return;
1132
1133

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

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

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

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

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

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

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

1195
1196
    if (n < 1) 
      return;
1197

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

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

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

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

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

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

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

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

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

1278
1279
    if (n < 1) 
      return;
1280

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

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

    (*drv)[cdof[2]] = 
1293
1294
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1295
    (*drv)[cdof[3]] = 
1296
1297
1298
1299
      (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]];
1300
    (*drv)[cdof[6]] =  
1301
1302
1303
1304
      (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]]));
1305
    (*drv)[cdof[9]] = 
1306
1307
1308
      (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]]);
1309
1310
1311
1312
1313
1314
1315
1316
1317
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1318
1319
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1320
    (*drv)[cdof[9]] =  
1321
1322
1323
      (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]]);
1324
1325
1326
1327
1328
1329
1330
1331
1332
1333
1334
1335
1336
1337
1338
1339
1340
1341
1342

    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]] = 
1343
1344
1345
1346
      (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]]));
1347
    (*drv)[cdof[9]] = 
1348
1349
1350
      (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]]);
1351
1352
1353
1354
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

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

    (*drv)[dof9] =  
1360
1361
1362
      (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]]);
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372

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

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

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

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

    (*drv)[cdof[2]] =
1390
1391
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1392
    (*drv)[cdof[3]] = 
1393
1394
1395
1396
      (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]];
1397
    (*drv)[cdof[6]] =  
1398
1399
1400
1401
      (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]]));
1402
    (*drv)[cdof[9]] = 
1403
1404
1405
      (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]]);
1406
1407
1408
1409
1410
1411
1412

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

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

1413
    (*drv)[cdof[5]] = (*drv)[pdof[8]];
1414
    (*drv)[cdof[6]] =  
1415
1416
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1417
    (*drv)[cdof[9]] =  
1418
1419
1420
      (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]]);
1421

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

    /****************************************************************************/
    /*  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]] =  
1439
1440
1441
1442
      (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]]));
1443
    (*drv)[cdof[9]] = 
1444
1445
1446
      (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]]);
1447
1448
1449
1450
    /****************************************************************************/
    /*  values on neigh's child[0]                                              */
    /****************************************************************************/

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

    (*drv)[dof9] =  
1456
1457
1458
      (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]]);
1459
1460
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

    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;

1584
	TEST_EXIT_DBG(lr_set)("no values set on both neighbours\n");
1585
1586
1587
1588
1589
1590
1591
1592
1593
1594
1595
1596
1597
1598
1599
1600
1601