Element.h 17.5 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
// ============================================================================
// ==                                                                        ==
// == 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
     */
87
    Element() {}
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
  public:
    /** \brief
     * constructs an Element which belongs to Mesh
     */
    Element(Mesh *);

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

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

104
105
    void deleteElementDOFs();

106
107
108
109
110
111
    /** \brief
     * Clone this Element and return a reference to it. Because also the DOFs
     * are cloned, \ref Mesh::serializedDOfs must be used.
     */
    Element* cloneWithDOFs();

112
113
114
115
116
117
118
119
120
    // ===== getting methods ======================================================

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

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

    /** \brief
     * Returns \ref child[1]
     */
128
    inline Element* getSecondChild() const {
129
      return child[1];
130
    }
131
132
133
134

    /** \brief
     * Returns \ref child[i], i=0,1
     */
135
    inline Element* getChild(int i) const {
136
      TEST_EXIT_DBG(i==0 || i==1)("i must be 0 or 1\n");
137
      return child[i];
138
    }
139
140
141
142
143

    /** \brief
     * Returns true if Element is a leaf element (\ref child[0] == NULL), returns
     * false otherwise.
     */
144
    inline const bool isLeaf() const { 
145
      return (child[0] == NULL); 
146
    }
147
148
149
150

    /** \brief
     * Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
     */
151
152
    const DegreeOfFreedom getDOF(int i, int j) const { 
      return dof[i][j];
153
    }
154
155
156
157

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

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

    /** \brief
     * Returns \ref mesh of Element
     */
172
173
    inline Mesh* getMesh() const { 
      return mesh; 
174
    }
175
176
177
178
179
180
181

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

187
188
189
190
	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
      }	
      
      return 0.0;
191
    }
192
193
194
195
196
197

    /** \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) {
198
      if (isLeaf()) {
199
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
200
	ElementData *ld = elementData->getElementData(COARSENABLE);
201
	TEST_EXIT_DBG(ld)("element data not coarsenable!\n");
202

203
	return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
204
      }
205
206
      
      return 0.0;
207
    }
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223

    /** \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,
224
225
				    int positionIndex,
				    int vertexIndex) const = 0;
226

227
228
229
    /** \brief
     *
     */
230
231
    virtual int getPositionOfVertex(int side, int vertex) const = 0;

232
233
234
    /** \brief
     *
     */
235
236
237
238
239
240
241
242
243
244
245
246
    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;
247
    }
248

249
250
251
252
253
    /// Returns \ref newCoord[i]
    inline double getNewCoord(int i) const {
	TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
	return (*newCoord)[i];
    }
254

255
    /// Returns Element's \ref index
256
257
    inline int getIndex() const { 
      return index; 
258
    }
259

260
    /// Returns \ref newCoord
261
262
    inline WorldVector<double>* getNewCoord() const { 
      return newCoord; 
263
    }
264
265
266
267
268
269
270
271
272
273
274
275

    /** \} */

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

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

    /** \brief
     * Sets \ref child[0]
     */
276
277
    virtual void setFirstChild(Element *aChild) {
      child[0] = aChild;
278
    }
279
280
281
282

    /** \brief
     * Sets \ref child[1]
     */
283
284
    virtual void setSecondChild(Element *aChild) {
      child[1] = aChild;
285
    }
286
287
288
289

    /** \brief
     * Sets \ref elementData of Element
     */
290
291
    void setElementData(ElementData* ed) {
      elementData = ed;
292
    }
293
294
295
296
297

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

    /** \brief
     * Sets \ref mesh.
     */
305
306
    inline void setMesh(Mesh *m) {
      mesh = m;
307
    }
308
309
310
311

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

    /** \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)
    {
323
      if (isLeaf()) {
324
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
325
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
326
	TEST_EXIT_DBG(ld)("leaf data not estimatable\n");
327
328
329

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
330
      } else {
331
332
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
333
    }
334
335
336
337
338
339
340

    /** \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)
    {
341
      if (isLeaf()) {
342
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
343
	ElementData *ld = elementData->getElementData(COARSENABLE);
344
	TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
345
346
347

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
348
      } else {
349
350
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
351
    }
352
353
354
355

    /** \brief
     * Sets Elements \ref mark = mark + 1;
     */
356
357
358
    inline void incrementMark() {
      mark++;
    }
359
360
361
362

    /** \brief
     * Sets Elements \ref mark = mark - 1;
     */
363
364
365
    inline void decrementMark() {
      if (0 < mark) 
	mark--;
366
    }
367
368
369
370

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

    /** \} */

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

    /** \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
     */
439
    Element& operator=(const Element& el);
440
441
442
443
444
445
446
447
448
449

    /** \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
     */
450
451
452
453
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) {
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
454
455
456
457
458
	  ElementData *tmp = elementData->getDecorated();
	  DELETE elementData;
	  elementData = tmp;
	}
      }
459
    }
460
461
462
463

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

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

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

    /** \brief
498
     * Deletes the \ref elementData with a specific typeID.
499
     */
500
    bool deleteElementData(int typeID);
501
502
503
504
505
506
507
508

    /** \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, 
509
			 unsigned char elementTyp = 255);
510
511
512
513

    /** \brief
     * Returns whether Element's \ref newCoord is set
     */
514
515
    inline bool isNewCoordSet() const { 
      return (newCoord != NULL);
516
    }
517
518
519
520
521
522
523
524

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

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

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

529
530
    int calcMemoryUsage();

531
532
533
534
535
536
537
538
539
540
    // ===== 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.
     */
541
542
    inline void setIndex(int i) {
      index = i;
543
    }
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559

    /** \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.
     */
560
    Element *child[2];
561
562
563
564
565

    /** \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 
566
567
568
     * 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 
569
570
     * barycenter, if present.
     */
571
    DegreeOfFreedom **dof;
572
573
574
575
576
577

    /** \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).
     */
578
    int index;
579
580
581
582
583
584

    /** \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.
     */
585
    signed char mark;
586
587
588
589
590
 
    /** \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 
591
592
     * information can be also produced by the traversal routines in the case of 
     * curved boundary.
593
594
595
596
597
598
     */
    WorldVector<double> *newCoord;

    /** \brief
     * Pointer to the Mesh this element belongs to
     */
599
    Mesh* mesh;
600
601
602
603

    /** \brief
     * Pointer to Element's leaf data
     */
604
    ElementData* elementData;
605
606


607
608
609
610
611
612
    /** \brief
     * This map is used for deletion of all DOFs of all elements of a mesh. Once
     * a DOF-vector (all DOFS at a node, edge, etc.) is deleted, its address is
     * added to this map to note not to delete it a second time.
     */
    static std::map<DegreeOfFreedom*, bool> deletedDOFs;
613
614
615
616
617
618
619
620

    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H