Element.h 19.4 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
    /// Clone this Element and return a reference to it. Because also the DOFs
    /// are cloned, \ref Mesh::serializedDOfs must be used.
87
88
    Element* cloneWithDOFs();

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

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

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

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

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

119
    /// Returns \ref dof[i][j] which is the j-th DOF of the i-th node of Element.
120
    inline DegreeOfFreedom getDof(int i, int j) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
121
    { 
122
      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
123

124
      return dof[i][j];
125
    }
126

127
    /// Returns \ref dof[i] which is a pointer to the DOFs of the i-th node.
128
    inline const DegreeOfFreedom* getDof(int i) const 
Thomas Witkowski's avatar
Thomas Witkowski committed
129
    {
130
      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
131

132
      return dof[i];
133
    }
134

135
    /// Returns a pointer to the DOFs of this Element
136
    inline const DegreeOfFreedom** getDof() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
137
    {
138
      TEST_EXIT_DBG(dof != NULL)("DOFs are not valid in element %d!\n", index);
139

140
      return const_cast<const DegreeOfFreedom**>(dof);
141
    }
142

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

149
150
    /// Returns \ref elementData's error estimation, if Element is a leaf element
    /// and has leaf data. 
151
152
    inline double getEstimation(int row) const
    {
153
      if (isLeaf()) {
154
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
155
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
156
	TEST_EXIT_DBG(ld)("leaf data not estimatable!\n");
157

158
159
160
161
	return dynamic_cast<LeafDataEstimatableInterface*>(ld)->getErrorEstimate(row);
      }	
      
      return 0.0;
162
    }
163

164
165
    /// 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
166
167
    inline double getCoarseningEstimation(int row) 
    {
168
      if (isLeaf()) {
169
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
170
	ElementData *ld = elementData->getElementData(COARSENABLE);
171
	TEST_EXIT_DBG(ld)("element data not coarsenable!\n");
172

173
	return dynamic_cast<LeafDataCoarsenableInterface*>(ld)->getCoarseningErrorEstimate(row);
174
      }
175
176
      
      return 0.0;
177
    }
178

179
    /// Returns region of element if defined, -1 else.
180
181
    int getRegion() const;

182
    /// Returns local vertex number of the j-th vertex of the i-th edge
183
184
    virtual int getVertexOfEdge(int i, int j) const = 0; 

185
186
    /// Returns local vertex number of the vertexIndex-th vertex of the
    /// positionIndex-th part of type position (vertex, edge, face)
187
    virtual int getVertexOfPosition(GeoIndex position,
188
189
				    int positionIndex,
				    int vertexIndex) const = 0;
190

191
    ///
192
193
    virtual int getPositionOfVertex(int side, int vertex) const = 0;

194
    ///
195
196
    virtual int getEdgeOfFace(int face, int edge) const = 0;

197
    ///
198
199
200
201
    virtual DofEdge getEdge(int localEdgeIndex) const = 0;

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

203
    /// Returns the number of parts of type i in this element
204
205
    virtual int getGeo(GeoIndex i) const = 0;

206
    /// Returns Element's \ref mark
207
    inline int getMark() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
208
    { 
209
      return mark;
210
    }
211

212
    /// Returns \ref newCoord[i]
Thomas Witkowski's avatar
Thomas Witkowski committed
213
214
    inline double getNewCoord(int i) const 
    {
215
216
217
	TEST_EXIT_DBG(newCoord)("newCoord = NULL\n");
	return (*newCoord)[i];
    }
218

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

225
    /// Returns \ref newCoord
Thomas Witkowski's avatar
Thomas Witkowski committed
226
227
    inline WorldVector<double>* getNewCoord() const 
    { 
228
      return newCoord; 
229
    }
230
231
232
233
234
235
236

    /** \} */

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

237
    /// Sets \ref child[0]
Thomas Witkowski's avatar
Thomas Witkowski committed
238
239
    virtual void setFirstChild(Element *aChild) 
    {
240
      child[0] = aChild;
241
    }
242

243
    /// Sets \ref child[1]
Thomas Witkowski's avatar
Thomas Witkowski committed
244
245
    virtual void setSecondChild(Element *aChild) 
    {
246
      child[1] = aChild;
247
    }
248

249
    /// Sets \ref elementData of Element
Thomas Witkowski's avatar
Thomas Witkowski committed
250
251
    void setElementData(ElementData* ed) 
    {
252
      elementData = ed;
253
    }
254

255
256
    /// 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
257
258
    inline void setNewCoord(WorldVector<double>* coord) 
    {
259
      newCoord = coord;
260
    }
261

262
    /// Sets \ref mesh.
Thomas Witkowski's avatar
Thomas Witkowski committed
263
264
    inline void setMesh(Mesh *m) 
    {
265
      mesh = m;
266
    }
267

268
    /// Sets the pointer to the DOFs of the i-th node of Element
269
    inline void setDof(int pos, DegreeOfFreedom* p) 
Thomas Witkowski's avatar
Thomas Witkowski committed
270
    {
271
      dof[pos] = p;
272
    }
273

274
275
276
277
278
279
280
281
282
    /// Delets the main DOF pointer. The caller must ensure that the DOFs are
    /// either freed correctly elsewhere or that they are still used by other
    /// elements.
    void delDofPtr()
    {
      delete [] dof;
      dof = NULL;
    }

283
284
    /// 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.
285
286
    inline void setEstimation(double est, int row)
    {
287
288
      FUNCNAME("Element::setEstimation()");

289
      if (isLeaf()) {
290
	TEST_EXIT_DBG(elementData)("Leaf element %d without leaf data!\n", index);
291
	ElementData *ld = elementData->getElementData(ESTIMATABLE);
292
	TEST_EXIT_DBG(ld)("Leaf data %d not estimatable!\n", index);
293
294
295

	dynamic_cast<LeafDataEstimatableInterface*>(ld)->
	  setErrorEstimate(row, est);
296
      } else {
297
298
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
299
    }
300

301
302
    /// 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.
303
304
    inline void setCoarseningEstimation(double est, int row)
    {
305
      if (isLeaf()) {
306
	TEST_EXIT_DBG(elementData)("leaf element without leaf data\n");
307
	ElementData *ld = elementData->getElementData(COARSENABLE);
308
	TEST_EXIT_DBG(ld)("leaf data not coarsenable\n");
309
310
311

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
312
      } else {
313
314
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
315
    }
316

317
    /// Sets Elements \ref mark = mark + 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
318
319
    inline void incrementMark() 
    {
320
321
      mark++;
    }
322

323
    /// Sets Elements \ref mark = mark - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
324
325
    inline void decrementMark() 
    {
326
327
      if (0 < mark) 
	mark--;
328
    }
329

330
    /// Sets Element's \ref mark
331
    inline void setMark(int m) 
Thomas Witkowski's avatar
Thomas Witkowski committed
332
    {
333
      mark = m;
334
    }
335
336
337
338
339
340
341
342

    /** \} */


    /** \name pure virtual methods 
     * \{ 
     */

343
344
345
    /// Orient the vertices of edges/faces. Used
    /// by Estimator for the jumps => same quadrature nodes from both sides!
    virtual void sortFaceIndices(int face, FixVec<int, WORLD> &vec) const = 0;
346

347
348
    /// Returns a copy of itself. Needed by Mesh to create Elements by 
    /// a prototype. 
349
350
    virtual Element *clone() = 0;

351
352
353
    /// 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.
354
355
    virtual int getSideOfChild(int childnr, int sidenr, int elType = 0) const = 0;

356
    /** \brief
357
358
359
     * 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.
360
361
362
363
364
365
366
367
368
     *
     * \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;

369
370
371
372
373
374
    /// 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;
375

376
    /// Returns whether Element is a Line
377
378
    virtual bool isLine() const = 0;

379
    /// Returns whether Element is a Triangle
380
381
    virtual bool isTriangle() const = 0;

382
    /// Returns whether Element is a Tetrahedron
383
384
    virtual bool isTetrahedron() const = 0;

385
    /// Returns whether Element has sideElem as one of its sides.
386
387
    virtual bool hasSide(Element *sideElem) const = 0;

388
389
390
391
392
393
394
    /** \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;

395
    /** \brief
396
397
398
399
     * Traverses a vertex/edge/face of a given element (this includes also all
     * children of the element having the same edge/face). All DOFs on mesh
     * nodes alonge this vertex/edge/face are assembled and put together to 
     * a list.
400
     *
401
402
403
404
405
406
407
408
     * \param[in]  feSpace     FE space which is used to get the dofs.
     * \param[in]  bound       Defines the vertex/edge/face of the element on
     *                         which all vertex dofs are assembled.
     * \param[out] dofs        List of dofs, where the result is stored.
     * \param[in]  baseDofPtr  If true, the base DOF pointes are stored. Thus,
     *                         dof* [\ref dof] of the element is inserted. If 
     *                         false, &(dof[.][n0]) is put to the result vector, 
     *                         with n0 beging the number of predofs.
409
     */
410
    virtual void getNodeDofs(const FiniteElemSpace* feSpace, 
411
			     BoundaryObject bound,
412
413
			     DofContainer& dofs,
			     bool baseDofPtr = false) const = 0;
414
415

    /** \brief
416
417
418
419
     * Traverses a vertex/edge/face of a given element (this includes also all
     * children of the element having the same edge/face). All DOFs belonging
     * to higher order basis functions alonge this vertex/edge/face are 
     * assembled and put together to a list.
420
     *
421
422
423
424
425
426
427
428
     * \param[in]  feSpace     FE space which is used to get the dofs.
     * \param[in]  bound       Defines the edge/face of the element on which
     *                         all non vertex dofs are assembled.
     * \param[out] dofs        All dofs are put to this dof list.
     * \param[in]  baseDofPtr  If true, the base DOF pointes are stored. Thus,
     *                         dof* [\ref dof] of the element is inserted. If 
     *                         false, &(dof[.][n0]) is put to the result vector, 
     *                         with n0 beging the number of predofs.
429
430
431
432
     * \param[out] dofGeoIndex Optional, the function can store to each DOF in
     *                         the DofContainer dofs the geometric index, thus
     *                         identifing the DOF to be a vertex, edge, face or
     *                         center DOF.
433
     */
434
435
436
437
438
439
440
    virtual void 
    getHigherOrderDofs(const FiniteElemSpace* feSpace, 
		       BoundaryObject bound,
		       DofContainer& dofs,
		       bool baseDofPtr = false,
		       vector<GeoIndex>* dofGeoIndex = NULL) const = 0;

441
442
    /// Combines \ref getNodeDofs and \ref getHigherOrderDofs to one function. 
    /// See parameter description there.
443
    void getAllDofs(const FiniteElemSpace* feSpace, 
444
		    BoundaryObject bound, 
445
		    DofContainer& dofs,
446
447
		    bool baseDofPtr = false,
		    vector<GeoIndex>* dofGeoIndex = NULL);
448

Thomas Witkowski's avatar
Thomas Witkowski committed
449
450
    virtual void getSubBoundary(BoundaryObject bound, 
				vector<BoundaryObject> &subBound) const = 0;
451

452
453
454
455
    /** \} */

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

456
    /// assignment operator
457
    Element& operator=(const Element& el);
458

Thomas Witkowski's avatar
Thomas Witkowski committed
459
460
    /// 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
461
462
    int oppVertex(FixVec<DegreeOfFreedom*, DIMEN> pdof) const;

463
    /// Refines Element's leaf data
Thomas Witkowski's avatar
Thomas Witkowski committed
464
465
    inline void refineElementData(Element* child1, Element* child2, int elType = 0) 
    {
466
467
468
      if (elementData) {
	bool remove = elementData->refineElementData(this, child1, child2, elType);
	if (remove) {
469
	  ElementData *tmp = elementData->getDecorated();
Thomas Witkowski's avatar
Thomas Witkowski committed
470
	  delete elementData;
471
472
473
	  elementData = tmp;
	}
      }
474
    }
475

476
    /// Coarsens Element's leaf data
477
478
    inline void coarsenElementData(Element* child1, Element* child2, 
				   int elType = 0) 
Thomas Witkowski's avatar
Thomas Witkowski committed
479
    {
480
481
      ElementData *childData;
      childData = child1->getElementData();
482
      if (childData) {
483
	childData->coarsenElementData(this, child1, child2, elType);
Thomas Witkowski's avatar
Thomas Witkowski committed
484
	delete childData;
485
486
487
	child1->setElementData(NULL);
      }
      childData = child2->getElementData();
488
      if (childData) {
489
	childData->coarsenElementData(this, child2, child1, elType);
Thomas Witkowski's avatar
Thomas Witkowski committed
490
	delete childData;
491
492
	child2->setElementData(NULL);
      }
493
    }
494

495
    /// Returns pointer to \ref elementData
Thomas Witkowski's avatar
Thomas Witkowski committed
496
497
    inline ElementData* getElementData() const 
    {
498
      return elementData;
499
    }
500

501
    ///
Thomas Witkowski's avatar
Thomas Witkowski committed
502
503
504
    inline ElementData* getElementData(int typeID) const 
    {
      if (elementData)
505
	return elementData->getElementData(typeID);
Thomas Witkowski's avatar
Thomas Witkowski committed
506

507
      return NULL;
508
    }
509

510
    /// Deletes the \ref elementData with a specific typeID.
511
    bool deleteElementData(int typeID);
512
513
514
515
516
517
518
519

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

522
    /// Returns whether Element's \ref newCoord is set
Thomas Witkowski's avatar
Thomas Witkowski committed
523
524
    inline bool isNewCoordSet() const 
    { 
525
      return (newCoord != NULL);
526
    }
527

528
    /// Frees memory for \ref newCoord
529
    void eraseNewCoord();
530
531
 
    /// Serialize the element to a file.
532
    void serialize(std::ostream &out);
533

534
    /// Deserialize an element from a file.
535
    void deserialize(std::istream &in);
536

537
    /// Sets Element's \ref dof pointer.
538
    void createNewDofPtrs(bool setDofs = false);
539
 
540
  protected:
541
    /// Sets Element's \ref index. Used by friend class Mesh.
Thomas Witkowski's avatar
Thomas Witkowski committed
542
543
    inline void setIndex(int i) 
    {
544
      index = i;
545
    }
546

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

550
    /// Used by friend class Mesh while dofCompress
551
552
553
    void newDofFct2(const DOFAdmin*);

    /// Changes old dofs to negative new dofs
554
555
    void changeDofs1(const DOFAdmin* admin, std::vector<int>& newDofIndex,
		     int n0, int nd0, int nd, int pos);
556
557
558

    /// Changes negative new dofs to positive
    void changeDofs2(int n0, int nd0, int nd, int pos);
559
560

  protected:
561
562
    /// Pointers to the two children of interior elements of the tree. Pointers
    /// to NULL for leaf elements.
563
    Element *child[2];
564

Thomas Witkowski's avatar
Thomas Witkowski committed
565
566
567
568
569
570
571
    /// 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.
572
    DegreeOfFreedom **dof;
573

Thomas Witkowski's avatar
Thomas Witkowski committed
574
575
576
    /// 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).
577
    int index;
578

Thomas Witkowski's avatar
Thomas Witkowski committed
579
580
581
    /// 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.
582
    int mark;
583
 
Thomas Witkowski's avatar
Thomas Witkowski committed
584
585
586
587
588
    /// 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.
589
590
    WorldVector<double> *newCoord;

591
    /// Pointer to the Mesh this element belongs to
592
    Mesh* mesh;
593

594
    /// Pointer to Element's leaf data
595
    ElementData* elementData;
596

597
598
599
    /// 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.
600
    static std::map<DegreeOfFreedom*, bool> deletedDOFs;
601
602
603

    friend class Mesh;
  };
604
605
606
607
608


  /// 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);
609
610
}

611
612
#include "Element.hh"

613
614
#endif  // AMDIS_ELEMENT_H