Dune Core Modules (2.3.1)

macrodata.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ALBERTA_MACRODATA_HH
4#define DUNE_ALBERTA_MACRODATA_HH
5
11#include <dune/common/fvector.hh>
12#include <dune/common/fmatrix.hh>
13
14#include <dune/grid/albertagrid/misc.hh>
15#include <dune/grid/albertagrid/algebra.hh>
16#include <dune/grid/albertagrid/albertaheader.hh>
17
18#if HAVE_ALBERTA
19
20namespace Dune
21{
22
23 namespace Alberta
24 {
25
26 template< int dim >
27 class MacroData
28 {
29 typedef MacroData< dim > This;
30
31 typedef ALBERTA MACRO_DATA Data;
32
33 static const int dimension = dim;
34 static const int numVertices = NumSubEntities< dimension, dimension >::value;
35 static const int numEdges = NumSubEntities< dimension, dimension-1 >::value;
36
37 static const int initialSize = 4096;
38
39 public:
40 template< int >
41 struct Library;
42
43 template< int > friend struct InstantiateMacroDataLibrary;
44
45 public:
46 typedef int ElementId[ numVertices ];
47
48 static const int supportPeriodicity = (DUNE_ALBERTA_VERSION >= 0x300);
49
50 MacroData ()
51 : data_( NULL ),
52 vertexCount_( -1 ),
53 elementCount_( -1 )
54 {}
55
56 operator Data * () const
57 {
58 return data_;
59 }
60
61 int vertexCount () const
62 {
63 return (vertexCount_ < 0 ? data_->n_total_vertices : vertexCount_);
64 }
65
66 int elementCount () const
67 {
68 return (elementCount_ < 0 ? data_->n_macro_elements : elementCount_);
69 }
70
71 ElementId &element ( int i ) const;
72 GlobalVector &vertex ( int i ) const;
73 int &neighbor ( int element, int i ) const;
74 BoundaryId &boundaryId ( int element, int i ) const;
75
80 void create ();
81
90 void finalize ();
91
100 void markLongestEdge ();
101
110 void setOrientation ( const Real orientation );
111
122 bool checkNeighbors () const;
123
125 void release ()
126 {
127 if( data_ != NULL )
128 {
129 ALBERTA free_macro_data( data_ );
130 data_ = NULL;
131 }
132 vertexCount_ = elementCount_ = -1;
133 }
134
140 int insertElement ( const ElementId &id );
141
147 int insertVertex ( const GlobalVector &coords )
148 {
149 assert( vertexCount_ >= 0 );
150 if( vertexCount_ >= data_->n_total_vertices )
151 resizeVertices( 2*vertexCount_ );
152 copy( coords, vertex( vertexCount_ ) );
153 return vertexCount_++;
154 }
155
161 int insertVertex ( const FieldVector< Real, dimWorld > &coords )
162 {
163 assert( vertexCount_ >= 0 );
164 if( vertexCount_ >= data_->n_total_vertices )
165 resizeVertices( 2*vertexCount_ );
166 copy( coords, vertex( vertexCount_ ) );
167 return vertexCount_++;
168 }
169
170 void insertWallTrafo ( const GlobalMatrix &m, const GlobalVector &t );
171 void insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
172 const FieldVector< Real, dimWorld > &shift );
173
174 void checkCycles ();
175
176 void read ( const std::string &filename, bool binary = false );
177
178 bool write ( const std::string &filename, bool binary = false ) const
179 {
180 if( binary )
181 return ALBERTA write_macro_data_xdr( data_, filename.c_str() );
182 else
183 return ALBERTA write_macro_data( data_, filename.c_str() );
184 }
185
186 private:
187 template< class Vector >
188 void copy ( const Vector &x, GlobalVector &y )
189 {
190 for( int i = 0; i < dimWorld; ++i )
191 y[ i ] = x[ i ];
192 }
193
194 void resizeElements ( const int newSize );
195
196 void resizeVertices ( const int newSize )
197 {
198 const int oldSize = data_->n_total_vertices;
199 data_->n_total_vertices = newSize;
200 data_->coords = memReAlloc< GlobalVector >( data_->coords, oldSize, newSize );
201 assert( (data_->coords != NULL) || (newSize == 0) );
202 }
203
204 private:
205 Data *data_;
206 int vertexCount_;
207 int elementCount_;
208 };
209
210
211
212 // MacroData::Library
213 // ------------------
214
215 template< int dim >
216 template< int >
217 struct MacroData< dim >::Library
218 {
219 typedef Alberta::MacroData< dim > MacroData;
220
221 static bool checkNeighbors ( const MacroData &macroData );
222 static void markLongestEdge ( MacroData &macroData );
223 static void setOrientation ( MacroData &macroData, const Real orientation );
224
225 private:
226 static Real edgeLength ( const MacroData &macroData, const ElementId &e, int edge );
227 static int longestEdge ( const MacroData &macroData, const ElementId &e );
228
229 template< class Type >
230 static void rotate ( Type *array, int i, int shift );
231
232 static void rotate ( MacroData &macroData, int i, int shift );
233 static void swap ( MacroData &macroData, int el, int v1, int v2 );
234 };
235
236
237
238 // Implementation of MacroData
239 // ---------------------------
240
241 template< int dim >
242 inline typename MacroData< dim >::ElementId &
243 MacroData< dim >::element ( int i ) const
244 {
245 assert( (i >= 0) && (i < data_->n_macro_elements) );
246 const int offset = i * numVertices;
247 return *reinterpret_cast< ElementId * >( data_->mel_vertices + offset );
248 }
249
250
251 template< int dim >
252 inline GlobalVector &MacroData< dim >::vertex ( int i ) const
253 {
254 assert( (i >= 0) && (i < data_->n_total_vertices) );
255 return data_->coords[ i ];
256 }
257
258
259 template< int dim >
260 inline int &MacroData< dim >::neighbor ( int element, int i ) const
261 {
262 assert( (element >= 0) && (element < data_->n_macro_elements) );
263 assert( (i >= 0) && (i < numVertices) );
264 return data_->neigh[ element*numVertices + i ];
265 }
266
267
268 template< int dim >
269 inline BoundaryId &MacroData< dim >::boundaryId ( int element, int i ) const
270 {
271 assert( (element >= 0) && (element < data_->n_macro_elements) );
272 assert( (i >= 0) && (i < numVertices) );
273 return data_->boundary[ element*numVertices + i ];
274 }
275
276
277#if DUNE_ALBERTA_VERSION >= 0x300
278 template< int dim >
279 inline void MacroData< dim >::create ()
280 {
281 release();
282 data_ = ALBERTA alloc_macro_data( dim, initialSize, initialSize );
283 data_->boundary = memAlloc< BoundaryId >( initialSize*numVertices );
284 if( dim == 3 )
285 data_->el_type = memAlloc< ElementType >( initialSize );
286 vertexCount_ = elementCount_ = 0;
287 elementCount_ = 0;
288 }
289#endif // #if DUNE_ALBERTA_VERSION >= 0x300
290
291#if DUNE_ALBERTA_VERSION == 0x200
292 template< int dim >
293 inline void MacroData< dim >::create ()
294 {
295 release();
296 data_ = ALBERTA alloc_macro_data( dim, initialSize, initialSize, 0 );
297 data_->boundary = memAlloc< BoundaryId >( initialSize*numVertices );
298 if( dim == 3 )
299 data_->el_type = memAlloc< ElementType >( initialSize );
300 vertexCount_ = elementCount_ = 0;
301 elementCount_ = 0;
302 }
303#endif // #if DUNE_ALBERTA_VERSION == 0x200
304
305
306 template< int dim >
307 inline void MacroData< dim >::finalize ()
308 {
309 if( (vertexCount_ >= 0) && (elementCount_ >= 0) )
310 {
311 resizeVertices( vertexCount_ );
312 resizeElements( elementCount_ );
313 ALBERTA compute_neigh_fast( data_ );
314
315 // assign default boundary id (if none is assigned)
316 for( int element = 0; element < elementCount_; ++element )
317 {
318 for( int i = 0; i < numVertices; ++i )
319 {
320 BoundaryId &id = boundaryId( element, i );
321 if( neighbor( element, i ) >= 0 )
322 {
323 assert( id == InteriorBoundary );
324 id = InteriorBoundary;
325 }
326 else
327 id = (id == InteriorBoundary ? DirichletBoundary : id);
328 }
329 }
330
331 vertexCount_ = elementCount_ = -1;
332 }
333 assert( (vertexCount_ < 0) && (elementCount_ < 0) );
334 }
335
336
337 template< int dim >
338 inline void MacroData< dim >::markLongestEdge ()
339 {
340 Library< dimWorld >::markLongestEdge( *this );
341 }
342
343
344 template< int dim >
345 inline void MacroData< dim >::setOrientation ( const Real orientation )
346 {
347 Library< dimWorld >::setOrientation( *this, orientation );
348 }
349
350
351 template< int dim >
352 inline bool MacroData< dim >::checkNeighbors () const
353 {
354 return Library< dimWorld >::checkNeighbors( *this );
355 }
356
357
358 template< int dim >
359 inline int MacroData< dim >::insertElement ( const ElementId &id )
360 {
361 assert( elementCount_ >= 0 );
362 if( elementCount_ >= data_->n_macro_elements )
363 resizeElements( 2*elementCount_ );
364
365 ElementId &e = element( elementCount_ );
366 for( int i = 0; i < numVertices; ++i )
367 {
368 e[ i ] = id[ i ];
369 boundaryId( elementCount_, i ) = InteriorBoundary;
370 }
371 if( dim == 3 )
372 data_->el_type[ elementCount_ ] = 0;
373
374 return elementCount_++;
375 }
376
377
378#if DUNE_ALBERTA_VERSION >= 0x300
379 template< int dim >
380 inline void MacroData< dim >
381 ::insertWallTrafo ( const GlobalMatrix &matrix, const GlobalVector &shift )
382 {
383 int &count = data_->n_wall_trafos;
384 AffineTransformation *&array = data_->wall_trafos;
385
386 // resize wall trafo array
387 array = memReAlloc< AffineTransformation >( array, count, count+1 );
388 assert( data_->wall_trafos != NULL );
389
390 // copy matrix and shift
391 for( int i = 0; i < dimWorld; ++i )
392 copy( matrix[ i ], array[ count ].M[ i ] );
393 copy( shift, array[ count ].t );
394 ++count;
395 }
396
397 template< int dim >
398 inline void MacroData< dim >
399 ::insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
400 const FieldVector< Real, dimWorld > &shift )
401 {
402 int &count = data_->n_wall_trafos;
403 AffineTransformation *&array = data_->wall_trafos;
404
405 // resize wall trafo array
406 array = memReAlloc< AffineTransformation >( array, count, count+1 );
407 assert( data_->wall_trafos != NULL );
408
409 // copy matrix and shift
410 for( int i = 0; i < dimWorld; ++i )
411 copy( matrix[ i ], array[ count ].M[ i ] );
412 copy( shift, array[ count ].t );
413 ++count;
414 }
415#endif // #if DUNE_ALBERTA_VERSION >= 0x300
416
417#if DUNE_ALBERTA_VERSION <= 0x200
418 template< int dim >
419 inline void MacroData< dim >
420 ::insertWallTrafo ( const GlobalMatrix &m, const GlobalVector &t )
421 {
422 DUNE_THROW( AlbertaError,
423 "Periodic grids are only supported in ALBERTA 3.0 or higher." );
424 }
425
426 template< int dim >
427 inline void MacroData< dim >
428 ::insertWallTrafo ( const FieldMatrix< Real, dimWorld, dimWorld > &matrix,
429 const FieldVector< Real, dimWorld > &shift )
430 {
431 DUNE_THROW( AlbertaError,
432 "Periodic grids are only supported in ALBERTA 3.0 or higher." );
433 }
434#endif // #if DUNE_ALBERTA_VERSION <= 0x200
435
436
437 template< int dim >
438 inline void MacroData< dim >::checkCycles ()
439 {
440 // ensure that the macro data has been finalized
441 finalize();
442 ALBERTA macro_test( data_, NULL );
443 }
444
445
446 template< int dim >
447 inline void MacroData< dim >::read ( const std::string &filename, bool binary )
448 {
449 release();
450 if( binary )
451 data_ = ALBERTA read_macro_xdr( filename.c_str() );
452 else
453 data_ = ALBERTA read_macro( filename.c_str() );
454 }
455
456
457 template< int dim >
458 inline void MacroData< dim >::resizeElements ( const int newSize )
459 {
460 const int oldSize = data_->n_macro_elements;
461 data_->n_macro_elements = newSize;
462 data_->mel_vertices = memReAlloc( data_->mel_vertices, oldSize*numVertices, newSize*numVertices );
463 data_->boundary = memReAlloc( data_->boundary, oldSize*numVertices, newSize*numVertices );
464 if( dim == 3 )
465 data_->el_type = memReAlloc( data_->el_type, oldSize, newSize );
466 assert( (newSize == 0) || (data_->mel_vertices != NULL) );
467 }
468
469 }
470
471}
472
473#endif // #if HAVE_ALBERTA
474
475#endif
fmatrix.hh
Implements a matrix constructed from a given type representing a field and compile-time given number ...
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
DUNE_THROW
#define DUNE_THROW(E, m)
Definition: exceptions.hh:244
Dune
Dune namespace.
Definition: alignment.hh:14
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)