Lagrange.cc 169 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


13
14
15
16
#include <stdio.h>
#include <algorithm>
#include <list>
#include "boost/lexical_cast.hpp"
17
18
19
20
21
22
23
24
25
26
27
28
29
30
#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 {

31
32
  using boost::lexical_cast;

33
34
35
36
37
38
  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];
39
  std::list<Lagrange*> Lagrange::allBasFcts;
40
41
42


  Lagrange::Lagrange(int dim, int degree)
43
    : BasisFunction(std::string("Lagrange"), dim, degree)
44
45
  {
    // set name
46
    name += lexical_cast<std::string>(dim) + " " + lexical_cast<std::string>(degree);
47
48
49
50
51
52
53
54
55
56
57
58
59

    // set nDOF
    setNDOF();

    // set barycentric coordinates
    setBary();

    // set function pointer
    setFunctionPointer();
  }

  Lagrange::~Lagrange()
  {
60
    for (int i = 0; i < static_cast<int>(bary->size()); i++)
61
      if ((*bary)[i]) {
62
	delete (*bary)[i];
63
64
	(*bary)[i] = NULL;
      }
65
66
67
68
  }

  Lagrange* Lagrange::getLagrange(int dim, int degree)
  {
69
    std::list<Lagrange*>::iterator it;
Thomas Witkowski's avatar
Thomas Witkowski committed
70
71
    for (it = allBasFcts.begin(); it != allBasFcts.end(); it++)
      if (((*it)->dim == dim) && ((*it)->degree == degree))
72
73
	return (*it);

74
    Lagrange* newLagrange = new Lagrange(dim, degree);
75
76
77
78
    allBasFcts.push_back(newLagrange);
    return newLagrange;
  }

79
80
  void Lagrange::clear()
  {    
81
82
    for (std::list<Lagrange*>::iterator it = allBasFcts.begin(); 
	 it != allBasFcts.end(); it++)
83
      if (*it) {
84
	delete *it;
85
86
	*it = NULL;
      }     
87
88
  }

89
90
  void Lagrange::setFunctionPointer()
  {
91
    if (static_cast<int>(phiDimDegree[dim][degree].size()) == 0) {
92
      // for all positions
93
      for (int i = 0; i < dim + 1; i++) {
94
	// no vertex dofs for degree 0 ?
95
96
	if (degree == 0 && i != dim) 
	  continue;
97
	// for all vertices/edges/...
98
	for (int j = 0; j < Global::getGeo(INDEX_OF_DIM(i, dim),dim); j++) {
99
	  // for all dofs at this position
100
	  for (int k = 0; k < (*nDOF)[INDEX_OF_DIM(i, dim)]; k++) {
101
	    // basis function
102
	    phiDimDegree[dim][degree].push_back(new Phi(this, 
103
							INDEX_OF_DIM(i ,dim), j, k));
104
	    // gradients
105
	    grdPhiDimDegree[dim][degree].push_back(new GrdPhi(this, 
106
107
							      INDEX_OF_DIM(i, dim),
							      j, k));
108
	    // D2-Matrices
109
	    D2PhiDimDegree[dim][degree].push_back(new D2Phi(this, 
110
111
							    INDEX_OF_DIM(i, dim),
							    j, k));
112
113
114
115
	  }
	}
      }
    }
116
    phi = &phiDimDegree[dim][degree];
117
    grdPhi = &grdPhiDimDegree[dim][degree];
118
    d2Phi = &D2PhiDimDegree[dim][degree];
119

Thomas Witkowski's avatar
Thomas Witkowski committed
120
    switch (degree) {
121
122
123
124
125
126
127
128
129
130
131
    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
132
      switch (dim) {
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
      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
152
      switch (dim) {
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
      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
172
      switch (dim) {
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
      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()
  {
198
    if (static_cast<int>(baryDimDegree[dim][degree].size()) == 0) {
199
      ndofDimDegree[dim][degree] = new DimVec<int>(dim, DEFAULT_VALUE, 0);
200

Thomas Witkowski's avatar
Thomas Witkowski committed
201
      if (degree != 0)
202
	(*ndofDimDegree[dim][degree])[VERTEX] = 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
203
204
      else
	(*ndofDimDegree[dim][degree])[VERTEX] = 0;      
205

206
      for (int i = 1; i < dim + 1; i++) {
207
	nBasFcts = getNumberOfDofs(i, degree);
208
	(*ndofDimDegree[dim][degree])[INDEX_OF_DIM(i, dim)] = nBasFcts;
209
	for (int j = 0; j < i; j++) {
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
	  (*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
230
231
    FUNCNAME("Lagrange::setVertices()");

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

234
235
    int dimOfPosition = DIM_OF_INDEX(position, dim);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
238
    if (degree == 4 && dimOfPosition == 1) {
239
240
241
242
243
244
245
246
      (*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
247
    } else if (degree==4 && dimOfPosition == 2) {
248
      for (int i = 0; i < dimOfPosition + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
249
	(*vertices)[(i + dimOfPosition*nodeIndex) % (dimOfPosition + 1)] = 
250
251
252
253
254
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }    
    } else {
255
      for (int i = 0; i < dimOfPosition + 1; i++) {
Thomas Witkowski's avatar
Thomas Witkowski committed
256
	(*vertices)[(i + nodeIndex) % (dimOfPosition + 1)] = 
257
258
259
260
261
262
263
	  Global::getReferenceElement(dim)->getVertexOfPosition(position,
								positionIndex,
								i);
      }
    }
  }

264
  Lagrange::Phi::Phi(Lagrange* owner, 
265
266
267
		     GeoIndex position, 
		     int positionIndex, 
		     int nodeIndex)
268
    : vertices(NULL)
269
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
270
271
272
273
274
275
276
277
278
    FUNCNAME("Lagrange::Phi::Phi()");

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

389
  Lagrange::Phi::~Phi() { 
390
    delete [] vertices; 
391
  }
392
393


394
  Lagrange::GrdPhi::GrdPhi(Lagrange* owner, 
395
396
397
			   GeoIndex position, 
			   int positionIndex, 
			   int nodeIndex)
398
    : vertices(NULL)
399
400
401
402
403
404
405
406
407
408
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
409
    switch (owner->getDegree()) {
410
    case 0:
411
      switch (position) {
412
413
414
415
416
417
418
419
      case CENTER:
	func = grdPhi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
420
      switch (position) {
421
422
423
424
425
426
427
428
      case VERTEX:
	func = grdPhi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
429
      switch (position) {
430
431
432
433
434
435
436
      case VERTEX:
	func = grdPhi2v;
	break;
      case EDGE:
	func = grdPhi2e;
	break;
      case CENTER:
437
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
438
439
440
441
442
443
444
	func = grdPhi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
445
      switch (position) {
446
447
448
449
450
451
452
      case VERTEX:
	func = grdPhi3v;
	break;
      case EDGE:
	func = grdPhi3e;
	break;
      case FACE:
453
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
454
455
456
	func = grdPhi3f;
	break;
      case CENTER:
457
	switch (owner->getDim()) {
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
	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:
474
      switch (position) {
475
476
477
478
      case VERTEX:
	func = grdPhi4v;
	break;
      case EDGE:
479
	if (nodeIndex == 1)
480
481
482
483
484
	  func = grdPhi4e1;
	else 
	  func = grdPhi4e0;
	break;
      case FACE:
485
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
486
487
488
	func = grdPhi4f;      
	break;
      case CENTER:
489
	switch (owner->getDim()) {
490
	case 1:
491
	  if (nodeIndex == 1)
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
	    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
516
517
  Lagrange::GrdPhi::~GrdPhi() 
  { 
518
    delete [] vertices; 
519
  }
520

521
  Lagrange::D2Phi::D2Phi(Lagrange* owner, 
522
523
524
			 GeoIndex position, 
			 int positionIndex, 
			 int nodeIndex)
525
    : vertices(NULL)
526
527
528
529
530
531
532
533
534
535
  {
    // get relevant vertices
    Lagrange::setVertices(owner->getDim(), 
			  owner->getDegree(), 
			  position, 
			  positionIndex, 
			  nodeIndex, 
			  &vertices);

    // set function pointer
536
    switch (owner->getDegree()) {
537
    case 0:
538
      switch (position) {
539
540
541
542
543
544
545
546
      case CENTER:
	func = D2Phi0c;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 1:
547
      switch (position) {
548
549
550
551
552
553
554
555
      case VERTEX:
	func = D2Phi1v;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 2:
556
      switch (position) {
557
558
559
560
      case VERTEX:
	func = D2Phi2v;
	break;
      case EDGE:
561
	TEST_EXIT_DBG(owner->getDim() > 1)("no edge in 1d\n");
562
563
564
	func = D2Phi2e;
	break;
      case CENTER:
565
	TEST_EXIT_DBG(owner->getDim() == 1)("no center dofs for dim != 1\n");
566
567
568
569
570
571
572
	func = D2Phi2e;
	break;
      default:
	ERROR_EXIT("invalid position\n");
      }
      break;
    case 3:
573
      switch (position) {
574
575
576
577
578
579
580
      case VERTEX:
	func = D2Phi3v;
	break;
      case EDGE:
	func = D2Phi3e;
	break;
      case FACE:
581
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
582
583
584
	func = D2Phi3f;
	break;
      case CENTER:
585
	switch (owner->getDim()) {
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
	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:
602
      switch (position) {
603
604
605
606
      case VERTEX:
	func = D2Phi4v;
	break;
      case EDGE:
Thomas Witkowski's avatar
Thomas Witkowski committed
607
	if (nodeIndex == 1)
608
609
610
611
612
	  func = D2Phi4e1;
	else 
	  func = D2Phi4e0;
	break;
      case FACE:
613
	TEST_EXIT_DBG(owner->getDim() >= 3)("no faces in dim < 3\n");
614
615
616
	func = D2Phi4f;      
	break;
      case CENTER:
617
	switch (owner->getDim()) {
618
	case 1:
619
	  if (nodeIndex == 1)
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
	    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
644
645
  Lagrange::D2Phi::~D2Phi() 
  { 
646
    delete [] vertices; 
647
  }
648

Thomas Witkowski's avatar
Thomas Witkowski committed
649
  void Lagrange::createCoords(int* coordInd, int numCoords, int dimIndex, int rest, 
650
651
			      DimVec<double>* vec)
  {
Thomas Witkowski's avatar
Thomas Witkowski committed
652
    if (vec == NULL)
653
      vec = new DimVec<double>(dim, DEFAULT_VALUE, 0.0);
654

655
656
    if (dimIndex == numCoords - 1) {
      (*vec)[coordInd[dimIndex]] = double(rest) / degree;
657
      DimVec<double>* newCoords = new DimVec<double>(*vec);
658
659
      bary->push_back(newCoords);
    } else {
660
      for (int i = rest - 1; i >= 1; i--) {
661
	(*vec)[coordInd[dimIndex]] = double(i) / degree;
662
	createCoords(coordInd, numCoords, dimIndex + 1, rest - i, vec);
663
664
665
666
667
668
      }
    }
  }

  void Lagrange::setBary()
  {
669
    bary = &baryDimDegree[dim][degree];
670
671
672

    if (static_cast<int>(bary->size()) == 0) {
      for (int i = 0; i <= dim; i++) { // for all positions
673
	int partsAtPos = Global::getGeo(INDEX_OF_DIM(i, dim), dim);
674
	for (int j = 0; j < partsAtPos; j++) { // for all vertices/edges/faces/...
Thomas Witkowski's avatar
Thomas Witkowski committed
675
	  int *coordInd = new int[i + 1];      // indices of relevant coords
Thomas Witkowski's avatar
Thomas Witkowski committed
676
	  for (int k = 0; k < i + 1; k++)
677
	    coordInd[k] = Global::getReferenceElement(dim)->
Thomas Witkowski's avatar
Thomas Witkowski committed
678
	      getVertexOfPosition(INDEX_OF_DIM(i, dim), j, k);
Thomas Witkowski's avatar
Thomas Witkowski committed
679
	  
680
	  createCoords(coordInd, i + 1, 0, degree);
Thomas Witkowski's avatar
Thomas Witkowski committed
681
682
683
	  delete [] coordInd;
	  if (static_cast<int>(bary->size()) == nBasFcts) 
	    return;
684
685
686
687
688
	}
      }
    }
  }

689
  int Lagrange::getNumberOfDofs(int dim, int degree)
690
691
  {
    int result = 0;
Thomas Witkowski's avatar
Thomas Witkowski committed
692
    for (int i = 0; i <= degree; i++)
693
      result += fac(dim - 1 + i) / (fac(i) * fac(dim - 1));
Thomas Witkowski's avatar
Thomas Witkowski committed
694

695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
    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
723
724
    if (dimOfPosition == 0)
      return &sortedVertex;   
725
726

    // edge
727
    if (dimOfPosition == 1 && degree == 2)
Thomas Witkowski's avatar
Thomas Witkowski committed
728
      return &sortedEdgeDeg2;    
729
730
      
    int vertex[3];
731
    int** dof = const_cast<int**>(el->getDof());
732
    int verticesOfPosition = dimOfPosition + 1;
733

734
735
736
    for (int i = 0; i < verticesOfPosition; i++)
      vertex[i] = Global::getReferenceElement(dim)->
	getVertexOfPosition(position, positionIndex, i);   
Thomas Witkowski's avatar
Thomas Witkowski committed
737
    
738
739
    if (dimOfPosition == 1) {
      if (degree == 3) {
Thomas Witkowski's avatar
Thomas Witkowski committed
740
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
741
	  return sortedEdgeDeg3[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
742
743
	else
	  return sortedEdgeDeg3[1];	
744
      } else { // degree == 4
Thomas Witkowski's avatar
Thomas Witkowski committed
745
	if (dof[vertex[0]][0] < dof[vertex[1]][0])
746
	  return sortedEdgeDeg4[0];
Thomas Witkowski's avatar
Thomas Witkowski committed
747
748
	else
	  return sortedEdgeDeg4[1];	
749
750
751
752
      }
    }

    // face
753
754
    if (dimOfPosition == 2) {
      if (degree == 3) {
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
	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
777
    if (dimOfPosition == 3 && degree == 4)
Thomas Witkowski's avatar
Thomas Witkowski committed
778
      return &sortedCenterDeg4;    
779
780
781
782
783

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

Thomas Witkowski's avatar
Thomas Witkowski committed
784
  void Lagrange::getBound(const ElInfo* el_info, BoundaryType* bound) const
785
786
787
788
  {
    el_info->testFlag(Mesh::FILL_BOUND);

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

795
    for (int i = 0; i < dim; i++) {
796
797
      int jto = offset + Global::getGeo(INDEX_OF_DIM(i, dim), dim);
      for (int j = offset; j < jto; j++) {
798
	boundaryType = el_info->getBoundary(j);
799
	int kto = (*nDOF)[INDEX_OF_DIM(i, dim)];
Thomas Witkowski's avatar
Thomas Witkowski committed
800
801
	for (int k = 0; k < kto; k++)
	  bound[index++] = boundaryType;	
802
      }
803
      offset -= Global::getGeo(INDEX_OF_DIM(i + 1, dim), dim);
804
805
806
    }

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

810
    TEST_EXIT_DBG(index == nBasFcts)("found not enough boundarys\n");
811
812
813
814
815
816
817
818
  }


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

821
822
    static double* localVec = NULL;
    static int localVecSize = 0;
823

824
    double *rvec = NULL;
825

826
    if (vec) {
827
828
      rvec = vec;
    } else {
829
830
831
832
833
834
835
836
837
      if (localVec && nBasFcts > localVecSize)  {
	delete [] localVec;
	localVec = new double[nBasFcts]; 
      }
      if (!localVec)
	localVec = new double[nBasFcts]; 

      localVecSize = nBasFcts;
      rvec = localVec;
838
839
    } 

840
    WorldVector<double> x;
841
842

    el_info->testFlag(Mesh::FILL_COORDS);
843
844
845
846
847
 
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	rvec[0] = 0.0;
848
	return(const_cast<const double *>(rvec));
849
      }
850
851
852
853
854
	
      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 {
855
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
856
857
858
859
860
861
862
863
	  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++) {
864
	el_info->coordToWorld(*(*bary)[i], x);
865
866
867
	rvec[i] = (*f)(x);
      }
    }  
868
869
870
871
    
    return(const_cast<const double *>( rvec));
  }

872

873
  const WorldVector<double>* 
874
875
  Lagrange::interpol(const ElInfo *el_info, 
		     int no, 
876
877
878
879
		     const int *b_no,
		     AbstractFunction<WorldVector<double>, WorldVector<double> > *f, 
		     WorldVector<double> *vec)
  {
880
881
    FUNCNAME("*Lagrange::interpol_d()");

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

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

891
892
893
894
895
896
    if (b_no) {
      if (no <= 0  ||  no > getNumber()) {
	ERROR("something is wrong, doing nothing\n");
	for (int i = 0; i < dow; i++)
	  rvec[0][i] = 0.0;
	return(const_cast<const WorldVector<double> *>( rvec));
897
      }
898
899
900
901
      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 {
902
	  el_info->coordToWorld(*(*bary)[b_no[i]], x);
903
904
905
906
907
908
909
	  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++) {
910
	el_info->coordToWorld(*(*bary)[i], x);
911
912
913
	rvec[i] = (*f)(x);
      }
    }  
914
915
916
917
918
919
920
921
922
923
924
925
    
    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;

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

928
    int nrDOFs, n0, node0, num = 0;
929
930
    GeoIndex posIndex;
    DegreeOfFreedom* result;
Thomas Witkowski's avatar
Thomas Witkowski committed
931

932
    if (indices) {
933
934
      result = indices;
    } else {
Thomas Witkowski's avatar
Thomas Witkowski committed
935
      if (localVec && nBasFcts > localVecSize) {
Thomas Witkowski's avatar
Thomas Witkowski committed
936
	delete [] localVec;
Thomas Witkowski's avatar
Thomas Witkowski committed
937
938
939
940
941
	localVec = new DegreeOfFreedom[nBasFcts];
      }
      if (!localVec)
	localVec = new DegreeOfFreedom[nBasFcts];

942
943
944
      localVecSize = nBasFcts;
      result = localVec;
    }
Thomas Witkowski's avatar
Thomas Witkowski committed
945

946
    for (int pos = 0, j = 0; pos <= dim; pos++) {
947
      posIndex = INDEX_OF_DIM(pos, dim);
948
      nrDOFs = admin->getNumberOfDofs(posIndex);
949

Thomas Witkowski's avatar
Thomas Witkowski committed
950
      if (nrDOFs) {
951
	n0 = admin->getNumberOfPreDofs(posIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
952
953
954
955
	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
956
	  const int *indi = orderOfPositionIndices(el, posIndex, i);
957

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

964
965
966
    return result;
  }

Thomas Witkowski's avatar
Thomas Witkowski committed
967
968
969
970
971
972
973
  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);

974
    const DegreeOfFreedom **dof = el->getDof();
Thomas Witkowski's avatar
Thomas Witkowski committed
975
976
977
978
979
    GeoIndex posIndex;
    int nrDOFs, n0, node0, num = 0;

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

      if (nrDOFs) {
983
	n0 = admin->getNumberOfPreDofs(posIndex);
Thomas Witkowski's avatar
Thomas Witkowski committed
984
985
986
987
988
989
990
991
992
993
994
995
996
	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]]);
	}
      }
    }
  }

997
998
999
1000
  void Lagrange::l2ScpFctBas(Quadrature *q,
			     AbstractFunction<WorldVector<double>, WorldVector<double> >* f,
			     DOFVector<WorldVector<double> >* fh)
  {
For faster browsing, not all history is shown. View entire blame