Mesh.h 26 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
67
68
    virtual ~Mesh();

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
79
80
    Mesh& operator=(const Mesh&);

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

    /** \brief
     * Returns geometric information about this mesh. With GeoIndex p it is 
Backofen, Rainer's avatar
Backofen, Rainer committed
81
     * specified which information is requested.
82
     */
83
84
    inline int getGeo(GeoIndex p) const 
    { 
85
      return Global::getGeo(p, dim); 
86
    }
87

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

185
186
187
188
    /// 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. 
189
    DegreeOfFreedom *getDof(GeoIndex position);
190

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

197
198
199
    /// 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.
200
    const DOFAdmin* createDOFAdmin(string lname, DimVec<int> nDof);
201

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

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

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

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

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

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

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

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

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

253
254
255
256
257
258
    /** \} */

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

393
394
395
    /// Recalculates the number of leave elements.
    void updateNumberOfLeaves();

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

405
    /* \brief
406
407
408
     * 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.
409
     */
410
    void removeMacroElements(std::set<MacroElement*>& macros,
411
			     vector<const FiniteElemSpace*>& feSpaces);
412

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

Thomas Witkowski's avatar
Thomas Witkowski committed
416
    /// Used by \ref findElementAtPoint. 
417
418
419
    bool findElInfoAtPoint(const WorldVector<double>& xy,
			   ElInfo *el_info,
			   DimVec<double>& bary,
420
			   const MacroElement *start_mel,
421
422
			   const WorldVector<double> *xy0,
			   double *sp);
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
455
456
457
458

    /** \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.
     */
459
460
461
    bool findElementAtPoint(const WorldVector<double>& xy,
			    Element **elp, 
			    DimVec<double>& bary,
462
			    const MacroElement *start_mel,
463
464
			    const WorldVector<double> *xy0,
			    double *sp);
465

466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
    /** \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,
481
482
483
484
			   WorldVector<double>& coords)
    {
      return getDofIndexCoords(*dof, feSpace, coords);
    }
485

486
487
488
489
490

    /** \brief
     * This function is equal to \ref getDofIndexCoords as defined above, but takes
     * a DOF index instead of a DOF pointer.
     */
491
492
493
    bool getDofIndexCoords(DegreeOfFreedom dof, 
			   const FiniteElemSpace* feSpace,
			   WorldVector<double>& coords);
494

495
    /** \brief
496
497
498
     * 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.
499
     *
500
501
     * @param[in]   feSpace   The FE space to be used for the search.
     * @param[out]  coords    DOF vector that stores the coordinates to each DOF.
502
503
504
505
     */
    void getDofIndexCoords(const FiniteElemSpace* feSpace,
			   DOFVector<WorldVector<double> >& coords);

506
507
508
509
510
511
    /** \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.
     */
512
    void getAllDofs(const FiniteElemSpace *feSpace, 
513
514
		    std::set<const DegreeOfFreedom*>& allDofs);

Thomas Witkowski's avatar
Thomas Witkowski committed
515
    /// Returns FILL_ANY_?D
516
517
    inline static const Flag& getFillAnyFlag(int dim) 
    {
518
      switch (dim) {
519
520
521
522
523
524
525
526
527
528
529
530
531
      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;
      }
532
    }
533

Thomas Witkowski's avatar
Thomas Witkowski committed
534
    /// Serialize the mesh to a file.
535
    void serialize(ostream &out);
536

Thomas Witkowski's avatar
Thomas Witkowski committed
537
    /// Deserialize a mesh from a file.
538
    void deserialize(istream &in);
539

Thomas Witkowski's avatar
Thomas Witkowski committed
540
    /// Returns \ref elementIndex and increments it by 1.
541
542
    inline int getNextElementIndex() 
    { 
543
      return elementIndex++; 
544
    }
545

Thomas Witkowski's avatar
Thomas Witkowski committed
546
    /// Returns \ref initialized.
547
548
    inline bool isInitialized() 
    {
549
550
      return initialized; 
    }
551
  
Thomas Witkowski's avatar
Thomas Witkowski committed
552
    ///
553
    inline map<BoundaryType, VertexVector*>& getPeriodicAssociations() 
554
    {
555
      return periodicAssociations;
556
    }
557

558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
    /// 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
577
    ///
578
579
    bool associated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

Thomas Witkowski's avatar
Thomas Witkowski committed
580
    ///
581
582
    bool indirectlyAssociated(DegreeOfFreedom dof1, DegreeOfFreedom dof2);

Thomas Witkowski's avatar
Thomas Witkowski committed
583
    /// Returns \macroFileInfo
584
585
    inline MacroInfo* getMacroFileInfo() 
    { 
586
      return macroFileInfo;
587
    }
588

589
590
591
592
593
594
595
596
597
598
599
600
    /// 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
601
    ///
602
603
    void clearMacroFileInfo();

Thomas Witkowski's avatar
Thomas Witkowski committed
604
    ///
605
606
    int calcMemoryUsage();

607
608
609
    ///
    void deleteMeshStructure();

610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
#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

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

629
  public:
Thomas Witkowski's avatar
Thomas Witkowski committed
630
    ///
631
632
    static const Flag FILL_NOTHING;

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
642
    ///
643
    static const Flag FILL_OPP_COORDS; 
644

Thomas Witkowski's avatar
Thomas Witkowski committed
645
    ///
646
647
    static const Flag FILL_ORIENTATION; 

Thomas Witkowski's avatar
Thomas Witkowski committed
648
    ///
649
    static const Flag FILL_ADD_ALL; 
650
  
Thomas Witkowski's avatar
Thomas Witkowski committed
651
    ///
652
    static const Flag FILL_ANY_1D; 
653

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

Thomas Witkowski's avatar
Thomas Witkowski committed
657
    ///
658
    static const Flag FILL_ANY_3D; 
659

Thomas Witkowski's avatar
Thomas Witkowski committed
660
    ///
661
    static const Flag FILL_DET;
662

Thomas Witkowski's avatar
Thomas Witkowski committed
663
    ///
664
665
666
667
668
669
    static const Flag FILL_GRD_LAMBDA;

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

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

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

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

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

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

Thomas Witkowski's avatar
Thomas Witkowski committed
685
    ///
686
    static const Flag CALL_EL_LEVEL;
687

Thomas Witkowski's avatar
Thomas Witkowski committed
688
    ///
689
690
    static const Flag CALL_MG_LEVEL;

691
692
693
    /// If set, left and right children are swapped in traverse.
    static const Flag CALL_REVERSE_MODE;

694
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
695
    ///
696
    bool findElementAtPointRecursive(ElInfo *elinfo,
697
				     const DimVec<double>& lambda,
698
				     int outside,
699
700
				     ElInfo *final_el_info);

701
702
703
704
705
706
707
708
#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.
     *
709
710
     * The function overwrittes the macro and periodic filenames, if a new macro
     * fule was created for the current parallel usage.
711
     *
712
713
714
715
716
717
718
     * \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.
719
     */
720
721
    void checkParallelMacroFile(string &macroFilename, 
				string &periodicFilename,
722
723
724
				int check);
#endif

725
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
726
    /// maximal number of DOFs at one position
727
728
    static const int MAX_DOF;

Thomas Witkowski's avatar
Thomas Witkowski committed
729
    /// Name of this Mesh
730
    string name;
731

Thomas Witkowski's avatar
Thomas Witkowski committed
732
    /// Dimension of this Mesh. Doesn't have to be equal to dimension of world.
733
734
    int dim;

Thomas Witkowski's avatar
Thomas Witkowski committed
735
    /// Number of vertices in this Mesh
736
737
    int nVertices;

Thomas Witkowski's avatar
Thomas Witkowski committed
738
    /// Number of Edges in this Mesh
739
740
    int nEdges;

Thomas Witkowski's avatar
Thomas Witkowski committed
741
    /// Number of leaf elements in this Mesh
742
743
    int nLeaves;

Thomas Witkowski's avatar
Thomas Witkowski committed
744
    /// Total number of elements in this Mesh
745
746
    int nElements;

Thomas Witkowski's avatar
Thomas Witkowski committed
747
    /// Number of faces in this Mesh
748
749
750
751
752
753
754
755
756
    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
757
    /// Diameter of the mesh in the DIM_OF_WORLD directions
758
759
760
761
762
763
764
765
766
767
    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
768
769
770
771
772
773
774
     * 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.
775
776
777
     */
    bool preserveCoarseDOFs;

Thomas Witkowski's avatar
Thomas Witkowski committed
778
    /// Number of all DOFs on a single element
779
    int nDofEl;
780
781
782
783
784

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

798
    /// Number of nodes on a single element where DOFs are located. Needed for 
799
    /// the (de-) allocation of the DOF-vector on the element (\ref Element::dof).
800
    /// Here "node" is equivalent to the number of basis functions on the element.
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
    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
821
    /// List of all DOFAdmins
822
    vector<DOFAdmin*> admin;
823

Thomas Witkowski's avatar
Thomas Witkowski committed
824
    /// List of all MacroElements of this Mesh
825
    deque<MacroElement*> macroElements;
826

Thomas Witkowski's avatar
Thomas Witkowski committed
827
    /// Used by check functions
828
    static vector<DegreeOfFreedom> dof_used;
829

830
831
832
833
834
835
836
837
838
839
    /** \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.
     */
840
    static map<pair<DegreeOfFreedom, int>, DegreeOfFreedom*> serializedDOFs;
841
842
843
844
845
846
847
848
849

    /** \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
850
851
    /// Prototype for leaf data. Used for creation of new leaf data while 
    /// refinement.
852
853
    ElementData* elementDataPrototype;

Thomas Witkowski's avatar
Thomas Witkowski committed
854
    /// Used for enumeration of all mesh elements
855
856
    int elementIndex;

Thomas Witkowski's avatar
Thomas Witkowski committed
857
    /// True if the mesh is already initialized, false otherwise.
858
859
    bool initialized;

Thomas Witkowski's avatar
Thomas Witkowski committed
860
    /// Map of managed periodic vertex associations.
861
    map<BoundaryType, VertexVector*> periodicAssociations;
862

Thomas Witkowski's avatar
Thomas Witkowski committed
863
864
    /// If the mesh has been created by reading a macro file, here the information
    /// are stored about the content of the file.
865
    MacroInfo *macroFileInfo;
866

Thomas Witkowski's avatar
Thomas Witkowski committed
867
868
869
870
    /// 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.
871
872
    long changeIndex;

873
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Thomas Witkowski's avatar
Thomas Witkowski committed
874
875
876
    /// 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.
877
878
879
    int nParallelPreRefinements;
#endif

880
  protected:
Thomas Witkowski's avatar
Thomas Witkowski committed
881
    /// for findElement-Fcts
882
    DimVec<double> final_lambda;
883
884

    /** \brief
885
886
     * Temporary variables that are used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
887
     */
888
    const WorldVector<double> *g_xy0, *g_xy;
889
890

    /** \brief
891
892
     * Temporary variable that is used in functions \ref fineElInfoatPoint and
     * \ref fineElementAtPointRecursive.
893
894
     */    
    double *g_sp;
895
896
897
898
899
900
901
902
903
904
905
906
907
   
    friend class MacroInfo;
    friend class MacroReader;
    friend class MacroWriter;
    friend class MacroElement;
    friend class Element;
  };

}

#endif  // AMDIS_MESH_H