Element.h 18.9 KB
Newer Older
1
2
3
4
5
6
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
7
// ==  TU Dresden                                                            ==
8
// ==                                                                        ==
9
10
11
// ==  Institut fr Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
12
13
14
15
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
16
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
17
18
19
20
21
22
23
24
// ==                                                                        ==
// ============================================================================

/** \file Element.h */

#ifndef AMDIS_ELEMENT_H
#define AMDIS_ELEMENT_H

25
#include "AMDiS_fwd.h"
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
#include "Global.h"
#include "RefinementManager.h"
#include "Serializable.h"
#include "ElementData.h"
#include "LeafData.h"

namespace AMDiS {

  template<typename T, GeoIndex d> class FixVec;

#define AMDIS_UNDEFINED  5

  /** \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:
68
    /// private standard constructor because an Element must know his Mesh
69
    Element() {}
70

71
  public:
72
    /// constructs an Element which belongs to Mesh
73
74
    Element(Mesh *);

75
    /// copy constructor
76
77
    Element(const Element& old);

78
    /// destructor
79
80
    virtual ~Element();

81
    ///
82
83
    void deleteElementDOFs();

84
85
86
87
88
89
    /** \brief
     * Clone this Element and return a reference to it. Because also the DOFs
     * are cloned, \ref Mesh::serializedDOfs must be used.
     */
    Element* cloneWithDOFs();

90
91
92
93
    /** \name getting methods
     * \{
     */

94
    /// Returns \ref child[0]
Thomas Witkowski's avatar
Thomas Witkowski committed
95
96
    inline Element* getFirstChild() const 
    {
97
      return child[0];
98
    }
99

100
    /// Returns \ref child[1]
Thomas Witkowski's avatar
Thomas Witkowski committed
101
102
    inline Element* getSecondChild() const 
    {
103
      return child[1];
104
    }
105

106
    /// Returns \ref child[i], i=0,1
Thomas Witkowski's avatar
Thomas Witkowski committed
107
108
    inline Element* getChild(int i) const 
    {
Thomas Witkowski's avatar
Thomas Witkowski committed
109
      TEST_EXIT_DBG(i == 0 || i == 1)("There is only child 0 or 1! (i = %d)\n", i);
110
      return child[i];
111
    }
112
113
114
115
116

    /** \brief
     * Returns true if Element is a leaf element (\ref child[0] == NULL), returns
     * false otherwise.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
117
118
    inline const bool isLeaf() const 
    { 
119
      return (child[0] == NULL); 
120
    }
121

122
    /// Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
123
    inline DegreeOfFreedom getDof(int i, int j) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
124
    { 
125
126
      TEST_EXIT_DBG(dofValid)("Should not happen!\n");

127
      return dof[i][j];
128
    }
129

130
    /// Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
131
    inline const DegreeOfFreedom* getDof(int i) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
132
    {
133
134
      TEST_EXIT_DBG(dofValid)("Should not happen!\n");

135
      return dof[i];
136
    }
137

138
    /// Returns a pointer to the DOFs of this Element
139
    inline const DegreeOfFreedom** getDof() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
140
    {
141
142
      TEST_EXIT_DBG(dofValid)("Should not happen!\n");

143
      return const_cast<const DegreeOfFreedom**>(dof);
144
    }
145

146
    /// Returns \ref mesh of Element
Thomas Witkowski's avatar
Thomas Witkowski committed
147
148
    inline Mesh* getMesh() const 
    { 
149
      return mesh; 
150
    }
151

152
153
154
155
156
    inline void setDofValid(bool b)
    {
      dofValid = b;
    }

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

168
169
170
171
	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
      }	
      
      return 0.0;
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.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
178
179
    inline double getCoarseningEstimation(int row) 
    {
180
      if (isLeaf()) {
181
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
182
	ElementData *ld = elementData->getElementData(COARSENABLE);
183
	TEST_EXIT_DBG(ld)("element data not coarsenable!\n");
184

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

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

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

205
    ///
206
207
    virtual int getPositionOfVertex(int side, int vertex) const = 0;

208
    ///
209
210
    virtual int getEdgeOfFace(int face, int edge) const = 0;

211
    ///
212
213
214
215
    virtual DofEdge getEdge(int localEdgeIndex) const = 0;

    ///
    virtual DofFace getFace(int localFaceIndex) const = 0;
216

217
    /// Returns the number of parts of type i in this element
218
219
    virtual int getGeo(GeoIndex i) const = 0;

220
    /// Returns Element's \ref mark
221
    inline const int getMark() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
222
    { 
223
      return mark;
224
    }
225

226
    /// Returns \ref newCoord[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
227
228
    inline double getNewCoord(int i) const 
    {
229
230
231
	TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
	return (*newCoord)[i];
    }
232

233
    /// Returns Element's \ref index
Thomas Witkowski's avatar
Thomas Witkowski committed
234
235
    inline int getIndex() const 
    { 
236
      return index; 
237
    }
238

239
    /// Returns \ref newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
240
241
    inline WorldVector<double>* getNewCoord() const 
    { 
242
      return newCoord; 
243
    }
244
245
246
247
248
249
250

    /** \} */

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

251
    /// Sets \ref child[0]
Thomas Witkowski's avatar
Thomas Witkowski committed
252
253
    virtual void setFirstChild(Element *aChild) 
    {
254
      child[0] = aChild;
255
    }
256

257
    /// Sets \ref child[1]
Thomas Witkowski's avatar
Thomas Witkowski committed
258
259
    virtual void setSecondChild(Element *aChild) 
    {
260
      child[1] = aChild;
261
    }
262

263
    /// Sets \ref elementData of Element
Thomas Witkowski's avatar
Thomas Witkowski committed
264
265
    void setElementData(ElementData* ed) 
    {
266
      elementData = ed;
267
    }
268
269
270
271
272

    /** \brief
     * Sets \ref newCoord of Element. Needed by refinement, if Element has a
     * boundary edge on a curved boundary.
     */
Thomas Witkowski's avatar
Thomas Witkowski committed
273
274
    inline void setNewCoord(WorldVector<double>* coord) 
    {
275
      newCoord = coord;
276
    }
277

278
    /// Sets \ref mesh.
Thomas Witkowski's avatar
Thomas Witkowski committed
279
280
    inline void setMesh(Mesh *m) 
    {
281
      mesh = m;
282
    }
283

284
    /// Sets the pointer to the DOFs of the i-th node of Element
285
    inline void setDof(int pos, DegreeOfFreedom* p) 
Thomas Witkowski's avatar
Thomas Witkowski committed
286
    {
287
      dof[pos] = p;
288
    }
289
290
291
292
293
294
295

    /** \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)
    {
296
297
      FUNCNAME("Element::setEstimation()");

298
      if (isLeaf()) {
299
	TEST_EXIT_DBG(elementData)("Leaf element %d without leaf data!\n", index);
300
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
301
	TEST_EXIT_DBG(ld)("Leaf data %d not estimatable!\n", index);
302
303
304

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
305
      } else {
306
307
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
308
    }
309
310
311
312
313
314
315

    /** \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)
    {
316
      if (isLeaf()) {
317
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
318
	ElementData *ld = elementData->getElementData(COARSENABLE);
319
	TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
320
321
322

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
323
      } else {
324
325
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
326
    }
327

328
    /// Sets Elements \ref mark = mark + 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
329
330
    inline void incrementMark() 
    {
331
332
      mark++;
    }
333

334
    /// Sets Elements \ref mark = mark - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
335
336
    inline void decrementMark() 
    {
337
338
      if (0 < mark) 
	mark--;
339
    }
340

341
    /// Sets Element's \ref mark
342
    inline void setMark(int m) 
Thomas Witkowski's avatar
Thomas Witkowski committed
343
    {
344
      mark = m;
345
    }
346
347
348
349
350
351
352
353
354
355
356
357
358

    /** \} */


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

361
    /// Returns a copy of itself. Needed by Mesh to create Elements by a prototype. 
362
363
364
    virtual Element *clone() = 0;

    /** \brief
Thomas Witkowski's avatar
Thomas Witkowski committed
365
366
     * 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. 
367
368
369
     */
    virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0;

370
371
372
373
374
375
376
377
378
379
380
381
    /** \brief
     * Generalization of \ref getSideOfChild to arbitrary subObject. Thus, e.g., in 3d
     * we can ask for the local id of a verte, edge or face on the elements children.
     *
     * \param[in]  childnr    Either 0 or 1 for the left or right children.
     * \param[in]  subObj     Defines whether we ask for VERTEX, EDGE or FACE.
     * \param[in]  ithObj     Number of the object on the parent.
     * \param[in]  elType     Type of the element. Important only in 3D.
     */
    virtual int getSubObjOfChild(int childnr, GeoIndex subObj, int ithObj, 
				 int elType = 0) const = 0;

382
383
384
385
386
387
388
    /** \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;

389
    /// Returns whether Element is a Line
390
391
    virtual bool isLine() const = 0;

392
    /// Returns whether Element is a Triangle
393
394
    virtual bool isTriangle() const = 0;

395
    /// Returns whether Element is a Tetrahedron
396
397
    virtual bool isTetrahedron() const = 0;

398
    /// Returns whether Element has sideElem as one of its sides.
399
400
    virtual bool hasSide(Element *sideElem) const = 0;

401
402
403
404
405
406
407
    /** \brief
     * Returns for a given element type number the element type number of the children.
     * For 1d and 2d this is always 0, because element type number are used in the 
     * 3d case only.
     */
    virtual int getChildType(int elType) const = 0;

408
    /** \brief
409
410
411
     * Traverses an edge/face of a given element (this includes also all children of 
     * the element having the same edge/face). All vertex dofs alonge this edge/face
     * are assembled and put together to a list. 
412
413
     *
     * \param[in]  feSpace         FE space which is used to get the dofs.
414
415
     * \param[in]  bound           Defines the edge/face of the element on which
     *                             all vertex dofs are assembled.
416
417
     * \param[out] dofs            List of dofs, where the result is stored.
     */
418
    virtual void getVertexDofs(FiniteElemSpace* feSpace, 
419
			       BoundaryObject bound,
420
			       DofContainer& dofs) const = 0;
421
422

    /** \brief
423
424
425
     * Traverses an edge/face of a given element (this includes also all children of 
     * the element having the same edge/face). All non vertex dofs alonge this edge/face 
     * are assembled and put together to a list.
426
427
     *
     * \param[in]  feSpace         FE space which is used to get the dofs.
428
429
     * \param[in]  bound           Defines the edge/face of the element on which
     *                             all non vertex dofs are assembled.
430
431
     * \param[out] dofs            All dofs are put to this dof list.
     */
432
    virtual void getNonVertexDofs(FiniteElemSpace* feSpace, 
433
				  BoundaryObject bound,
434
435
436
				  DofContainer& dofs) const = 0;


437
438
439
440
    /** \} */

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

441
    /// assignment operator
442
    Element& operator=(const Element& el);
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;

450
    /// Refines Element's leaf data
Thomas Witkowski's avatar
Thomas Witkowski committed
451
452
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) 
    {
453
454
455
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
456
	  ElementData *tmp = elementData->getDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
457
	  delete elementData;
458
459
460
	  elementData = tmp;
	}
      }
461
    }
462

463
    /// Coarsens Element's leaf data
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
    inline void coarsenElementData(Element* child1, Element* child2, int elType = 0) 
    {
466
467
      ElementData *childData;
      childData = child1->getElementData();
468
      if (childData) {
469
	childData->coarsenElementData(this, child1, child2, elType);
Thomas Witkowski's avatar
Thomas Witkowski committed
470
	delete childData;
471
472
473
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
474
      if (childData) {
475
	childData->coarsenElementData(this, child2, child1, elType);
Thomas Witkowski's avatar
Thomas Witkowski committed
476
	delete childData;
477
478
	child2->setElementData(NULL);
      }
479
    }
480

481
    /// Returns pointer to \ref elementData
Thomas Witkowski's avatar
Thomas Witkowski committed
482
483
    inline ElementData* getElementData() const 
    {
484
      return elementData;
485
    }
486

487
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
488
489
490
    inline ElementData* getElementData(int typeID) const 
    {
      if (elementData)
491
	return elementData->getElementData(typeID);
Thomas Witkowski's avatar
Thomas Witkowski committed
492

493
      return NULL;
494
    }
495

496
    /// Deletes the \ref elementData with a specific typeID.
497
    bool deleteElementData(int typeID);
498
499
500
501
502
503
504
505

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

508
    /// Returns whether Element's \ref newCoord is set
Thomas Witkowski's avatar
Thomas Witkowski committed
509
510
    inline bool isNewCoordSet() const 
    { 
511
      return (newCoord != NULL);
512
    }
513

514
    /// Frees memory for \ref newCoord
515
    void eraseNewCoord();
516
517
 
    /// Serialize the element to a file.
518
    void serialize(std::ostream &out);
519

520
    /// Deserialize an element from a file.
521
    void deserialize(std::istream &in);
522

523
524
525
    /// Sets Element's \ref dof pointer.
    void createNewDofPtrs();
 
526
  protected:
527
    /// Sets Element's \ref index. Used by friend class Mesh.
Thomas Witkowski's avatar
Thomas Witkowski committed
528
529
    inline void setIndex(int i) 
    {
530
      index = i;
531
    }
532

533
    /// Used by friend class Mesh while dofCompress
534
    void newDofFct1(const DOFAdmin*, std::vector<int>&);
535

536
    /// Used by friend class Mesh while dofCompress
537
538
539
    void newDofFct2(const DOFAdmin*);

    /// Changes old dofs to negative new dofs
540
541
    void changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
		     int n0, int nd0, int nd, int pos);
542
543
544

    /// Changes negative new dofs to positive
    void changeDofs2(int n0, int nd0, int nd, int pos);
545
546
547
548
549
550

  protected:
    /** \brief
     * Pointers to the two children of interior elements of the tree. Pointers
     * to NULL for leaf elements.
     */
551
    Element *child[2];
552
553
554
555
556

    /** \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 
557
558
559
     * 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 
560
561
     * barycenter, if present.
     */
562
    DegreeOfFreedom **dof;
563
564
565
566
567
568

    /** \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).
     */
569
    int index;
570
571
572
573
574
575

    /** \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.
     */
576
    int mark;
577
578
579
580
581
 
    /** \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 
582
583
     * information can be also produced by the traversal routines in the case of 
     * curved boundary.
584
585
586
     */
    WorldVector<double> *newCoord;

587
    /// Pointer to the Mesh this element belongs to
588
    Mesh* mesh;
589

590
    /// Pointer to Element's leaf data
591
    ElementData* elementData;
592

593
594
595
596
597
598
599
600
    /** \brief
     * Is true, if the DOF pointers are valid. In sequential computations this
     * is always the case. In parallel computations all macro elements are stored
     * in memory though not all of them are part of the mesh (because they are owned
     * by another rank). In this case, there are no valid DOFs on the element. 
     */
    bool dofValid;

601
602
603
604
605
606
    /** \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;
607
608
609

    friend class Mesh;
  };
610
611
612
613
614


  /// Writes the element hierarchie to a Graphviz dot-file. Using the dot-tool from
  /// Graphvis, this dot-file can be converted to a ps-file. Useful for debugging!
  void writeDotFile(Element *el, std::string filename, int maxLevels = -1);
615
616
617
618
}

#endif  // AMDIS_ELEMENT_H