Element.h 17.4 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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
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
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
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
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
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
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
// ============================================================================
// ==                                                                        ==
// == 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 {
      FUNCNAME("Element::getChild");
      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.
     */
    inline const bool isLeaf() const { return (child[0]==NULL); };

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

    /** \brief
     * Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
     */
    const DegreeOfFreedom* getDOF(int i) const {return dof[i];};

    /** \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
     */
    inline Mesh* getMesh() const { return mesh; };

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

	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  getErrorEstimate(row);
      }
      else return 0.0;
    };

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

	return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  getCoarseningErrorEstimate(row);
      }
      else return 0.0;
    };

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

    virtual int getPositionOfVertex(int side, int vertex) const = 0;

    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]
     */
    virtual void setFirstChild(Element *aChild) { child[0] = aChild; };

    /** \brief
     * Sets \ref child[1]
     */
    virtual void setSecondChild(Element *aChild) { child[1] = aChild; };

    /** \brief
     * Sets \ref elementData of Element
     */
    void setElementData(ElementData* ed) { elementData = ed; };

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

    /** \brief
     * Sets \ref mesh.
     */
    inline void setMesh(Mesh *m) { mesh = m; };

    /** \brief
     * Sets the pointer to the DOFs of the i-th node of Element
     */
    DegreeOfFreedom* setDOF(int i,DegreeOfFreedom* p) {dof[i]=p;return dof[i];};

    /** \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)
    {
      if(isLeaf()) {
	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);
      }
      else {
	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)
    {
      if(isLeaf()) {
	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);
      }
      else {
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
    };

    /** \brief
     * Sets Elements \ref mark = mark + 1;
     */
    inline void incrementMark() {mark++;}

    /** \brief
     * Sets Elements \ref mark = mark - 1;
     */
    inline void decrementMark() {if (0<mark) mark--;};

    /** \brief
     * Sets Element's \ref mark
     */
    inline void setMark(signed char m) {mark=m;};

    /** \} */

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

    /** \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
     */
    inline void refineElementData(Element* child1, Element* child2, int elType=0) {
      if(elementData) {
	bool remove = 
	  elementData->refineElementData(this, child1, child2, elType);
	if(remove) {
	  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();
      if(childData) {
	childData->coarsenElementData(this, child1, child2, elType);
	DELETE childData;
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
      if(childData) {
	childData->coarsenElementData(this, child2, child1, elType);
	DELETE childData;
	child2->setElementData(NULL);
      }
    };

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

    inline ElementData* getElementData(int typeID) const {
      if(elementData) {
	return elementData->getElementData(typeID);
      }
      return NULL;
    };

    /** \brief
     * kills \ref elementData
     */
    bool deleteElementData(int typeID) {
      FUNCNAME("Element::deleteElementData()");
      if(elementData) {
	if(elementData->isOfType(typeID)) {
	  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, 
			 unsigned char elementTyp=255);

    /** \brief
     * Returns whether Element's \ref newCoord is set
     */
    inline bool isNewCoordSet() const { return (newCoord != NULL);};

    /** \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.
     */
    inline void setIndex(int i) {index=i ; };

    /** \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.
     */
    Element          *child[2];

    /** \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).
     */
    int              index;

    /** \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.
     */
    signed char      mark;
 
    /** \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
     */
    Mesh*            mesh;

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

    //  struct ElementConnection
    //  {
    //    Element *connectedElement;
    //    int connectedToPart;   // number of edge/face of this element
    //    ElementConnection* virtualConnection;
    //  };

    //  ::std::list<ElementConnection*> connections;


    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H