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

namespace AMDiS {

19
20
  using boost::lexical_cast;

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


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

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();
  }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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


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

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

    int dimOfPosition = DIM_OF_INDEX(position, dim);

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

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

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

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

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

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

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

	return sortedFaceDeg4[no];
      }
    }

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

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

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

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

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

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

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


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

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

812
    double *rvec = NULL;
813

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

      localVecSize = nBasFcts;
      rvec = localVec;
826
827
    } 

828
    WorldVector<double> x;
829
830

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

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

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

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

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


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

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

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

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

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

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

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

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

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

953
954
955
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
956
957
958
959
  void Lagrange::getLocalIndicesVec(const Element* el,
				    const DOFAdmin *admin,
				    Vector<DegreeOfFreedom> *indices) const
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
960
    if (indices->getSize() < nBasFcts)
Thomas Witkowski's avatar
Thomas Witkowski committed
961
962
963
964
965
      indices->resize(nBasFcts);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
  void Lagrange::getLocalDofPtrVec(const Element *el, 
				   const DOFAdmin *admin,
				   std::vector<const DegreeOfFreedom*>& vec) const
  {
    if (static_cast<int>(vec.size()) < nBasFcts)
      vec.resize(nBasFcts);

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

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

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

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

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

996
997
998
999
1000
1001
1002
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
    ERROR_EXIT("not yet\n");
  }

1003
  void Lagrange::l2ScpFctBas(Quadrature *quad,
1004
1005
1006
			     AbstractFunction<double, WorldVector<double> >* f,
			     DOFVector<double>* fh)
  {
1007
1008
    FUNCNAME("Lagrange::l2ScpFctBas()");

1009
    TEST_EXIT_DBG(fh)("no DOF_REAL_VEC fh\n");
Thomas Witkowski's avatar
Thomas Witkowski committed
1010
1011
1012
1013
    TEST_EXIT_DBG(fh->getFESpace())
      ("no fe_space in DOF_REAL_VEC %s\n", fh->getName().c_str());
    TEST_EXIT_DBG(fh->getFESpace()->getBasisFcts() == this)
      ("wrong basis fcts for fh\n");
1014
1015
    TEST_EXIT_DBG(!fh->getFESpace()->getMesh()->getParametric())
      ("Not yet implemented!");
1016

1017
1018
    if (!quad)
      quad = Quadrature::provideQuadrature(dim, 2 * degree - 2);
1019

Thomas Witkowski's avatar
Thomas Witkowski committed
1020
    const FastQuadrature *quad_fast = 
1021
1022
1023
1024
1025
      FastQuadrature::provideFastQuadrature(this, *quad, INIT_PHI);
    double *wdetf_qp = new double[quad->getNumPoints()];
    int nPoints = quad->getNumPoints();
    DOFAdmin *admin = fh->getFESpace()->getAdmin();
    WorldVector<double> x;
1026

1027
1028
    TraverseStack stack;
    ElInfo *el_info = stack.traverseFirst(fh->getFESpace()->getMesh(), -1, 
Thomas Witkowski's avatar
Thomas Witkowski committed
1029
					  Mesh::CALL_LEAF_EL | Mesh::FILL_COORDS);
1030
    while (el_info) {
1031
1032
      const DegreeOfFreedom  *dof = getLocalIndices(el_info->getElement(), admin, NULL);
      double det = el_info->getDet();
1033

1034
1035
1036
      for (int iq = 0; iq < nPoints; iq++) {
	el_info->coordToWorld(quad->getLambda(iq), x);
	wdetf_qp[iq] = quad->getWeight(iq) * det * ((*f)(x));
1037
1038
      }
      
1039
1040
1041
1042
      for (int j = 0; j < nBasFcts; j++) {
	double val = 0.0;
	for (int iq = 0; iq < nPoints; iq++)
	  val += quad_fast->getPhi(iq, j) * wdetf_qp[iq];
Thomas Witkowski's avatar
Thomas Witkowski committed
1043

1044
	(*fh)[dof[j]] += val;
1045
1046
      }

1047
1048
1049
      el_info = stack.traverseNext(el_info);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
1050
    delete [] wdetf_qp;
1051
1052
1053
1054
1055
1056
1057
  }


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

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

1060
1061
1062
1063
1064
1065
    if (n < 1) 
      return;

    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(CENTER);
    Element *el = list->getElement(0);
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);
1066
1067
1068
1069
    // Parent center
    DegreeOfFreedom dof0 = el->getDOF(node, n0);           
    // Newest vertex is center 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(node, n0);  
1070
1071

    (*drv)[dof_new] = (*drv)[dof0];
1072
1073
    // Newest vertex is center 
    dof_new = el->getChild(1)->getDOF(node, n0);  
1074
1075
1076
    (*drv)[dof_new] = (*drv)[dof0];
  }

1077
1078
  void Lagrange::refineInter1(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
			      BasisFunction* basFct)
1079
  {
1080
    FUNCNAME("Lagrange::refineInter1_1d()");
1081

1082
1083
    if (n < 1) 
      return;
1084

1085
1086
1087
    int dim = drv->getFESpace()->getMesh()->getDim();
    int n0 = drv->getFESpace()->getAdmin()->getNumberOfPreDOFs(VERTEX);
    Element *el = list->getElement(0);
1088
1089
1090
1091
1092
1093
    // 1st endpoint of refinement edge
    DegreeOfFreedom dof0 = el->getDOF(0, n0);
    // 2nd endpoint of refinement edge          
    DegreeOfFreedom dof1 = el->getDOF(1, n0);
    // newest vertex is DIM 
    DegreeOfFreedom dof_new = el->getChild(0)->getDOF(dim, n0); 
1094
    (*drv)[dof_new] = 0.5 * ((*drv)[dof0] + (*drv)[dof1]);
1095
1096
  }

1097
1098
  void Lagrange::refineInter2_1d(DOFIndexed<double> *drv, RCNeighbourList* list, int n, 
				 BasisFunction* basFct)
1099
  {
1100
    FUNCNAME("Lagrange::refineInter2_1d()");
1101

1102
1103
1104
1105
1106
1107
    if (n < 1) 
      return;
    
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1108
  
1109
1110
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1111
1112
1113
1114
1115

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

1116
1117
    // newest vertex is DIM
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node + 1, n0);
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
    (*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] = 
1129
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1130
1131
1132
1133
1134
1135
1136

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node, n0); 
    (*drv)[cdof] = 
1137
      -0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[2]];
1138
1139
1140
1141
1142
1143
  }

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

1146
1147
    if (n < 1) 
      return;
1148
1149

    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1150
1151
    Element *el = list->getElement(0);
    const DegreeOfFreedom *pdof = basFct->getLocalIndices(el, admin, NULL);
1152
  
1153
1154
    int node = drv->getFESpace()->getMesh()->getNode(VERTEX);        
    int n0 = admin->getNumberOfPreDOFs(VERTEX);
1155
1156
1157
1158
1159

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

1160
    DegreeOfFreedom cdof = el->getChild(0)->getDOF(node+2, n0);  /*      newest vertex is DIM */
1161
1162
1163
1164
1165
1166
1167
1168
1169
1170
1171
    (*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] = 
1172
      0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] + 0.75 * (*drv)[pdof[5]];
1173
1174
1175
1176
1177
1178
1179

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

    cdof = el->getChild(0)->getDOF(node+1, n0); 
    (*drv)[cdof] =
1180
1181
      -0.125 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[5]]
      + 0.5 * ((*drv)[pdof[3]] + (*drv)[pdof[4]]);
1182
1183
1184
1185
1186
1187
1188

    /****************************************************************************/
    /*  midpoint of edge on child[1] at the refinement edge                     */
    /****************************************************************************/
  
    cdof = el->getChild(1)->getDOF(node+1, n0); 
    (*drv)[cdof] = 
1189
1190
1191
1192
1193
1194
1195
1196
1197
1198
1199
1200
1201
1202
      -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]]);
    }
1203
1204
1205
1206
1207
1208
  }

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

1211
1212
    if (n < 1) 
      return;
1213

1214
1215
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1216

1217
    DegreeOfFreedom pdof[10];
1218
1219
    basFct->getLocalIndices(el, admin, pdof);

1220
1221
    int node0 = drv->getFESpace()->getMesh()->getNode(EDGE);
    int n0 = admin->getNumberOfPreDOFs(EDGE);
1222
1223
1224
1225
1226

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

1227
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1228
1229
1230

    (*drv)[cdof[3]] = ((*drv)[pdof[4]]);
    (*drv)[cdof[6]] = 
1231
1232
      (0.375 * (*drv)[pdof[0]] - 0.125 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1233
    (*drv)[cdof[8]] = 
1234
1235
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[5]] + (*drv)[pdof[7]]));
1236
    (*drv)[cdof[9]] = 
1237
1238
      (0.125 * (-(*drv)[pdof[0]] - (*drv)[pdof[1]]) + 0.25 * (*drv)[pdof[4]]
       + 0.5 * ((*drv)[pdof[6]] + (*drv)[pdof[8]]));
1239
1240
1241
1242
1243

    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/
  
1244
    DegreeOfFreedom cdofi = el->getChild(1)->getDOF(node0 + 2, n0);
1245
    (*drv)[cdofi] = 
1246
1247
      (-0.125 * (*drv)[pdof[0]] + 0.375 * (*drv)[pdof[1]] 
       + 0.75 * (*drv)[pdof[4]]);
1248
1249
1250
1251
1252

    /****************************************************************************/
    /*   adjust neighbour values                                                */
    /****************************************************************************/
  
1253
1254
1255
1256
1257
1258
1259
1260
1261
1262
1263
1264
1265
1266
1267
1268
1269
1270
1271
1272
1273
1274
1275
1276
1277
1278
1279
1280
1281
1282
    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]]));
	
1283
      }
1284
    }
1285
1286
1287
1288
1289
1290
  }

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

1293
1294
    if (n < 1) 
      return;
1295

1296
1297
    Element *el = list->getElement(0);
    const DOFAdmin *admin = drv->getFESpace()->getAdmin();
1298

Thomas Witkowski's avatar
Thomas Witkowski committed
1299
    DegreeOfFreedom *pdof = new DegreeOfFreedom[basFct->getNumber()];
1300
1301
1302
1303
1304
    basFct->getLocalIndices(el, admin, pdof);
  
    /****************************************************************************/
    /*  values on child[0]                                                      */
    /****************************************************************************/
1305
    const DegreeOfFreedom *cdof = basFct->getLocalIndices(el->getChild(0), admin, NULL);
1306
1307

    (*drv)[cdof[2]] = 
1308
1309
      (-0.0625 * ((*drv)[pdof[0]] + (*drv)[pdof[1]]) 
       + 0.5625 * ((*drv)[pdof[7]] + (*drv)[pdof[8]]));
1310
    (*drv)[cdof[3]] = 
1311
1312
1313
1314
      (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]];
1315
    (*drv)[cdof[6]] =  
1316
1317
1318
1319
      (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]]));
1320
    (*drv)[cdof[9]] = 
1321
1322
1323
      (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]]);
1324
1325
1326
1327
1328
1329
1330
1331
1332
  
    /****************************************************************************/
    /*  values on child[1]                                                      */
    /****************************************************************************/

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

    (*drv)[cdof[5]] =  (*drv)[pdof[8]];
    (*drv)[cdof[6]] =  
1333
1334
      (0.0625 * (*drv)[pdof[0]] + 0.9375 * (*drv)[pdof[8]]
       + 0.3125 * ((*drv)[pdof[1]] - (*drv)[pdof[7]]));
1335
    (*drv)[cdof[9]] =  
1336
1337
1338
      (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]]);
1339
1340

    if (n <= 1) {
Thomas Witkowski's avatar
Thomas Witkowski committed
1341
      delete [] pdof;
1342
1343
1344
1345
1346
1347
1348
1349
1350
1351
1352
1353
1354
1355
1356
1357
      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]] = 
1358
1359
1360
1361
      (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]]));
1362
    (*drv)[cdof[9]] = 
1363
1364
1365
      (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]]);
1366
1367
1368
1369
    /****************************************************************************/
    /*  (*drv)alues on neigh's child[0]                                         */
    /****************************************************************************/

1370
1371
1372
    int node = drv->getFESpace()->getMesh()->getNode(CENTER);  
    int n0 = admin->getNumberOfPreDOFs(CENTER);
    int dof9 = el->getChild(1)->getDOF(node, n0);
1373
1374

    (*drv)[dof9] =  
1375
1376
1377
      (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]]);
1378

Thomas Witkowski's avatar
Thomas Witkowski committed
1379
    delete [] pdof;
1380
1381
1382
1383
1384
1385
1386
  }

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