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
254
255
256
257
258

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

    /** \brief
     * Returns Element's \ref index
     */
    inline int getIndex() const { 
      return index; 
259
    }
260
261
262
263
264
265

    /** \brief
     * Returns \ref newCoord
     */
    inline WorldVector<double>* getNewCoord() const { 
      return newCoord; 
266
    }
267
268
269
270
271
272
273
274
275
276
277
278

    /** \} */

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

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

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

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

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

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

    /** \brief
     * Sets \ref mesh.
     */
308
309
    inline void setMesh(Mesh *m) {
      mesh = m;
310
    }
311
312
313
314

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

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

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

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

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

    /** \brief
     * Sets Elements \ref mark = mark + 1;
     */
359
360
361
    inline void incrementMark() {
      mark++;
    }
362
363
364
365

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

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

    /** \} */

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

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

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

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

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

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

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

    /** \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, 
512
			 unsigned char elementTyp = 255);
513
514
515
516

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

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

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

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

532
533
    int calcMemoryUsage();

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

    /** \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.
     */
563
    Element *child[2];
564
565
566
567
568

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

    /** \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).
     */
581
    int index;
582
583
584
585
586
587

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

    /** \brief
     * Pointer to the Mesh this element belongs to
     */
602
    Mesh* mesh;
603
604
605
606

    /** \brief
     * Pointer to Element's leaf data
     */
607
    ElementData* elementData;
608
609


610
611
612
613
614
615
    /** \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;
616
617
618
619
620
621
622
623

    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H