Am Montag, 13. Mai 2022, finden Wartungsarbeiten am Gitlab-Server (Update auf neue Version statt). Der Dienst wird daher am Montag für einige Zeit nicht verfügbar sein.
On Monday, May 13th 2022, the Gitlab server will be updated. The service will therefore not be accessible for some time on Monday.

Element.h 18.3 KB
Newer Older
1
2
3
4
// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
5
// ==  http://www.amdis-fem.org                                              ==
6
7
// ==                                                                        ==
// ============================================================================
8
9
10
11
12
13
14
15
16
17
18
19
//
// Software License for AMDiS
//
// Copyright (c) 2010 Dresden University of Technology 
// All rights reserved.
// Authors: Simon Vey, Thomas Witkowski et al.
//
// This file is part of AMDiS
//
// See also license.opensource.txt in the distribution.


20
21
22
23
24
25

/** \file Element.h */

#ifndef AMDIS_ELEMENT_H
#define AMDIS_ELEMENT_H

26
#include "AMDiS_fwd.h"
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
#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:
69
    /// private standard constructor because an Element must know his Mesh
70
    Element() {}
71

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

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

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

82
    ///
83
84
    void deleteElementDOFs();

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

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

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

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

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

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

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

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

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

136
      return dof[i];
137
    }
138

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

    /** \} */

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

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

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

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

    /** \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
274
275
    inline void setNewCoord(WorldVector<double>* coord) 
    {
276
      newCoord = coord;
277
    }
278

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

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

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

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

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

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

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

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

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

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

    /** \} */


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

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

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

371
372
373
374
375
376
377
378
379
380
381
382
    /** \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;

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

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

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

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

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

402
403
404
405
406
407
408
    /** \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;

409
    /** \brief
410
411
412
     * 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. 
413
414
     *
     * \param[in]  feSpace         FE space which is used to get the dofs.
415
416
     * \param[in]  bound           Defines the edge/face of the element on which
     *                             all vertex dofs are assembled.
417
418
     * \param[out] dofs            List of dofs, where the result is stored.
     */
419
    virtual void getVertexDofs(FiniteElemSpace* feSpace, 
420
			       BoundaryObject bound,
421
			       DofContainer& dofs) const = 0;
422
423

    /** \brief
424
425
426
     * 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.
427
428
     *
     * \param[in]  feSpace         FE space which is used to get the dofs.
429
430
     * \param[in]  bound           Defines the edge/face of the element on which
     *                             all non vertex dofs are assembled.
431
432
     * \param[out] dofs            All dofs are put to this dof list.
     */
433
    virtual void getNonVertexDofs(FiniteElemSpace* feSpace, 
434
				  BoundaryObject bound,
435
436
437
				  DofContainer& dofs) const = 0;


438
439
440
441
    /** \} */

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

442
    /// assignment operator
443
    Element& operator=(const Element& el);
444
445
446
447
448
449
450

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

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

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

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

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

494
      return NULL;
495
    }
496

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

594
595
596
597
598
599
600
601
    /** \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;

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

    friend class Mesh;
  };
611
612
613
614
615


  /// 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);
616
617
618
619
}

#endif  // AMDIS_ELEMENT_H