Mesh.h 25.9 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
26
27
28
29
30
31
32
33
34
35
36

/** \file Mesh.h */

/** \defgroup Triangulation Triangulation module
 * @{ <img src="triangulation.png"> @}
 *
 * Example:
 *
 * @{ <img src="hierarchicalMesh.png"> @}
 *
 * \brief
 * Contains all triangulation classes.
 */

#ifndef AMDIS_MESH_H
#define AMDIS_MESH_H

Thomas Witkowski's avatar
Thomas Witkowski committed
37
38
39
40
#include <deque>
#include <set>
#include <stdio.h>
#include "AMDiS_fwd.h"
41
42
43
44
45
46
47
48
49
50
51
52
#include "DOFAdmin.h"
#include "Line.h"
#include "Triangle.h"
#include "Tetrahedron.h"
#include "Element.h"
#include "ElInfo.h"
#include "FixVec.h"
#include "Serializable.h"
#include "BoundaryCondition.h"

namespace AMDiS {

53
54
55
  using namespace std;


56
57
58
59
60
61
62
  /** \ingroup Triangulation 
   * \brief
   * A Mesh holds all information about a triangulation. 
   */
  class Mesh : public Serializable
  {
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
63
    /// Creates a mesh with the given name of dimension dim
64
    Mesh(string name, int dim);
65

Thomas Witkowski's avatar
Thomas Witkowski committed
66
    /// Destructor
Thomas Witkowski's avatar
Thomas Witkowski committed
67
    ~Mesh();
68

Thomas Witkowski's avatar
Thomas Witkowski committed
69
    /// Reads macro triangulation.
70
71
    void initialize();

Thomas Witkowski's avatar
Thomas Witkowski committed
72
    /// Assignment operator
73
74
75
76
77
78
    Mesh& operator=(const Mesh&);

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

Thomas Witkowski's avatar
Thomas Witkowski committed
79
80
    /// Returns geometric information about this mesh. With GeoIndex p it is 
    /// specified which information is requested.
81
82
    inline int getGeo(GeoIndex p) const 
    { 
83
      return Global::getGeo(p, dim); 
84
    }
85

Thomas Witkowski's avatar
Thomas Witkowski committed
86
    /// Returns \ref name of the mesh
87
    inline string getName() const 
88
    { 
89
      return name; 
90
    }
91

Thomas Witkowski's avatar
Thomas Witkowski committed
92
    /// Returns \ref dim of the mesh
93
    inline int getDim() const
94
    { 
95
      return dim; 
96
    }
97

98
99
    /// Returns \ref nDofEl of the mesh
    inline const int getNumberOfAllDofs() const 
100
    { 
101
      return nDofEl; 
102
    }
103

Thomas Witkowski's avatar
Thomas Witkowski committed
104
    /// Returns \ref nNodeEl of the mesh
105
106
    inline const int getNumberOfNodes() const 
    { 
107
      return nNodeEl; 
108
    }
109

Thomas Witkowski's avatar
Thomas Witkowski committed
110
    /// Returns \ref nVertices of the mesh
111
112
    inline const int getNumberOfVertices() const 
    { 
113
      return nVertices; 
114
    }
115

Thomas Witkowski's avatar
Thomas Witkowski committed
116
    /// Returns \ref nEdges of the mesh 
117
118
    inline const int getNumberOfEdges() const 
    { 
119
      return nEdges; 
120
    }
121

Thomas Witkowski's avatar
Thomas Witkowski committed
122
    /// Returns \ref nFaces of the mesh 
123
124
    inline const int getNumberOfFaces() const 
    { 
125
      return nFaces; 
126
    }
127

Thomas Witkowski's avatar
Thomas Witkowski committed
128
    /// Returns \ref nLeaves of the mesh 
129
130
    inline const int getNumberOfLeaves() const 
    { 
131
      return nLeaves; 
132
    }
133

Thomas Witkowski's avatar
Thomas Witkowski committed
134
    /// Returns \ref nElements of the mesh
135
136
    inline const int getNumberOfElements() const 
    { 
137
      return nElements; 
138
    }
139

Thomas Witkowski's avatar
Thomas Witkowski committed
140
    /// Returns \ref maxEdgeNeigh of the mesh
141
142
    inline const int getMaxEdgeNeigh() const 
    { 
143
      return maxEdgeNeigh; 
144
    }
145

Thomas Witkowski's avatar
Thomas Witkowski committed
146
    /// Returns \ref parametric of the mesh
147
148
    inline Parametric *getParametric() const 
    { 
149
      return parametric; 
150
    }
151

Thomas Witkowski's avatar
Thomas Witkowski committed
152
    /// Returns \ref diam of the mesh
153
154
    inline const WorldVector<double>& getDiameter() const 
    { 
155
      return diam; 
156
    }
157

158
159
    /// Returns nDof[i] of the mesh
    inline const int getNumberOfDofs(int i) const 
160
    { 
161
      TEST_EXIT_DBG(i <= dim)("Wrong index: %d %d\n", i, dim);
162
      return nDof[i]; 
163
    }
164

Thomas Witkowski's avatar
Thomas Witkowski committed
165
    /// Returns \ref elementPrototype of the mesh
166
167
    inline Element* getElementPrototype() 
    { 
168
      return elementPrototype; 
169
    }
170

Thomas Witkowski's avatar
Thomas Witkowski committed
171
    /// Returns \ref leafDataPrototype of the mesh
172
173
    inline ElementData* getElementDataPrototype() 
    { 
174
      return elementDataPrototype; 
175
    }
176

Thomas Witkowski's avatar
Thomas Witkowski committed
177
    /// Returns node[i] of the mesh 
178
179
    inline int getNode(int i) const 
    { 
180
      return node[i]; 
181
    }
182

183
184
185
186
    /// Allocates the number of DOFs needed at position and registers the DOFs
    /// at the DOFAdmins. The number of needed DOFs is the sum over the needed
    /// DOFs of all DOFAdmin objects belonging to this mesh. 
    /// The return value is a pointer to the first allocated DOF. 
187
    DegreeOfFreedom *getDof(GeoIndex position);
188

Thomas Witkowski's avatar
Thomas Witkowski committed
189
    /// Returns *(\ref admin[i]) of the mesh
190
    inline const DOFAdmin& getDofAdmin(int i) const 
191
    {
192
      return *(admin[i]);
193
    }
194

195
196
197
    /// Creates a DOFAdmin with name lname. nDof specifies how many DOFs 
    /// are needed at the different positions (see \ref DOFAdmin::nrDOF).
    /// A pointer to the created DOFAdmin is returned.
198
    const DOFAdmin* createDOFAdmin(string lname, DimVec<int> nDof);
199

200
201
    /// Returns the size of \ref admin which is the number of the DOFAdmins
    /// belonging to this mesh
202
203
    const int getNumberOfDOFAdmin() const 
    {
204
      return admin.size();
205
    }
206

207
208
    /// Returns the size of \ref macroElements which is the number of
    /// of macro elements of this mesh
209
210
    const int getNumberOfMacros() const 
    {
211
      return macroElements.size();
212
    }
213

Thomas Witkowski's avatar
Thomas Witkowski committed
214
    /// Returns a DOFAdmin which at least manages vertex DOFs
215
216
    const DOFAdmin* getVertexAdmin() const;

Thomas Witkowski's avatar
Thomas Witkowski committed
217
218
    /// Allocates an array of DOF pointers. The array holds one pointer for 
    /// each node.
Thomas Witkowski's avatar
Thomas Witkowski committed
219
    DegreeOfFreedom **createDofPtrs();
220

Thomas Witkowski's avatar
Thomas Witkowski committed
221
    /// Returns \ref preserveCoarseDOFs of the mesh
222
223
    inline bool queryCoarseDOFs() const 
    { 
224
      return preserveCoarseDOFs;
225
    }
226

Thomas Witkowski's avatar
Thomas Witkowski committed
227
    /// Returns an iterator to the begin of \ref macroElements
228
    inline deque<MacroElement*>::iterator firstMacroElement() 
229
    {
230
      return macroElements.begin();
231
    }
232

Thomas Witkowski's avatar
Thomas Witkowski committed
233
    /// Returns macroElements[i].
234
235
    inline MacroElement *getMacroElement(int i) 
    { 
236
      return macroElements[i]; 
237
    }
238

Thomas Witkowski's avatar
Thomas Witkowski committed
239
    /// Returns an iterator to the end of \ref macroElements
240
    inline deque<MacroElement*>::iterator endOfMacroElements() 
241
    {
242
      return macroElements.end();
243
    }
244

245
    /// Returns \ref macroElements, the list of all macro elements in the mesh.
246
    deque<MacroElement*>& getMacroElements()
247
248
249
250
    {
      return macroElements;
    }

251
252
253
254
255
256
    /** \} */

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

Thomas Witkowski's avatar
Thomas Witkowski committed
257
    /// Sets \ref name of the mesh
258
    inline void setName(string aName) 
259
    { 
260
      name = aName;
261
    }
262

Thomas Witkowski's avatar
Thomas Witkowski committed
263
    /// Sets \ref nVertices of the mesh
264
265
    inline void setNumberOfVertices(int n) 
    { 
266
      nVertices = n; 
267
    }
268

Thomas Witkowski's avatar
Thomas Witkowski committed
269
    /// Sets \ref nFaces of the mesh
270
271
    inline void setNumberOfFaces(int n) 
    { 
272
      nFaces = n; 
273
    }
274

Thomas Witkowski's avatar
Thomas Witkowski committed
275
    /// Increments \ref nVertices by inc
276
277
    inline void incrementNumberOfVertices(int inc) 
    { 
278
      nVertices += inc; 
279
    }
280
 
Thomas Witkowski's avatar
Thomas Witkowski committed
281
    /// Sets \ref nEdges of the mesh
282
283
    inline void setNumberOfEdges(int n) 
    { 
284
      nEdges = n; 
285
    }
286

Thomas Witkowski's avatar
Thomas Witkowski committed
287
    /// Increments \ref nEdges by inc
288
289
    inline void incrementNumberOfEdges(int inc) 
    { 
290
      nEdges += inc; 
291
    }
292

Thomas Witkowski's avatar
Thomas Witkowski committed
293
    /// Increments \ref nFaces by inc
294
295
    inline void incrementNumberOfFaces(int inc) 
    { 
296
      nFaces += inc; 
297
    }
298

Thomas Witkowski's avatar
Thomas Witkowski committed
299
    /// Sets \ref nLeaves of the mesh
300
301
    inline void setNumberOfLeaves(int n) 
    { 
302
      nLeaves = n; 
303
    }
304

Thomas Witkowski's avatar
Thomas Witkowski committed
305
    /// Increments \ref nLeaves by inc
306
307
    inline void incrementNumberOfLeaves(int inc) 
    { 
308
      nLeaves += inc; 
309
    }
310

Thomas Witkowski's avatar
Thomas Witkowski committed
311
    /// Sets \ref nElements of the mesh
312
313
    inline void setNumberOfElements(int n) 
    { 
314
      nElements = n; 
315
    }
316

Thomas Witkowski's avatar
Thomas Witkowski committed
317
    /// Increments \ref nElements by inc
318
319
    inline void incrementNumberOfElements(int inc) 
    { 
320
      nElements += inc; 
321
    }
322

Thomas Witkowski's avatar
Thomas Witkowski committed
323
    /// Sets *\ref diam to w
324
325
    void setDiameter(const WorldVector<double>& w);

Thomas Witkowski's avatar
Thomas Witkowski committed
326
    /// Sets (*\ref diam)[i] to d
327
328
    void setDiameter(int i, double d);

Thomas Witkowski's avatar
Thomas Witkowski committed
329
    /// Sets \ref preserveCoarseDOFs = true
330
331
    inline void retainCoarseDOFs() 
    {
332
      preserveCoarseDOFs = true;
333
    }
334

Thomas Witkowski's avatar
Thomas Witkowski committed
335
    /// Sets \ref preserveCoarseDOFs = b
336
337
    inline void setPreserveCoarseDOFs(bool b) 
    {
338
      preserveCoarseDOFs = b;
339
    }
340

Thomas Witkowski's avatar
Thomas Witkowski committed
341
    /// Sets \ref preserveCoarseDOFs = false
342
343
    inline void noCoarseDOFs() 
    {
344
      preserveCoarseDOFs = false;
345
    }
346

Thomas Witkowski's avatar
Thomas Witkowski committed
347
    /// Sets \ref elementPrototype of the mesh
348
349
    inline void setElementPrototype(Element* prototype) 
    {
350
      elementPrototype = prototype;
351
352
    }
    
Thomas Witkowski's avatar
Thomas Witkowski committed
353
    /// Sets \ref elementDataPrototype of the mesh
354
355
    inline void setElementDataPrototype(ElementData* prototype) 
    {
356
      elementDataPrototype = prototype;
357
    }
358

Thomas Witkowski's avatar
Thomas Witkowski committed
359
    ///
360
361
    inline void setParametric(Parametric *param) 
    {
362
      parametric = param;
363
    }
364

Thomas Witkowski's avatar
Thomas Witkowski committed
365
    ///
366
367
    inline void setMaxEdgeNeigh(int m) 
    { 
368
      maxEdgeNeigh = m; 
369
    }
370
371
372
  
    /** \} */

Thomas Witkowski's avatar
Thomas Witkowski committed
373
    /// Creates a new Element by cloning \ref elementPrototype
374
375
    Element* createNewElement(Element *parent = NULL);

Thomas Witkowski's avatar
Thomas Witkowski committed
376
    /// Creates a new ElInfo dependent of \ref dim of the mesh
377
378
    ElInfo* createNewElInfo();

Thomas Witkowski's avatar
Thomas Witkowski committed
379
    /// Frees DOFs at the given position pointed by dof 
380
    void freeDof(DegreeOfFreedom* dof, GeoIndex position);
381

Thomas Witkowski's avatar
Thomas Witkowski committed
382
    /// Frees memory for the given element el
383
384
    void freeElement(Element* el);

Thomas Witkowski's avatar
Thomas Witkowski committed
385
    /// Performs DOF compression for all DOFAdmins (see \ref DOFAdmin::compress)
386
387
    void dofCompress();

Thomas Witkowski's avatar
Thomas Witkowski committed
388
    /// Adds a DOFAdmin to the mesh
Thomas Witkowski's avatar
Thomas Witkowski committed
389
    void addDOFAdmin(DOFAdmin *admin);
390

391
392
393
    /// Recalculates the number of leave elements.
    void updateNumberOfLeaves();

Thomas Witkowski's avatar
Thomas Witkowski committed
394
    /// Clears \ref macroElements
395
396
    inline void clearMacroElements() 
    { 
397
      macroElements.clear();
398
    }
399
  
Thomas Witkowski's avatar
Thomas Witkowski committed
400
    /// Adds a macro element to the mesh
401
402
    void addMacroElement(MacroElement* me);

403
404
405
    /// Removes a set of macro elements from the mesh. This works only for the 
    /// case, that there are no global or local refinements, i.e., all macro 
    /// elements have no children.
406
    void removeMacroElements(std::set<MacroElement*>& macros,
407
			     vector<const FiniteElemSpace*>& feSpaces);
408

Thomas Witkowski's avatar
Thomas Witkowski committed
409
    /// Frees the array of DOF pointers (see \ref createDofPtrs)
410
    void freeDofPtrs(DegreeOfFreedom **ptrs);
411

Thomas Witkowski's avatar
Thomas Witkowski committed
412
    /// Used by \ref findElementAtPoint. 
413
414
415
    bool findElInfoAtPoint(const WorldVector<double>& xy,
			   ElInfo *el_info,
			   DimVec<double>& bary,
416
			   const MacroElement *start_mel,
417
418
			   const WorldVector<double> *xy0,
			   double *sp);
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454

    /** \brief
     * Access to an element at world coordinates xy. Some applications need the 
     * access to elements at a special location in world coordinates. Examples 
     * are characteristic methods for convection problems, or the implementation
     * of a special right hand side like point evaluations or curve integrals.
     * For such purposes, a routine is available which returns an element pointer
     * and corresponding barycentric coordinates.
     *
     * \param xy world coordinates of point
     * \param elp return address for a pointer to the element at xy
     * \param pary returns barycentric coordinates of xy
     * \param start_mel initial guess for the macro element containing xy or NULL
     * \param xy0 start point from a characteristic method, see below, or NULL
     * \param sp return address for relative distance to domain boundary in a 
     *        characteristic method, see below, or NULL
     * \return true is xy is inside the domain , false otherwise
     * 
     * For a characteristic method, where \f$ xy = xy_0 - V\tau \f$, it may be 
     * convenient to know the point on the domain's boundary which lies on the 
     * line segment between the old point xy0 and the new point xy, in case that 
     * xy is outside the domain. Such information is returned when xy0 and a 
     * pointer sp!=NULL are supplied: *sp is set to the value s such that 
     * \f$ xy_0 +s (xy -xy_0) \in \partial Domain \f$, and the element and local 
     * coordinates corresponding to that boundary point will be returned via elp 
     * and bary.
     *
     * The implementation of findElementAtPoint() is based on the transformation 
     * from world to local coordinates, available via the routine worldToCoord(),
     * At the moment, findElementAtPoint() works correctly only for domains with 
     * non-curved boundary. This is due to the fact that the implementation first
     * looks for the macro-element containing xy and then finds its path through 
     * the corresponding element tree based on the macro barycentric coordinates.
     * For non-convex domains, it is possible that in some cases a point inside
     * the domain is considered as external.
     */
455
456
457
    bool findElementAtPoint(const WorldVector<double>& xy,
			    Element **elp, 
			    DimVec<double>& bary,
458
			    const MacroElement *start_mel,
459
460
			    const WorldVector<double> *xy0,
			    double *sp);
461

462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
    /** \brief
     * Returns for a given dof its world coordinates in this mesh. Because we do
     * not have any direct connection between dofs and coordinates, this function
     * has to search for the element in this mesh, that contains the dof. Than the
     * coordinates can be computed. Therefore, this function is very costly and
     * should be used for debugging purpose only.
     *
     * @param[in]    dof       A pointer to the dof we have to search for.
     * @param[in]    feSpace   The fe space to be used for the search.
     * @param[out]   coords    World vector that stores the coordinates of the dof.
     *
     * The function returns true, if the dof was found, otherwise false.
     */
    bool getDofIndexCoords(const DegreeOfFreedom* dof, 
			   const FiniteElemSpace* feSpace,
477
478
479
480
			   WorldVector<double>& coords)
    {
      return getDofIndexCoords(*dof, feSpace, coords);
    }
481

482

483
484
    /// This function is equal to \ref getDofIndexCoords as defined above, but
    /// takes a DOF index instead of a DOF pointer.
485
486
487
    bool getDofIndexCoords(DegreeOfFreedom dof, 
			   const FiniteElemSpace* feSpace,
			   WorldVector<double>& coords);
488

489
    /** \brief
490
491
492
     * Traverse the whole mesh and stores to each DOF of the given finite
     * element space the coordinates in a given DOFVector. Works in the same
     * way as the function \ref getDofIndexCoords defined above.
493
     *
494
495
     * @param[in]   feSpace   The FE space to be used for the search.
     * @param[out]  coords    DOF vector that stores the coordinates to each DOF.
496
497
498
499
     */
    void getDofIndexCoords(const FiniteElemSpace* feSpace,
			   DOFVector<WorldVector<double> >& coords);

500
501
502
503
504
505
    /** \brief
     * Traverse the mesh and get all DOFs in this mesh for a given FE space.
     *
     * @param[in]   feSpace   The FE space to be used for collecting DOFs.
     * @param[out]  allDofs   The set which is filled with all DOFs.
     */
506
    void getAllDofs(const FiniteElemSpace *feSpace, 
507
508
		    std::set<const DegreeOfFreedom*>& allDofs);

Thomas Witkowski's avatar
Thomas Witkowski committed
509
    /// Returns FILL_ANY_?D
510
511
    inline static const Flag& getFillAnyFlag(int dim) 
    {
512
      switch (dim) {
513
514
515
516
517
518
519
520
521
522
523
524
525
      case 1:
	return FILL_ANY_1D;
	break;
      case 2:
	return FILL_ANY_2D;
	break;
      case 3:
	return FILL_ANY_3D;
	break;
      default:
	ERROR_EXIT("invalid dim\n");
	return FILL_ANY_1D;
      }
526
    }
527

Thomas Witkowski's avatar
Thomas Witkowski committed
528
    /// Serialize the mesh to a file.
529
    void serialize(ostream &out);
530

Thomas Witkowski's avatar
Thomas Witkowski committed
531
    /// Deserialize a mesh from a file.
532
    void deserialize(istream &in);
533

Thomas Witkowski's avatar
Thomas Witkowski committed
534
    /// Returns \ref elementIndex and increments it by 1.
535
536
    inline int getNextElementIndex() 
    { 
537
      return elementIndex++; 
538
    }
539

Thomas Witkowski's avatar
Thomas Witkowski committed
540
    /// Returns \ref initialized.
541
542
    inline bool isInitialized() 
    {
543
544
      return initialized; 
    }
545
  
Thomas Witkowski's avatar
Thomas Witkowski committed
546
    ///
547
    inline map<BoundaryType, VertexVector*>& getPeriodicAssociations() 
548
    {
549
      return periodicAssociations;
550
    }
551

552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
    /// Returns the periodic association for a specific boundary type.
    inline VertexVector& getPeriodicAssociations(BoundaryType b)
    {
      FUNCNAME("Mesh::getPeriodicAssociations()");

      TEST_EXIT_DBG(periodicAssociations.count(b) == 1)
	("There are no periodic assoications for boundary type %d!\n", b);

      return (*(periodicAssociations[b]));
    }

    
    /// Returns whether the given boundary type is periodic, i.e., if there is
    /// a periodic association for this boundary type.
    inline bool isPeriodicAssociation(BoundaryType b)
    {
      return (periodicAssociations.count(b) == 1 ? true : false);
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
571
    ///
572
573
    bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

Thomas Witkowski's avatar
Thomas Witkowski committed
574
    ///
575
576
    bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

Thomas Witkowski's avatar
Thomas Witkowski committed
577
    /// Returns \macroFileInfo
578
579
    inline MacroInfo* getMacroFileInfo() 
    { 
580
      return macroFileInfo;
581
    }
582

583
584
585
586
587
588
589
590
591
592
593
594
    /// Increment the value of mesh change index, see \ref changeIndex.
    inline void incChangeIndex()
    {
      changeIndex++;
    }

    /// Returns the mesh change index, see \ref changeIndex.
    inline long getChangeIndex()
    {
      return changeIndex;
    }

Thomas Witkowski's avatar
Thomas Witkowski committed
595
    ///
596
597
    void clearMacroFileInfo();

Thomas Witkowski's avatar
Thomas Witkowski committed
598
    ///
599
600
    int calcMemoryUsage();

601
602
603
    ///
    void deleteMeshStructure();

604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    /// In parallel computations the level of all macro elements is equal to the number
    /// of global pre refinements, \ref nParallelPreRefinements.
    inline int getMacroElementLevel()
    {
      return nParallelPreRefinements;
    }
#else
    /// In sequentiel computations the level of all macro elements is always 0.
    inline int getMacroElementLevel()
    {
      return 0;
    }
#endif

619
620
621
622
    /// Creates a map for all elements in mesh that maps from element indices
    /// to the corresponding pointers.
    void getElementIndexMap(map<int, Element*> &elIndexMap);

623
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
624
    ///
625
626
    static const Flag FILL_NOTHING;

Thomas Witkowski's avatar
Thomas Witkowski committed
627
    ///
628
    static const Flag FILL_COORDS; 
629

Thomas Witkowski's avatar
Thomas Witkowski committed
630
    ///
631
    static const Flag FILL_BOUND; 
632

Thomas Witkowski's avatar
Thomas Witkowski committed
633
    ///
634
    static const Flag FILL_NEIGH; 
635

Thomas Witkowski's avatar
Thomas Witkowski committed
636
    ///
637
    static const Flag FILL_OPP_COORDS; 
638

Thomas Witkowski's avatar
Thomas Witkowski committed
639
    ///
640
641
    static const Flag FILL_ORIENTATION; 

Thomas Witkowski's avatar
Thomas Witkowski committed
642
    ///
643
    static const Flag FILL_ADD_ALL; 
644
  
Thomas Witkowski's avatar
Thomas Witkowski committed
645
    ///
646
    static const Flag FILL_ANY_1D; 
647

Thomas Witkowski's avatar
Thomas Witkowski committed
648
    ///
649
    static const Flag FILL_ANY_2D; 
650

Thomas Witkowski's avatar
Thomas Witkowski committed
651
    ///
652
    static const Flag FILL_ANY_3D; 
653

Thomas Witkowski's avatar
Thomas Witkowski committed
654
    ///
655
    static const Flag FILL_DET;
656

Thomas Witkowski's avatar
Thomas Witkowski committed
657
    ///
658
659
660
661
662
663
    static const Flag FILL_GRD_LAMBDA;

    //**************************************************************************
    //  flags for Mesh traversal                                                
    //**************************************************************************

Thomas Witkowski's avatar
Thomas Witkowski committed
664
    ///
665
    static const Flag CALL_EVERY_EL_PREORDER;
666

Thomas Witkowski's avatar
Thomas Witkowski committed
667
    ///
668
    static const Flag CALL_EVERY_EL_INORDER;
669

Thomas Witkowski's avatar
Thomas Witkowski committed
670
    ///
671
    static const Flag CALL_EVERY_EL_POSTORDER;
672

Thomas Witkowski's avatar
Thomas Witkowski committed
673
    ///
674
    static const Flag CALL_LEAF_EL;
675

Thomas Witkowski's avatar
Thomas Witkowski committed
676
    ///
677
    static const Flag CALL_LEAF_EL_LEVEL;
678

Thomas Witkowski's avatar
Thomas Witkowski committed
679
    ///
680
    static const Flag CALL_EL_LEVEL;
681

Thomas Witkowski's avatar
Thomas Witkowski committed
682
    ///
683
684
    static const Flag CALL_MG_LEVEL;

685
686
687
    /// If set, left and right children are swapped in traverse.
    static const Flag CALL_REVERSE_MODE;

688
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
689
    ///
690
    bool findElementAtPointRecursive(ElInfo *elinfo,
691
				     const DimVec<double>& lambda,
692
				     int outside,
693
694
				     ElInfo *final_el_info);

695
696
697
698
699
700
701
702
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
    /** \brief
     * This functions is called in parallel computations by the function \ref
     * Mesh::initialize(). It checks that the macro file has enough macro elements
     * for the number of used processors and that all macro elements are of type 0.
     * If this is not the case, that macro mesh is globally refined in an
     * apropriate way and is written to a new macro file.
     *
703
704
     * The function overwrittes the macro and periodic filenames, if a new macro
     * fule was created for the current parallel usage.
705
     *
706
707
708
709
710
711
712
     * \param[in/out]  macroFilename      Name of the macro mesh file.
     * \param[in/out]  periodicFilename   If periodic boundaries are used, name of the
     *                                    periodicity file. Otherwise, the string must
     *                                    be empty.
     * \param[in]      check              If the mesh should be checked to be a correct
     *                                    AMDiS macro mesh, the value must be 1 and 0
     *                                    otherwise.
713
     */
714
715
    void checkParallelMacroFile(string &macroFilename, 
				string &periodicFilename,
716
717
718
				int check);
#endif

719
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
720
    /// maximal number of DOFs at one position
721
722
    static const int MAX_DOF;

Thomas Witkowski's avatar
Thomas Witkowski committed
723
    /// Name of this Mesh
724
    string name;
725

Thomas Witkowski's avatar
Thomas Witkowski committed
726
    /// Dimension of this Mesh. Doesn't have to be equal to dimension of world.
727
728
    int dim;

Thomas Witkowski's avatar
Thomas Witkowski committed
729
    /// Number of vertices in this Mesh
730
731
    int nVertices;

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    /// Number of Edges in this Mesh
733
734
    int nEdges;

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    /// Number of leaf elements in this Mesh
736
737
    int nLeaves;

Thomas Witkowski's avatar
Thomas Witkowski committed
738
    /// Total number of elements in this Mesh
739
740
    int nElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
741
    /// Number of faces in this Mesh
742
743
744
745
746
747
748
749
750
    int nFaces;

    /** \brief
     * Maximal number of elements that share one edge; used to allocate memory 
     * to store pointers to the neighbour at the refinement/coarsening edge 
     * (only 3d);
     */
    int maxEdgeNeigh;

Thomas Witkowski's avatar
Thomas Witkowski committed
751
    /// Diameter of the mesh in the DIM_OF_WORLD directions
752
753
754
755
756
757
758
759
760
761
    WorldVector<double> diam;

    /** \brief
     * Is pointer to NULL if mesh contains no parametric elements else pointer 
     * to a Parametric object containing coefficients of the parameterization 
     * and related information
     */
    Parametric *parametric;

    /** \brief
762
763
764
765
766
767
768
     * When an element is refined, not all dofs of the coarse element must be 
     * part of the new elements. An example are centered dofs when using higher
     * lagrange basis functions. The midpoint dof of the parents element is not
     * a dof of the both children elements. Therefore, the dof can be deleted. In
     * some situation, e.g., when using multigrid techniques, it can be necessary to
     * store this coarse dofs. Then this variable must be set to true. If false, the
     * not required coarse dofs will be deleted.
769
770
771
     */
    bool preserveCoarseDOFs;

Thomas Witkowski's avatar
Thomas Witkowski committed
772
    /// Number of all DOFs on a single element
773
    int nDofEl;
774
775
776
777
778

    /** \brief
     * Number of DOFs at the different positions VERTEX, EDGE, (FACE,) CENTER on
     * an element:
     *
779
     * - nDof[VERTEX]: number of DOFs at a vertex (>= 1)
780
     *
781
     * - nDof[EDGE]: number of DOFs at an edge; if no DOFs are associated to
782
783
     *   edges, then this value is 0
     *
784
     * - nDof[FACE]: number of DOFs at a face; if no DOFs are associated to
785
786
     *   faces, then this value is 0 (only 3d)
     *
787
     * - nDof[CENTER]: number of DOFs at the barycenter; if no DOFs are 
788
789
     *   associated to the barycenter, then this value is 0
     */
790
    DimVec<int> nDof;
791

792
    /// Number of nodes on a single element where DOFs are located. Needed for 
793
    /// the (de-) allocation of the DOF-vector on the element (\ref Element::dof).
794
    /// Here "node" is equivalent to the number of basis functions on the element.
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
    int nNodeEl;

    /** \brief
     * Gives the index of the first node at vertex, edge, face (only 3d), and 
     * barycenter:
     *
     * - node[VERTEX]: has always value 0; dof[0],...,dof[N_VERTICES-1] are 
     *   always DOFs at the vertices;
     *
     * - node[EDGE]: dof[node[EDGE]],..., dof[node[EDGE]+N_EDGES-1] are the DOFs
     *   at the N_EDGES edges, if DOFs are located at edges;
     *
     * - node[FACE]: dof[node[FACE]],..., dof[node[FACE]+N_FACES-1] are the DOFs
     *   at the N_FACES faces, if DOFs are located at faces (only 3d);
     *
     * - node[CENTER]: dof[node[CENTER]] are the DOFs at the barycenter, if DOFs
     *   are located at the barycenter;
     */
    DimVec<int> node;

Thomas Witkowski's avatar
Thomas Witkowski committed
815
    /// List of all DOFAdmins
816
    vector<DOFAdmin*> admin;
817

Thomas Witkowski's avatar
Thomas Witkowski committed
818
    /// List of all MacroElements of this Mesh
819
    deque<MacroElement*> macroElements;
820

Thomas Witkowski's avatar
Thomas Witkowski committed
821
    /// Used by check functions
822
    static vector<DegreeOfFreedom> dof_used;
823

824
825
826
827
828
829
830
831
832
833
    /** \brief
     * This map is used for serialization and deserialization of mesh elements. 
     * During the serialization process, all elements are visited and their dof indices
     * are written to the file. If a dof index at a position, i.e. vertex, line or face,
     * was written to file, the combination of dof index and position is inserted to
     * this map. That ensures that the same dof at the same position, but being part of
     * another element, is not written twice to the file. 
     * When a state should be deserialized, the information can be used to construct
     * exactly the same dof structure.
     */
834
    static map<pair<DegreeOfFreedom, int>, DegreeOfFreedom*> serializedDOFs;
835
836
837
838
839
840
841
842
843

    /** \brief
     * Used while mesh refinement. To create new elements 
     * elementPrototype->clone() is called, which returns a Element of the
     * same type as elementPrototype. So e.g. Elements of the different
     * dimensions can be created in a uniform way. 
     */
    Element* elementPrototype;

Thomas Witkowski's avatar
Thomas Witkowski committed
844
845
    /// Prototype for leaf data. Used for creation of new leaf data while 
    /// refinement.
846
847
    ElementData* elementDataPrototype;

Thomas Witkowski's avatar
Thomas Witkowski committed
848
    /// Used for enumeration of all mesh elements
849
850
    int elementIndex;

Thomas Witkowski's avatar
Thomas Witkowski committed
851
    /// True if the mesh is already initialized, false otherwise.
852
853
    bool initialized;

Thomas Witkowski's avatar
Thomas Witkowski committed
854
    /// Map of managed periodic vertex associations.
855
    map<BoundaryType, VertexVector*> periodicAssociations;
856

Thomas Witkowski's avatar
Thomas Witkowski committed
857
858
    /// If the mesh has been created by reading a macro file, here the information
    /// are stored about the content of the file.
859
    MacroInfo *macroFileInfo;
860

Thomas Witkowski's avatar
Thomas Witkowski committed
861
862
863
864
    /// This index is incremented every time the mesh is changed, e.g. by the 
    /// refinement or the coarsening manager. It can be used by other object if 
    /// the mesh has been changed by first copying this variable elsewhere and 
    /// comparing its values.
865
866
    long changeIndex;

867
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
868
869
870
    /// In parallel computations the mesh may be globally prerefined to achieve a
    /// fine enought starting mesh for the given number of ranks. The value of the
    /// variable will be defined in function \ref checkParallelMacroFile.
871
872
873
    int nParallelPreRefinements;
#endif

874
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
875
    /// for findElement-Fcts
876
    DimVec<double> final_lambda;
877
878

    /** \brief
879
880
     * Temporary variables that are used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
881
     */
882
    const WorldVector<double> *g_xy0, *g_xy;
883
884

    /** \brief
885
886
     * Temporary variable that is used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
887
888
     */    
    double *g_sp;
889
890
891
892
893
894
895
896
897
898
899
900
901
   
    friend class MacroInfo;
    friend class MacroReader;
    friend class MacroWriter;
    friend class MacroElement;
    friend class Element;
  };

}

#endif  // AMDIS_MESH_H