Element.h 16.7 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
// ============================================================================
// ==                                                                        ==
// == 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:
84
    /// private standard constructor because an Element must know his Mesh
85
    Element() {}
86

87
  public:
88
    /// constructs an Element which belongs to Mesh
89
90
    Element(Mesh *);

91
    /// copy constructor
92
93
    Element(const Element& old);

94
    /// destructor
95
96
    virtual ~Element();

97
    ///
98
99
    void deleteElementDOFs();

100
101
102
103
104
105
    /** \brief
     * Clone this Element and return a reference to it. Because also the DOFs
     * are cloned, \ref Mesh::serializedDOfs must be used.
     */
    Element* cloneWithDOFs();

106
107
108
109
110
111
    // ===== getting methods ======================================================

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

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

117
    /// Returns \ref child[1]
118
    inline Element* getSecondChild() const {
119
      return child[1];
120
    }
121

122
    /// Returns \ref child[i], i=0,1
123
    inline Element* getChild(int i) const {
124
      TEST_EXIT_DBG(i==0 || i==1)("i must be 0 or 1\n");
125
      return child[i];
126
    }
127
128
129
130
131

    /** \brief
     * Returns true if Element is a leaf element (\ref child[0] == NULL), returns
     * false otherwise.
     */
132
    inline const bool isLeaf() const { 
133
      return (child[0] == NULL); 
134
    }
135

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

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

146
    /// Returns a pointer to the DOFs of this Element
147
148
    const DegreeOfFreedom** getDOF() const {
      return const_cast<const DegreeOfFreedom**>(dof);
149
    }
150

151
    /// Returns \ref mesh of Element
152
153
    inline Mesh* getMesh() const { 
      return mesh; 
154
    }
155
156
157
158
159
160
161

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

167
168
169
170
	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
      }	
      
      return 0.0;
171
    }
172
173
174
175
176
177

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

183
	return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
184
      }
185
186
      
      return 0.0;
187
    }
188

189
    /// Returns region of element if defined, -1 else.
190
191
    int getRegion() const;

192
    /// Returns local vertex number of the j-th vertex of the i-th edge
193
194
195
196
197
198
199
    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,
200
201
				    int positionIndex,
				    int vertexIndex) const = 0;
202

203
    ///
204
205
    virtual int getPositionOfVertex(int side, int vertex) const = 0;

206
    ///
207
208
    virtual int getEdgeOfFace(int face, int edge) const = 0;

209
    /// Returns the number of parts of type i in this element
210
211
    virtual int getGeo(GeoIndex i) const = 0;

212
    /// Returns Element's \ref mark
213
214
    inline const signed char getMark() const { 
      return mark;
215
    }
216

217
218
219
220
221
    /// Returns \ref newCoord[i]
    inline double getNewCoord(int i) const {
	TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
	return (*newCoord)[i];
    }
222

223
    /// Returns Element's \ref index
224
225
    inline int getIndex() const { 
      return index; 
226
    }
227

228
    /// Returns \ref newCoord
229
230
    inline WorldVector<double>* getNewCoord() const { 
      return newCoord; 
231
    }
232
233
234
235
236
237
238
239
240

    /** \} */

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

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

241
    /// Sets \ref child[0]
242
243
    virtual void setFirstChild(Element *aChild) {
      child[0] = aChild;
244
    }
245

246
    /// Sets \ref child[1]
247
248
    virtual void setSecondChild(Element *aChild) {
      child[1] = aChild;
249
    }
250

251
    /// Sets \ref elementData of Element
252
253
    void setElementData(ElementData* ed) {
      elementData = ed;
254
    }
255
256
257
258
259

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

264
    /// Sets \ref mesh.
265
266
    inline void setMesh(Mesh *m) {
      mesh = m;
267
    }
268

269
    /// Sets the pointer to the DOFs of the i-th node of Element
270
271
272
    DegreeOfFreedom* setDOF(int pos, DegreeOfFreedom* p) {
      dof[pos] = p;
      return dof[pos];
273
    }
274
275
276
277
278
279
280

    /** \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)
    {
281
      if (isLeaf()) {
282
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
283
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
284
	TEST_EXIT_DBG(ld)("leaf data not estimatable\n");
285
286
287

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
288
      } else {
289
290
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
291
    }
292
293
294
295
296
297
298

    /** \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)
    {
299
      if (isLeaf()) {
300
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
301
	ElementData *ld = elementData->getElementData(COARSENABLE);
302
	TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
303
304
305

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
306
      } else {
307
308
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
309
    }
310

311
    /// Sets Elements \ref mark = mark + 1;
312
313
314
    inline void incrementMark() {
      mark++;
    }
315

316
    /// Sets Elements \ref mark = mark - 1;
317
318
319
    inline void decrementMark() {
      if (0 < mark) 
	mark--;
320
    }
321

322
    /// Sets Element's \ref mark
323
324
    inline void setMark(signed char m) {
      mark = m;
325
    }
326
327
328
329
330
331
332
333
334
335
336
337
338
339

    /** \} */

    // ===== 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>& 
340
      sortFaceIndices(int face, FixVec<int,WORLD> *vec) const = 0;
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363

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

364
    /// Returns whether Element is a Line
365
366
    virtual bool isLine() const = 0;

367
    /// Returns whether Element is a Triangle
368
369
    virtual bool isTriangle() const = 0;

370
    /// Returns whether Element is a Tetrahedron
371
372
    virtual bool isTetrahedron() const = 0;

373
    /// Returns whether Element has sideElem as one of its sides.
374
375
376
377
378
379
    virtual bool hasSide(Element *sideElem) const = 0;

    /** \} */

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

380
    /// assignment operator
381
    Element& operator=(const Element& el);
382
383
384
385
386
387
388

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

389
    /// Refines Element's leaf data
390
391
392
393
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) {
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
394
395
396
397
398
	  ElementData *tmp = elementData->getDecorated();
	  DELETE elementData;
	  elementData = tmp;
	}
      }
399
    }
400

401
    /// Coarsens Element's leaf data
402
    inline void coarsenElementData(Element* child1, Element* child2, int elType = 0) {
403
404
      ElementData *childData;
      childData = child1->getElementData();
405
      if (childData) {
406
407
408
409
410
	childData->coarsenElementData(this, child1, child2, elType);
	DELETE childData;
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
411
      if (childData) {
412
413
414
415
	childData->coarsenElementData(this, child2, child1, elType);
	DELETE childData;
	child2->setElementData(NULL);
      }
416
    }
417

418
    /// Returns pointer to \ref elementData
419
420
    inline ElementData* getElementData() const {
      return elementData;
421
    }
422

423
    ///
424
    inline ElementData* getElementData(int typeID) const {
425
      if (elementData) {
426
427
428
	return elementData->getElementData(typeID);
      }
      return NULL;
429
    }
430

431
    /// Deletes the \ref elementData with a specific typeID.
432
    bool deleteElementData(int typeID);
433
434
435
436
437
438
439
440

    /** \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, 
441
			 unsigned char elementTyp = 255);
442

443
    /// Returns whether Element's \ref newCoord is set
444
445
    inline bool isNewCoordSet() const { 
      return (newCoord != NULL);
446
    }
447

448
    /// Frees memory for \ref newCoord
449
    void eraseNewCoord();
450
451
 
    /// Serialize the element to a file.
452
    void serialize(std::ostream &out);
453

454
    /// Deserialize an element from a file.
455
    void deserialize(std::istream &in);
456

457
458
    int calcMemoryUsage();

459
460
    // ===== protected methods ====================================================
  protected:
461
    /// Sets Element's \ref dof pointer. Used by friend class Mesh.
462
463
    void setDOFPtrs();
  
464
    /// Sets Element's \ref index. Used by friend class Mesh.
465
466
    inline void setIndex(int i) {
      index = i;
467
    }
468

469
    /// Used by friend class Mesh while dofCompress
470
471
    void newDOFFct1(const DOFAdmin*);

472
    /// Used by friend class Mesh while dofCompress
473
474
475
476
477
478
479
    void newDOFFct2(const DOFAdmin*);

  protected:
    /** \brief
     * Pointers to the two children of interior elements of the tree. Pointers
     * to NULL for leaf elements.
     */
480
    Element *child[2];
481
482
483
484
485

    /** \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 
486
487
488
     * 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 
489
490
     * barycenter, if present.
     */
491
    DegreeOfFreedom **dof;
492
493
494
495
496
497

    /** \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).
     */
498
    int index;
499
500
501
502
503
504

    /** \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.
     */
505
    signed char mark;
506
507
508
509
510
 
    /** \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 
511
512
     * information can be also produced by the traversal routines in the case of 
     * curved boundary.
513
514
515
     */
    WorldVector<double> *newCoord;

516
    /// Pointer to the Mesh this element belongs to
517
    Mesh* mesh;
518

519
    /// Pointer to Element's leaf data
520
    ElementData* elementData;
521

522
523
524
525
526
527
    /** \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;
528
529
530
531
532
533
534
535

    friend class Mesh;
  };

}

#endif  // AMDIS_ELEMENT_H