DUNE-FEM (unstable)

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