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