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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
// ============================================================================
// ==                                                                        ==
// == 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]
     */
    virtual Element* getFirstChild() const {
      return child[0];
    };

    /** \brief
     * Returns \ref child[1]
     */
    virtual Element* getSecondChild() const {
      return child[1];
    };

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

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

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

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

    /** \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
     */
165
166
167
    inline Mesh* getMesh() const { 
      return mesh; 
    };
168
169
170
171
172
173
174

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

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

    /** \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) {
191
      if (isLeaf()) {
192
193
194
195
	TEST_EXIT(elementData)("leaf element without leaf data\n");
	ElementData *ld = elementData->getElementData(COARSENABLE);
	TEST_EXIT(ld)("element data not coarsenable!\n");

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

    /** \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;

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

225
226
227
    /** \brief
     *
     */
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
271
    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]
     */
272
273
274
    virtual void setFirstChild(Element *aChild) {
      child[0] = aChild;
    };
275
276
277
278

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

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

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

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

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

    /** \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)
    {
319
      if (isLeaf()) {
320
321
322
323
324
325
	TEST_EXIT(elementData)("leaf element without leaf data\n");
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
	TEST_EXIT(ld)("leaf data not estimatable\n");

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
326
      } else {
327
328
329
330
331
332
333
334
335
336
	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)
    {
337
      if (isLeaf()) {
338
339
340
341
342
343
	TEST_EXIT(elementData)("leaf element without leaf data\n");
	ElementData *ld = elementData->getElementData(COARSENABLE);
	TEST_EXIT(ld)("leaf data not coarsenable\n");

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

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

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

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

    /** \} */

    // ===== 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>& 
384
      sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
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
445

    /** \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
     */
446
447
448
449
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) {
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
450
451
452
453
454
455
456
457
458
459
460
461
462
	  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();
463
      if (childData) {
464
465
466
467
468
	childData->coarsenElementData(this, child1, child2, elType);
	DELETE childData;
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
469
      if (childData) {
470
471
472
473
474
475
476
477
478
479
480
481
482
	childData->coarsenElementData(this, child2, child1, elType);
	DELETE childData;
	child2->setElementData(NULL);
      }
    };

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

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

    /** \brief
     * kills \ref elementData
     */
    bool deleteElementData(int typeID) {
      FUNCNAME("Element::deleteElementData()");
498
499
      if (elementData) {
	if (elementData->isOfType(typeID)) {
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
	  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, 
518
			 unsigned char elementTyp = 255);
519
520
521
522

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

    /** \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.
     */
548
549
550
    inline void setIndex(int i) {
      index = i;
    };
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566

    /** \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.
     */
567
    Element *child[2];
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585

    /** \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).
     */
586
    int index;
587
588
589
590
591
592

    /** \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.
     */
593
    signed char mark;
594
595
596
597
598
599
600
601
602
603
604
605
606
607
 
    /** \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
     */
608
    Mesh* mesh;
609
610
611
612

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



    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H