Liebe Gitlab-Nutzerin, lieber Gitlab-Nutzer,
es ist nun möglich sich mittels des ZIH-Logins/LDAP an unserem Dienst anzumelden. Die Konten der externen Nutzer:innen sind über den Reiter "Standard" erreichbar.
Die Administratoren


Dear Gitlab user,
it is now possible to log in to our service using the ZIH login/LDAP. The accounts of external users can be accessed via the "Standard" tab.
The administrators

Element.h 18.7 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
      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
127

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
      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
135

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
      TEST_EXIT_DBG(dofValid)("DOFs are not valid in element %d!\n", index);
143

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 either vertex, edge or face of the element.
219 220
/*     template<typename T> */
/*     T getObject(int index);       */
221

222
    /// Returns the number of parts of type i in this element
223 224
    virtual int getGeo(GeoIndex i) const = 0;

225
    /// Returns Element's \ref mark
226
    inline const int getMark() const 
Thomas Witkowski's avatar
Thomas Witkowski committed
227
    { 
228
      return mark;
229
    }
230

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

238
    /// Returns Element's \ref index
Thomas Witkowski's avatar
Thomas Witkowski committed
239 240
    inline int getIndex() const 
    { 
241
      return index; 
242
    }
243

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

    /** \} */

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

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

262
    /// Sets \ref child[1]
Thomas Witkowski's avatar
Thomas Witkowski committed
263 264
    virtual void setSecondChild(Element *aChild) 
    {
265
      child[1] = aChild;
266
    }
267

268
    /// Sets \ref elementData of Element
Thomas Witkowski's avatar
Thomas Witkowski committed
269 270
    void setElementData(ElementData* ed) 
    {
271
      elementData = ed;
272
    }
273 274 275 276 277

    /** \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
278 279
    inline void setNewCoord(WorldVector<double>* coord) 
    {
280
      newCoord = coord;
281
    }
282

283
    /// Sets \ref mesh.
Thomas Witkowski's avatar
Thomas Witkowski committed
284 285
    inline void setMesh(Mesh *m) 
    {
286
      mesh = m;
287
    }
288

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

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

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

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

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

	dynamic_cast<LeafDataCoarsenableInterface*>(ld)->
	  setCoarseningErrorEstimate(row, est);
328
      } else {
329 330
	ERROR_EXIT("setEstimation only for leaf elements!\n");
      }
331
    }
332

333
    /// Sets Elements \ref mark = mark + 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
334 335
    inline void incrementMark() 
    {
336 337
      mark++;
    }
338

339
    /// Sets Elements \ref mark = mark - 1;
Thomas Witkowski's avatar
Thomas Witkowski committed
340 341
    inline void decrementMark() 
    {
342 343
      if (0 < mark) 
	mark--;
344
    }
345

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

    /** \} */


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

366
    /// Returns a copy of itself. Needed by Mesh to create Elements by a prototype. 
367 368 369
    virtual Element *clone() = 0;

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

375 376 377 378 379 380 381 382 383 384 385 386
    /** \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;

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

394
    /// Returns whether Element is a Line
395 396
    virtual bool isLine() const = 0;

397
    /// Returns whether Element is a Triangle
398 399
    virtual bool isTriangle() const = 0;

400
    /// Returns whether Element is a Tetrahedron
401 402
    virtual bool isTetrahedron() const = 0;

403
    /// Returns whether Element has sideElem as one of its sides.
404 405
    virtual bool hasSide(Element *sideElem) const = 0;

406 407 408 409 410 411 412
    /** \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;

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

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

443
    void getAllDofs(const FiniteElemSpace* feSpace, 
444 445
		    BoundaryObject bound, 
		    DofContainer& dofs);
446 447


448 449 450 451
    /** \} */

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

452
    /// assignment operator
453
    Element& operator=(const Element& el);
454 455 456 457 458 459 460

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

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

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

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

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

504
      return NULL;
505
    }
506

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

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

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

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

531
    /// Deserialize an element from a file.
532
    void deserialize(std::istream &in);
533

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

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

547
    /// Used by friend class Mesh while dofCompress
548 549 550
    void newDofFct2(const DOFAdmin*);

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

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

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

    /** \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 
568 569 570
     * 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 
571 572
     * barycenter, if present.
     */
573
    DegreeOfFreedom **dof;
574 575 576 577 578 579

    /** \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).
     */
580
    int index;
581 582 583 584 585 586

    /** \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.
     */
587
    int mark;
588 589 590 591 592
 
    /** \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 
593 594
     * information can be also produced by the traversal routines in the case of 
     * curved boundary.
595 596 597
     */
    WorldVector<double> *newCoord;

598
    /// Pointer to the Mesh this element belongs to
599
    Mesh* mesh;
600

601
    /// Pointer to Element's leaf data
602
    ElementData* elementData;
603

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

612 613 614 615 616 617
    /** \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;
618 619 620

    friend class Mesh;
  };
621 622 623 624 625


  /// 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);
626 627
}

628 629
#include "Element.hh"

630 631
#endif  // AMDIS_ELEMENT_H