Package sfc :: Package dofcode :: Module DofT
[hide private]
[frames] | no frames]

Source Code for Module sfc.dofcode.DofT

  1   
  2  header = """ 
  3  #ifndef SFC_DOFT_IS_INCLUDED 
  4  #define SFC_DOFT_IS_INCLUDED 
  5   
  6   
  7  #include <map> 
  8  #include <vector> 
  9  #include <iterator> 
 10  #include <iostream> 
 11  #include <stdexcept> 
 12  #include <tr1/memory> 
 13   
 14   
 15  namespace sfc 
 16  { 
 17    typedef std::pair< unsigned int, unsigned int > pair_ii; 
 18    typedef std::vector< unsigned int > vector_i; 
 19    typedef std::vector< std::pair<unsigned int, unsigned int> > vector_ii; 
 20   
 21   
 22    /* 
 23      A general degree of freedom (dof) builder for finite element applications, 
 24      counting global indices based on a given dof representation object of type D. 
 25      The type D must have an "bool D::operator<(const D&) const" implemented 
 26      to work as key in the std::map used here. 
 27   
 28      The structure dof2glob_map is the central structure, and is always constructed. 
 29      The structures loc2glob_array, loc2glob_map, glob2loc_map, and glob2dof_map are  
 30      optional, and enabled by the bools named create_loc2glob_array etc. 
 31      They are all initialized and updated by insert_dof(...). 
 32    */ 
 33    template <class D> 
 34    class DofT 
 35    { 
 36    protected: 
 37   
 38      unsigned int _global_dimension; 
 39      unsigned int _num_elements; 
 40      unsigned int _local_dimension; 
 41   
 42    public: // TODO: remove, temporary hack while changing sfc 
 43   
 44      /* 
 45         Central dynamic structure: based on dof representations 
 46         with type D, store the corresponding global indices. 
 47         (D Lj) -> (uint j) 
 48      */ 
 49      std::map< D, unsigned int > dof2glob_map; 
 50   
 51      /* 
 52         Static and efficient local to global mapping. 
 53         (uint e, uint i) -> (uint j) 
 54         Index like this: 
 55         global_index = loc2glob[element*local_dimension + local_index] 
 56      */ 
 57      bool create_loc2glob_array; 
 58      std::tr1::shared_ptr<unsigned int> loc2glob_array; 
 59   
 60      /* 
 61         Dynamic but slower local to global mapping. 
 62         (uint e, uint i) -> (uint j) 
 63      */ 
 64      bool create_loc2glob_map; 
 65      std::map< pair_ii, unsigned int > loc2glob_map; 
 66      // std::map< unsigned int e, vector_i > loc2glob_map; // TODO: this should be a bit more efficient for large local_dimension 
 67   
 68      /* 
 69         Dynamic global to local mapping. Given the global index of a dof, 
 70         provides a vector with the local index of the dof for each element 
 71         it occurs in. 
 72         (uint j) -> vector< pair<e1, i1>, ..  pair<en, in> > 
 73      */ 
 74      bool create_glob2loc_map; 
 75      std::map< unsigned int, vector_ii > glob2loc_map; 
 76   
 77      /* 
 78         Dynamic mapping from global index to dof representation object of type D. 
 79         (uint j) -> (D Lj) 
 80      */ 
 81      bool create_glob2dof_map; 
 82      std::map< unsigned int, D > glob2dof_map; 
 83   
 84   
 85    public: 
 86   
 87      DofT(bool create_glob2dof_map = false, bool create_glob2loc_map = false, 
 88           bool create_loc2glob_map = false, bool create_loc2glob_array = true): 
 89        _global_dimension(0), 
 90        _num_elements(0), 
 91        _local_dimension(0), 
 92        create_loc2glob_array( create_loc2glob_array ), 
 93        create_loc2glob_map( create_loc2glob_map ), 
 94        create_glob2loc_map( create_glob2loc_map ), 
 95        create_glob2dof_map( create_glob2dof_map ) 
 96      { 
 97      } 
 98   
 99      ~DofT() 
100      { 
101      } 
102   
103      // Clear all internal structures. 
104      void clear(); 
105   
106   
107      // Call to allocate static structures before using insert_dof(...). 
108      void init(unsigned int num_elements, unsigned int local_dimension); 
109   
110      // Call to allocate static structures before using insert_dof(...). 
111      void init(unsigned int num_elements, unsigned int local_dimension, std::tr1::shared_ptr<unsigned int> loc2glob); 
112   
113      // Update the internal structures with a new dof. 
114      unsigned int insert_dof(unsigned int e, unsigned int i, D Li); 
115   
116      // Get the local to global mapping array. 
117      std::tr1::shared_ptr<unsigned int> get_loc2glob_array() const 
118      { 
119        return loc2glob_array; 
120      } 
121   
122   
123      // Return the number of elements inserted. 
124      unsigned int num_elements() const 
125      { return _num_elements; } 
126   
127      // Return the number of global dofs inserted. 
128      unsigned int global_dimension() const 
129      { return _global_dimension; } 
130   
131      // Return the number of dofs per elements. 
132      unsigned int local_dimension() const 
133      { return _local_dimension; } 
134   
135   
136      // Return the global index for local dof i in element e. 
137      unsigned int glob_index(unsigned int e, unsigned int i) const; 
138   
139      // Return the global index for dof Lj represented with the templated type D. 
140      unsigned int glob_index(D Lj) const; 
141   
142      // Return the dof (represented with the templated type D) for global index j. 
143      D glob_dof(unsigned int j) const; 
144   
145      // Fill a vector with all the (element, index) pairs corresponding the dof with to this global index. 
146      void glob2loc(unsigned int j, vector_ii & loc) const; 
147   
148   
149      // Build loc2glob_array from loc2glob_map for future efficient lookup. 
150      void build_loc2glob(); 
151   
152    }; 
153   
154   
155    template <class D> 
156    void 
157    DofT<D>::clear() 
158    { 
159      _global_dimension = 0; 
160      _num_elements     = 0; 
161      _local_dimension  = 0; 
162   
163      dof2glob_map.clear(); 
164      loc2glob_array.reset(); 
165      loc2glob_map.clear(); 
166      glob2dof_map.clear(); 
167      glob2loc_map.clear(); 
168    } 
169   
170   
171    template<class T> 
172    class array_deleter 
173    { 
174    public: 
175      void operator() (T *p_) 
176      { 
177        delete [] p_; 
178      } 
179    }; 
180     
181    // constructing a shared_ptr to an array: 
182    template<class T> 
183    std::tr1::shared_ptr<T> create_shared_array(unsigned int n) 
184    { 
185      T * ptr = new T[n]; 
186      std::tr1::shared_ptr<T> sptr(ptr, array_deleter<T>()); 
187      return sptr; 
188    } 
189   
190   
191    template <class D> 
192    void 
193    DofT<D>::init(unsigned int num_elements, unsigned int local_dimension) 
194    { 
195      if( !create_loc2glob_array ) 
196      { 
197        std::cerr << "WARNING: Calling init with create_loc2glob_array == false makes little sense." << std::endl; 
198        std::cerr << "         Check your program logic." << std::endl; 
199      } 
200      create_loc2glob_array = true; 
201   
202      _num_elements    = num_elements; 
203      _local_dimension = local_dimension; 
204   
205      loc2glob_array = create_shared_array<unsigned int>(_num_elements*_local_dimension); 
206    } 
207   
208    template <class D> 
209    void 
210    DofT<D>::init(unsigned int num_elements, unsigned int local_dimension, std::tr1::shared_ptr<unsigned int> loc2glob) 
211    { 
212      if( !create_loc2glob_array ) 
213      { 
214        std::cerr << "WARNING: Calling init with create_loc2glob_array == false makes little sense." << std::endl; 
215        std::cerr << "         Check your program logic." << std::endl; 
216      } 
217      create_loc2glob_array = true; 
218   
219      _num_elements    = num_elements; 
220      _local_dimension = local_dimension; 
221   
222      loc2glob_array   = loc2glob; 
223    } 
224   
225    template <class D> 
226    unsigned int 
227    DofT<D>::insert_dof(unsigned int e, unsigned int i, D Li) 
228    { 
229      if (e+1 > _num_elements) 
230      { 
231          _num_elements = e+1; 
232          if( create_loc2glob_array ) 
233          { 
234              throw std::runtime_error("In this version of DofT we assume that the number of elements are predefined!"); 
235          } 
236      } 
237      if (i+1 > _local_dimension) 
238      { 
239          _local_dimension = i+1; 
240          if( create_loc2glob_array ) 
241          { 
242              throw std::runtime_error("In this version of DofT we assume that the local dimension is predefined!"); 
243          } 
244      } 
245   
246      // the return value 
247      unsigned int return_dof; 
248   
249      // check if the dof is new 
250      typename std::map< D, unsigned int >::iterator index_iter = dof2glob_map.find(Li); 
251   
252      if( index_iter != dof2glob_map.end() ) 
253      { 
254        // reuse global index 
255        return_dof = index_iter->second; 
256      } 
257      else 
258      { 
259        // pick a new global index 
260        return_dof = _global_dimension; 
261   
262        // count inserted global indices 
263        _global_dimension++; 
264   
265        // the central "D -> global index" map 
266        dof2glob_map[Li] = return_dof; 
267   
268        if ( create_glob2dof_map ) 
269        { 
270          std::pair<unsigned int, D> p(return_dof, Li); 
271          glob2dof_map.insert(p); 
272          //glob2dof_map[return_dof] = Li; 
273        } 
274   
275        if ( create_glob2loc_map ) 
276        { 
277          // initialize with empty vector 
278          glob2loc_map[return_dof] = vector_ii(); 
279          glob2loc_map[return_dof].reserve(_local_dimension); 
280        } 
281      } 
282   
283      if ( create_glob2loc_map ) 
284      { 
285        glob2loc_map[return_dof].push_back(pair_ii(e, i)); 
286      } 
287   
288      if( create_loc2glob_map ) 
289      { 
290        loc2glob_map[pair_ii(e, i)] = return_dof; 
291      } 
292   
293      if( create_loc2glob_array ) 
294      { 
295        unsigned int * l2g = loc2glob_array.get(); 
296        unsigned int k = e*_local_dimension + i; 
297        l2g[k] = return_dof; 
298      } 
299   
300      return return_dof; 
301    } 
302   
303   
304   
305    template <class D> 
306    unsigned int 
307    DofT<D>::glob_index(unsigned int e, unsigned int i) const 
308    { 
309      if ( create_loc2glob_array ) 
310      { 
311        if( e >= _num_elements ) 
312        { 
313          throw std::runtime_error("Invalid element index."); 
314        } 
315        if( i >= _local_dimension ) 
316        { 
317          throw std::runtime_error("Invalid local index."); 
318        } 
319        unsigned int * l2g = loc2glob_array.get(); 
320        return l2g[e*_local_dimension + i]; 
321      } 
322      else 
323      { 
324        if ( !create_loc2glob_map ) 
325        { 
326          throw std::runtime_error("This structure has not been created, you must turn on the create_loc2glob_map flag before initialization!"); 
327        } 
328       
329        pair_ii index(e, i); 
330        std::map< pair_ii, unsigned int >::const_iterator res = loc2glob_map.find(index); 
331         
332        if ( res == loc2glob_map.end() ) 
333        { 
334          throw std::runtime_error("In glob_index(e,i): Not found"); 
335        } 
336   
337        return res->second; 
338      } 
339    } 
340   
341   
342    template <class D> 
343    unsigned int 
344    DofT<D>::glob_index(D Lj) const 
345    { 
346      typename std::map< D, unsigned int >::const_iterator res = dof2glob_map.find(Lj); 
347   
348      if ( res == dof2glob_map.end() ) 
349      { 
350        throw std::runtime_error("In glob_index(Lj): Not found"); 
351      } 
352   
353      return res->second; 
354    } 
355   
356   
357    template <class D> 
358    D 
359    DofT<D>::glob_dof(unsigned int j) const 
360    { 
361      if ( !create_glob2dof_map ) 
362      { 
363        throw std::runtime_error("This structure has not been created, you must turn on the create_glob2dof_map flag before initialization!"); 
364      } 
365   
366      typename std::map<unsigned int, D>::const_iterator iter = glob2dof_map.find(j); 
367   
368      if ( iter == glob2dof_map.end() ) 
369      { 
370        //throw std::runtime_error("In glob_dof(j): Not found"); 
371        std::cerr << "In glob_dof(j): Not found" << std::endl; 
372        return D(); 
373      } 
374   
375      return iter->second; 
376    } 
377   
378   
379    template <class D> 
380    void 
381    DofT<D>::glob2loc(unsigned int j, vector_ii & loc) const 
382    { 
383      if ( !create_glob2loc_map ) 
384      { 
385        throw std::runtime_error("This structure has not been created, you must turn on the create_glob2loc_map flag before initialization!"); 
386      } 
387   
388      typename std::map<unsigned int, vector_ii>::const_iterator it = glob2loc_map.find(j); 
389      vector_ii::const_iterator b = it->second.begin(); 
390      vector_ii::const_iterator e = it->second.end(); 
391      loc.assign( b, e ); 
392    } 
393   
394   
395    template <class D> 
396    void 
397    DofT<D>::build_loc2glob() 
398    { 
399      if(create_loc2glob_array) 
400      { 
401        std::cerr << "WARNING: Calling build_loc2glob with create_loc2glob_array == true makes little sense." << std::endl; 
402        std::cerr << "         Doing nothing now. Check your program logic." << std::endl; 
403        return; 
404      } 
405   
406      if(!create_loc2glob_map) 
407      { 
408        throw std::runtime_error("This structure has not been created, you must turn on the create_loc2glob_map flag before initialization!"); 
409      } 
410   
411      loc2glob_array = create_shared_array<unsigned int>(_num_elements * _local_dimension); 
412   
413      typename std::map< pair_ii, unsigned int >::iterator iter; 
414      unsigned int * l2g = loc2glob_array.get(); 
415      for(iter = loc2glob_map.begin(); iter != loc2glob_map.end(); iter++) 
416      { 
417        unsigned int e = iter->first.first; 
418        unsigned int j = iter->first.second; 
419        unsigned int k = e * _local_dimension + j; 
420   
421        l2g[k] = iter->second; 
422      } 
423    } 
424   
425  } // namespace sfc 
426   
427  #endif 
428   
429  """ 
430   
431  implementation = "" 
432