Element.h 17.3 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
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  crystal growth group                                                  ==
// ==                                                                        ==
// ==  Stiftung caesar                                                       ==
// ==  Ludwig-Erhard-Allee 2                                                 ==
// ==  53175 Bonn                                                            ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  http://www.caesar.de/cg/AMDiS                                         ==
// ==                                                                        ==
// ============================================================================

/** \file Element.h */

#ifndef AMDIS_ELEMENT_H
#define AMDIS_ELEMENT_H

// ============================================================================
// ===== includes =============================================================
// ============================================================================

#include "Global.h"
#include "RefinementManager.h"
#include "Serializable.h"
#include "ElementData.h"
#include "LeafData.h"

namespace AMDiS {

  // ============================================================================
  // ===== forward declarations =================================================
  // ============================================================================

  class Mesh;
  class DOFAdmin;
  template<typename T> class WorldVector;
  class CoarseningManager;

  template<typename T, GeoIndex d> class FixVec;

#define AMDIS_UNDEFINED  5

  // ============================================================================
  // ===== class Element ========================================================
  // ============================================================================

  /** \ingroup Triangulation 
   * \brief
   * Base class for Line, Triangle, Tetrahedron
   *
   * Elements in AMDiS are always simplices (a simplex is a Line in 1d, a 
   * Triangle in 2d and a Tetrahedron in 3d). 
   * We restrict ourselves here to simplicial meshes, for several reasons:
   * -# A simplex is one of the most simple geometric types and complex domains 
   *    may be approximated by a set of simplices quite easily.
   * -# Simplicial meshes allow local refinement without the need of 
   *    nonconforming meshes (hanging nodes), parametric elements, or mixture of
   *    element types (which is the case for quadrilateral meshes).
   * -# Polynomials of any degree are easily represented on a simplex using 
   *    local (barycentric) coordinates.
   *
   * A Line element and its refinement:
   *
   * <img src="line.png">
   *
   * A Triangle element and its refinement:
   *
   * <img src="triangle.png">
   *
   * A Tetrahedron element and its refinements:
   *
   * <img src="tetrahedron.png">
   */
  class Element : public Serializable
  {
  private:
    /** \brief
     * private standard constructor because an Element must know his Mesh
     */
    Element() {};
  public:
    /** \brief
     * constructs an Element which belongs to Mesh
     */
    Element(Mesh *);

    /** \brief
     * copy constructor
     */
    Element(const Element& old);

    /** \brief
     * destructor
     */ 
    virtual ~Element();

    // ===== getting methods ======================================================

    /** \name getting methods
     * \{
     */

    /** \brief
     * Returns \ref child[0]
     */
113
    inline Element* getFirstChild() const {
114
115
116
117
118
119
      return child[0];
    };

    /** \brief
     * Returns \ref child[1]
     */
120
    inline Element* getSecondChild() const {
121
122
123
124
125
126
      return child[1];
    };

    /** \brief
     * Returns \ref child[i], i=0,1
     */
127
    inline Element* getChild(int i) const {
128
      TEST_EXIT_DBG(i==0 || i==1)("i must be 0 or 1\n");
129
130
131
132
133
134
135
      return child[i];
    };

    /** \brief
     * Returns true if Element is a leaf element (\ref child[0] == NULL), returns
     * false otherwise.
     */
136
    inline const bool isLeaf() const { 
137
      return (child[0] == NULL); 
138
    };
139
140
141
142

    /** \brief
     * Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
     */
143
144
145
    const DegreeOfFreedom getDOF(int i, int j) const { 
      return dof[i][j];
    };
146
147
148
149

    /** \brief
     * Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
     */
150
151
152
    const DegreeOfFreedom* getDOF(int i) const {
      return dof[i];
    };
153
154
155
156
157
158
159
160
161
162
163

    /** \brief
     * Returns a pointer to the DOFs of this Element
     */
    const DegreeOfFreedom** getDOF() const {
      return const_cast<const DegreeOfFreedom**>(dof);
    };

    /** \brief
     * Returns \ref mesh of Element
     */
164
165
166
    inline Mesh* getMesh() const { 
      return mesh; 
    };
167
168
169
170
171
172
173

    /** \brief
     * Returns \ref elementData's error estimation, if Element is a leaf element
     * and has leaf data. 
     */
    inline double getEstimation(int row) const
    {
174
      if (isLeaf()) {
175
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
176
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
177
	TEST_EXIT_DBG(ld)("leaf data not estimatable!\n");
178

179
180
181
182
	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
      }	
      
      return 0.0;
183
184
185
186
187
188
189
    };

    /** \brief
     * Returns Element's coarsening error estimation, if Element is a leaf 
     * element and if it has leaf data and if this leaf data are coarsenable.
     */
    inline double getCoarseningEstimation(int row) {
190
      if (isLeaf()) {
191
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
192
	ElementData *ld = elementData->getElementData(COARSENABLE);
193
	TEST_EXIT_DBG(ld)("element data not coarsenable!\n");
194

195
	return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
196
      }
197
198
      
      return 0.0;
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
    };

    /** \brief
     * Returns region of element if defined, -1 else.
     */
    int getRegion() const;

    /** \brief
     * Returns local vertex number of the j-th vertex of the i-th edge
     */
    virtual int getVertexOfEdge(int i, int j) const = 0; 

    /** \brief
     * Returns local vertex number of the vertexIndex-th vertex of the
     * positionIndex-th part of type position (vertex, edge, face)
     */
    virtual int getVertexOfPosition(GeoIndex position,
				    int      positionIndex,
				    int      vertexIndex) const = 0;

219
220
221
    /** \brief
     *
     */
222
223
    virtual int getPositionOfVertex(int side, int vertex) const = 0;

224
225
226
    /** \brief
     *
     */
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
    virtual int getEdgeOfFace(int face, int edge) const = 0;

    /** \brief
     * Returns the number of parts of type i in this element
     */
    virtual int getGeo(GeoIndex i) const = 0;

    /** \brief
     * Returns Element's \ref mark
     */
    inline const signed char getMark() const { 
      return mark;
    };

    /** \brief
     * Returns \ref newCoord[i]
     */
    double getNewCoord(int j) const;

    /** \brief
     * Returns Element's \ref index
     */
    inline int getIndex() const { 
      return index; 
    };

    /** \brief
     * Returns \ref newCoord
     */
    inline WorldVector<double>* getNewCoord() const { 
      return newCoord; 
    };

    /** \} */

    // ===== setting methods ======================================================

    /** \name setting methods
     * \{
     */

    /** \brief
     * Sets \ref child[0]
     */
271
272
273
    virtual void setFirstChild(Element *aChild) {
      child[0] = aChild;
    };
274
275
276
277

    /** \brief
     * Sets \ref child[1]
     */
278
279
280
    virtual void setSecondChild(Element *aChild) {
      child[1] = aChild;
    };
281
282
283
284

    /** \brief
     * Sets \ref elementData of Element
     */
285
286
287
    void setElementData(ElementData* ed) {
      elementData = ed;
    };
288
289
290
291
292

    /** \brief
     * Sets \ref newCoord of Element. Needed by refinement, if Element has a
     * boundary edge on a curved boundary.
     */
293
294
295
    inline void setNewCoord(WorldVector<double>* coord) {
      newCoord = coord;
    };
296
297
298
299

    /** \brief
     * Sets \ref mesh.
     */
300
301
302
    inline void setMesh(Mesh *m) {
      mesh = m;
    };
303
304
305
306

    /** \brief
     * Sets the pointer to the DOFs of the i-th node of Element
     */
307
308
309
    DegreeOfFreedom* setDOF(int pos, DegreeOfFreedom* p) {
      dof[pos] = p;
      return dof[pos];
310
    };
311
312
313
314
315
316
317

    /** \brief
     * Checks whether Element is a leaf element and whether it has leaf data.
     * If the checks don't fail, leaf data's error estimation is set to est.
     */
    inline void setEstimation(double est, int row)
    {
318
      if (isLeaf()) {
319
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
320
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
321
	TEST_EXIT_DBG(ld)("leaf data not estimatable\n");
322
323
324

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
325
      } else {
326
327
328
329
330
331
332
333
334
335
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
    };

    /** \brief
     * Sets Element's coarsening error estimation, if Element is a leaf element
     * and if it has leaf data and if this leaf data are coarsenable.
     */
    inline void setCoarseningEstimation(double est, int row)
    {
336
      if (isLeaf()) {
337
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
338
	ElementData *ld = elementData->getElementData(COARSENABLE);
339
	TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
340
341
342

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
343
      } else {
344
345
346
347
348
349
350
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
    };

    /** \brief
     * Sets Elements \ref mark = mark + 1;
     */
351
352
353
    inline void incrementMark() {
      mark++;
    }
354
355
356
357

    /** \brief
     * Sets Elements \ref mark = mark - 1;
     */
358
359
360
361
    inline void decrementMark() {
      if (0 < mark) 
	mark--;
    };
362
363
364
365

    /** \brief
     * Sets Element's \ref mark
     */
366
367
368
    inline void setMark(signed char m) {
      mark = m;
    };
369
370
371
372
373
374
375
376
377
378
379
380
381
382

    /** \} */

    // ===== pure virtual methods =================================================

    /** \name pure virtual methods 
     * \{ 
     */

    /** \brief
     * orient the vertices of edges/faces.
     * Used by Estimator for the jumps => same quadrature nodes from both sides!
     */
    virtual const FixVec<int,WORLD>& 
383
      sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444

    /** \brief
     * Returns a copy of itself. Needed by Mesh to create Elements by a
     * prototype. 
     */ 
    virtual Element *clone() = 0;

    /** \brief
     * Returns which side of child[childnr] corresponds to side sidenr of 
     * this Element. If the child has no corresponding
     * side, the return value is negative. *isBisected is true after the
     * function call, if the side of the child is only a part of element's 
     * side, false otherwise. 
     */
    virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0;

    /** \brief
     * Returns which vertex of elements parent corresponds to the vertexnr of
     * the element, if the element is the childnr-th child of the parent.
     * If the vertex is the ner vertex at the refinement edge, -1 is returned.
     */
    virtual int getVertexOfParent(int childnr, int vertexnr, int elType = 0) const = 0;

    /** \brief
     * Returns whether Element is a Line
     */
    virtual bool isLine() const = 0;

    /** \brief
     * Returns whether Element is a Triangle
     */
    virtual bool isTriangle() const = 0;

    /** \brief
     * Returns whether Element is a Tetrahedron
     */
    virtual bool isTetrahedron() const = 0;

    /** \brief
     * Returns whether Element has sideElem as one of its sides.
     */
    virtual bool hasSide(Element *sideElem) const = 0;

    /** \} */

    // ===== other public methods =================================================

    /** \brief
     * assignment operator
     */
    Element& operator=(const Element& old);

    /** \brief
     * Checks whether the face with vertices dof[0],..,dof[DIM-1] is
     * part of mel's boundary. returns the opposite vertex if true, -1 else
     */
    int oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const;

    /** \brief
     * Refines Element's leaf data
     */
445
446
447
448
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) {
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
449
450
451
452
453
454
455
456
457
458
459
460
461
	  ElementData *tmp = elementData->getDecorated();
	  DELETE elementData;
	  elementData = tmp;
	}
      }
    };

    /** \brief
     * Coarsens Element's leaf data
     */
    inline void coarsenElementData(Element* child1, Element* child2, int elType=0) {
      ElementData *childData;
      childData = child1->getElementData();
462
      if (childData) {
463
464
465
466
467
	childData->coarsenElementData(this, child1, child2, elType);
	DELETE childData;
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
468
      if (childData) {
469
470
471
472
473
474
475
476
477
478
479
480
481
	childData->coarsenElementData(this, child2, child1, elType);
	DELETE childData;
	child2->setElementData(NULL);
      }
    };

    /** \brief
     * Returns pointer to \ref elementData
     */
    inline ElementData* getElementData() const {
      return elementData;
    };

482
483
484
    /** \brief
     *
     */
485
    inline ElementData* getElementData(int typeID) const {
486
      if (elementData) {
487
488
489
490
491
492
493
494
495
496
	return elementData->getElementData(typeID);
      }
      return NULL;
    };

    /** \brief
     * kills \ref elementData
     */
    bool deleteElementData(int typeID) {
      FUNCNAME("Element::deleteElementData()");
497
498
      if (elementData) {
	if (elementData->isOfType(typeID)) {
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
	  ElementData *tmp = elementData;
	  elementData = elementData->getDecorated();
	  DELETE tmp;
	  return true;
	} else {
	  return elementData->deleteDecorated(typeID);
	}
      }
      return false;
    };

    /** \brief
     * Returns whether element is refined at side side
     * el1, el2 are the corresponding children. 
     * (not neccessarly the direct children!)
     * elementTyp is the type of this element (comes from ElInfo)
     */
    bool isRefinedAtSide(int side, Element *el1, Element *el2, 
517
			 unsigned char elementTyp = 255);
518
519
520
521

    /** \brief
     * Returns whether Element's \ref newCoord is set
     */
522
523
524
    inline bool isNewCoordSet() const { 
      return (newCoord != NULL);
    };
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546

    /** \brief
     * Frees memory for \ref newCoord
     */
    void eraseNewCoord();

    // ===== Serializable implementation =====
  
    void serialize(::std::ostream &out);

    void deserialize(::std::istream &in);

    // ===== protected methods ====================================================
  protected:
    /** \brief
     * Sets Element's \ref dof pointer. Used by friend class Mesh.
     */
    void setDOFPtrs();
  
    /** \brief
     * Sets Element's \ref index. Used by friend class Mesh.
     */
547
548
549
    inline void setIndex(int i) {
      index = i;
    };
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565

    /** \brief
     * Used by friend class Mesh while dofCompress
     */
    void newDOFFct1(const DOFAdmin*);

    /** \brief
     * Used by friend class Mesh while dofCompress
     */
    void newDOFFct2(const DOFAdmin*);

  protected:
    /** \brief
     * Pointers to the two children of interior elements of the tree. Pointers
     * to NULL for leaf elements.
     */
566
    Element *child[2];
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584

    /** \brief
     * Vector of pointers to DOFs. These pointers must be available for elements
     * vertices (for the geometric description of the mesh). There my be pointers
     * for the edges, for faces and for the center of an element. They are 
     * ordered
     * the following way: The first N_VERTICES entries correspond to the DOFs at
     * the vertices of the element. The next ones are those at the edges, if 
     * present, then those at the faces, if present, and then those at the 
     * barycenter, if present.
     */
    DegreeOfFreedom  **dof;

    /** \brief
     * Unique global index of the element. these indices are not strictly ordered
     * and may be larger than the number of elements in the binary tree (the list
     * of indices may have holes after coarsening).
     */
585
    int index;
586
587
588
589
590
591

    /** \brief
     * Marker for refinement and coarsening. if mark is positive for a leaf
     * element, this element is refined mark times. if mark is negative for
     * a leaf element, this element is coarsened -mark times.
     */
592
    signed char mark;
593
594
595
596
597
598
599
600
601
602
603
604
605
606
 
    /** \brief
     * If the element has a boundary edge on a curved boundary, this is a pointer
     * to the coordinates of the new vertex that is created due to the refinement
     * of the element, otherwise it is a NULL pointer. Thus coordinate 
     * information
     * can be also produced by the traversal routines in the case of curved 
     * boundary.
     */
    WorldVector<double> *newCoord;

    /** \brief
     * Pointer to the Mesh this element belongs to
     */
607
    Mesh* mesh;
608
609
610
611

    /** \brief
     * Pointer to Element's leaf data
     */
612
    ElementData* elementData;
613
614
615
616
617
618
619
620
621
622



    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H