1 /*--------------------------------------------------------------------------*/
2 /* ALBERTA: an Adaptive multi Level finite element toolbox using */
3 /* Bisectioning refinement and Error control by Residual */
4 /* Techniques for scientific Applications */
5 /* */
6 /* file: memory.c */
7 /* */
8 /* description: routines for memory handling */
9 /* */
10 /*--------------------------------------------------------------------------*/
11 /* */
12 /* authors: Alfred Schmidt */
13 /* Zentrum fuer Technomathematik */
14 /* Fachbereich 3 Mathematik/Informatik */
15 /* Universitaet Bremen */
16 /* Bibliothekstr. 2 */
17 /* D-28359 Bremen, Germany */
18 /* */
19 /* Kunibert G. Siebert */
20 /* Institut fuer Mathematik */
21 /* Universitaet Augsburg */
22 /* Universitaetsstr. 14 */
23 /* D-86159 Augsburg, Germany */
24 /* */
25 /* Daniel Koester */
26 /* Institut fuer Mathematik */
27 /* Albert-Ludwigs-Universitaet Freiburg */
28 /* Hermann-Herder-Str. 10 */
29 /* D-79104 Freiburg */
30 /* */
31 /* Claus-Justus Heine */
32 /* Abteilung fuer Angewandte Mathematik */
33 /* Albert-Ludwigs-Universitaet Freiburg */
34 /* Hermann-Herder-Str. 10 */
35 /* D-79104 Freiburg im Breisgau, Germany */
36 /* */
37 /* http://www.alberta-fem.de */
38 /* */
39 /* (c) by A. Schmidt and K.G. Siebert (1996-2003) */
40 /* D. Koester (2005), */
41 /* C.-J. Heine (2006-2007) */
42 /* */
43 /*--------------------------------------------------------------------------*/
44
45 #ifdef HAVE_CONFIG_H
46 # include "config.h"
47 #endif
48
49 #include "alberta_intern.h"
50 #include "alberta.h"
51
52 /*--------------------------------------------------------------------------*/
53 /* some routines for allocation and deallocation */
54 /*--------------------------------------------------------------------------*/
55
56 #define ALBERTA_DEBUG_MEMORY 0
57
58 /* #define ALBERTA_USE_MALLOC 1 */
59
60 #ifndef ALBERTA_USE_MALLOC
61 # define ALBERTA_USE_MALLOC \
62 (!ALBERTA_DEBUG_MEMORY && ALBERTA_DEBUG && ALBERTA_ALLOC_RECORD)
63 #endif
64
65 #if ALBERTA_DEBUG_MEMORY
66 # define MEMMAGIC "ABMM"
67 #endif
68
69 /* default value for increasing size of memory in an memoryadmin. see
70 * also setCapacityIncrement().
71 */
72 #if ALBERTA_DEBUG_MEMORY > 1
73 # define DEFAULTCAPACITYINCREMENT 5
74 #else
75 # define DEFAULTCAPACITYINCREMENT 1000
76 #endif
77
78
79 /* Alignment by integer division and multiplication */
80 #define ALIGNMEM(size, align) \
81 (((size_t)(size) + (size_t)(align) - 1) / (size_t)(align) * (size_t)(align))
82
83 #if ALBERTA_DEBUG_MEMORY
84 /* During memory debugging every memory block also either belongs to a
85 * second doubly linked free-list, or to a doubly link alloc-list. We
86 * maintain additionally a status flags (allocated or not), and a 4
87 * byte magic number just below and above the memory block to catch
88 * array bound violation.
89 */
90 typedef enum MemoryStatusEnum {
91 MemFree = 1,
92 MemAllocated = 2,
93 MemChecked = 4, /* used to tag list members while performing sanity checks */
94 } MEMORY_STATUS;
95
96 typedef struct MemoryInfo MEMORY_INFO;
97 struct MemoryInfo
98 {
99 MEMORY_INFO *prev;
100 MEMORY_INFO *next;
101 void *admin;
102 MEMORY_STATUS status;
103 };
104
105 /* Two global lists for all memory blocks */
106 static struct {
107 MEMORY_INFO *next;
108 MEMORY_INFO *prev;
109 } freeList = { (void *)&freeList, (void *)&freeList };
110
111 static struct {
112 MEMORY_INFO *next;
113 MEMORY_INFO *prev;
114 } allocList = { (void *)&allocList, (void *)&allocList };
115
116 #endif
117
118 /* Linked list of free memory blocks of size
119 * "MEMORYADMIN->objectsSize". The list pointer "next" occupies the
120 * first sizeof(void *) bytes of an allocated block (so it is
121 * destroyed during allocation).
122 */
123 typedef struct FreeMemory FREEMEMORY;
124 struct FreeMemory
125 {
126 FREEMEMORY *next;
127 };
128
129 /* Linked list of blocks. "next" and "end" are in front of the actual
130 * memory region which starts just behing "&end". "end" itself is the
131 * address just after the memory region attached behing "&end".
132 */
133 typedef struct Block BLOCK;
134 struct Block
135 {
136 BLOCK *next;
137 void *start; /* The actual start of the _aligned_ data area */
138 void *end; /* Points past the end of the data area. */
139 size_t size; /* The acutal _UNaligned_ size */
140 };
141
142 /* The free-lits "freeMem" may contain blocks from different
143 * memory-block regions.
144 */
145 typedef struct MemoryAdmin MEMORYADMIN;
146 struct MemoryAdmin
147 {
148 const char *name;
149 unsigned int capacity;
150 unsigned int capacityIncrement;
151 size_t alignment; /* the alignment of each object */
152 size_t objectSize; /* the aligned size of each object */
153 #if ALBERTA_DEBUG_MEMORY
154 size_t prefixSize; /* aligned size of data stored in front */
155 size_t suffixOffset; /* offset to trailing data */
156 #endif
157 BLOCK *blocks;
158 FREEMEMORY *freeMem;
159 };
160
161 #if ALBERTA_DEBUG_MEMORY
memListAdd(MEMORY_INFO * head,MEMORY_INFO * new)162 static void memListAdd(MEMORY_INFO *head, MEMORY_INFO *new)
163 {
164 head->next->prev = new;
165 new->next = head->next;
166 new->prev = head;
167 head->next = new;
168 }
169
memListDel(MEMORY_INFO * old)170 static void memListDel(MEMORY_INFO *old)
171 {
172 old->prev->next = old->next;
173 old->next->prev = old->prev;
174 }
175
adminDelFromMemList(MEMORY_INFO * head,MEMORYADMIN * adm)176 static void adminDelFromMemList(MEMORY_INFO *head, MEMORYADMIN *adm)
177 {
178 MEMORY_INFO *memInfo;
179
180 for (memInfo = head->next; memInfo != head; memInfo = memInfo->next) {
181 if (memInfo->admin == adm) {
182 memListDel(memInfo);
183 }
184 }
185 }
186
_AI_memListCheck(MEMORY_INFO * head,MEMORY_STATUS status)187 void _AI_memListCheck(MEMORY_INFO *head,
188 MEMORY_STATUS status)
189 {
190 FUNCNAME("_AI_memListCheck");
191 MEMORY_INFO *memInfo;
192 MEMORYADMIN *adm;
193 int garbled = 0;
194
195 for (memInfo = head->next; memInfo != head; memInfo = memInfo->next) {
196 adm = memInfo->admin;
197 if (memInfo->status != status) {
198 MSG("Block: %p (%p): Wrong status: %d, should be %d.\n",
199 memInfo, (char *)memInfo + adm->prefixSize,
200 memInfo->status, status);
201 garbled = 1;
202 }
203 if (memcmp((char *)memInfo + adm->prefixSize - strlen(MEMMAGIC),
204 MEMMAGIC, strlen(MEMMAGIC)) != 0) {
205 MSG("Block: %p (%p): Wrong magic in front of data area.\n",
206 memInfo, (char *)memInfo + adm->prefixSize);
207 garbled = 1;
208 }
209 if (memcmp((char *)memInfo + adm->suffixOffset,
210 MEMMAGIC, strlen(MEMMAGIC)) != 0) {
211 MSG("Block: %p (%p): Wrong magic after data area.\n",
212 memInfo, (char *)memInfo + adm->prefixSize);
213 garbled = 1;
214 }
215 if (status == MemFree) {
216 FREEMEMORY *freeMem = (FREEMEMORY *)((char *)memInfo + adm->prefixSize);
217 FREEMEMORY *pos;
218
219 for (pos = adm->freeMem; pos; pos = pos->next) {
220 if (pos == freeMem) {
221 break;
222 }
223 }
224 if (pos == NULL) {
225 MSG("Block: %p (%p): Free block not on admin's free-list.\n",
226 memInfo, (char *)memInfo + adm->prefixSize);
227 garbled = 1;
228 }
229 }
230 if (garbled) {
231 ERROR_EXIT("Memory management has been garbled.\n");
232 }
233 }
234 }
235
_AI_freeListCheck(void)236 void _AI_freeListCheck(void)
237 {
238 _AI_memListCheck((void *)&freeList, MemFree);
239 }
_AI_allocListCheck(void)240 void _AI_allocListCheck(void)
241 {
242 _AI_memListCheck((void *)&allocList, MemAllocated);
243 }
244
245 #endif
246
newBlock(MEMORYADMIN * ma,unsigned int capacityIncrement)247 static int newBlock(MEMORYADMIN* ma, unsigned int capacityIncrement)
248 {
249 FUNCNAME("newBlock");
250 #if ALBERTA_DEBUG_MEMORY
251 char *thisChunk;
252 #endif
253 BLOCK *block;
254 FREEMEMORY *freeMem;
255 int i;
256
257 #if ALBERTA_DEBUG_MEMORY
258 _AI_freeListCheck();
259 _AI_allocListCheck();
260 #endif
261
262 /* In the worst case we need ma->alignment-1 additional bytes to
263 * satisfy the memery alignment requirement.
264 */
265 block =
266 (BLOCK *)alberta_alloc(sizeof(BLOCK) +
267 ma->alignment - 1 +
268 capacityIncrement * ma->objectSize,
269 funcName, __FILE__, __LINE__);
270
271 /* compute the bounds for the _aligned_ data region */
272 block->size =
273 sizeof(BLOCK) + ma->alignment - 1 + capacityIncrement * ma->objectSize;
274 block->start = (void *)ALIGNMEM(block+1, ma->alignment);
275 block->end =
276 (void *)((char *)block->start + capacityIncrement * ma->objectSize);
277
278 /* now build up the free-list */
279 #if !ALBERTA_DEBUG_MEMORY
280 for (freeMem = (FREEMEMORY *)block->start, i = 0;
281 i < (int)(capacityIncrement - 1);
282 freeMem = (FREEMEMORY *)((char *)freeMem + ma->objectSize), i++) {
283 freeMem->next = (FREEMEMORY *)((char *)freeMem + ma->objectSize);
284 }
285 freeMem->next = ma->freeMem;
286 ma->freeMem = (FREEMEMORY *)block->start;
287 #else
288 /* build up the doubly linked free-mem list */
289 for (thisChunk = block->start;
290 thisChunk < (char *)block->end;
291 thisChunk += ma->objectSize) {
292 MEMORY_INFO *memInfo;
293
294 memInfo = (MEMORY_INFO *)thisChunk;
295 memInfo->admin = (void *)ma;
296 memInfo->status = MemFree;
297
298 /* add the magic numbers in front and after the data region */
299 memcpy(thisChunk + ma->prefixSize - strlen(MEMMAGIC),
300 MEMMAGIC, strlen(MEMMAGIC));
301 memcpy(thisChunk + ma->suffixOffset, MEMMAGIC, strlen(MEMMAGIC));
302
303 memListAdd((void *)&freeList, memInfo);
304 }
305
306 /* build up the ordinary singly linked list */
307 for (freeMem = (FREEMEMORY *)((char *)block->start + ma->prefixSize), i = 0;
308 i < capacityIncrement - 1;
309 freeMem = (FREEMEMORY *)((char *)freeMem + ma->objectSize), i++) {
310 freeMem->next = (FREEMEMORY *)((char *)freeMem + ma->objectSize);
311 }
312 freeMem->next = ma->freeMem;
313 ma->freeMem = (FREEMEMORY *)((char *)block->start + ma->prefixSize);
314
315 _AI_freeListCheck();
316 _AI_allocListCheck();
317 #endif
318
319 /* hook the new block into the block list of the memory admin */
320 ma->capacity += capacityIncrement;
321 block->next = ma->blocks;
322 ma->blocks = block;
323
324 return 0;
325 }
326
327 /*
328 * newObject(..) returns a void-pointer to a memoryadmin for objects of size
329 * objectSize. Use the returned value as parameter for getMemory.
330 */
_newObject(size_t objectSize,size_t alignment,unsigned int initialCapacity,const char * name)331 static void *_newObject(size_t objectSize,
332 size_t alignment,
333 unsigned int initialCapacity,
334 const char *name)
335 {
336 FUNCNAME("newObject");
337 MEMORYADMIN *ma;
338
339 TEST_EXIT(objectSize, "Attempted to allocate a zero length object!\n");
340
341 ma = (MEMORYADMIN *)alberta_alloc(sizeof(MEMORYADMIN),
342 funcName, __FILE__, __LINE__);
343
344 ma->name = name ? strdup(name) : NULL;
345
346 ma->capacity = 0;
347 ma->capacityIncrement =
348 (initialCapacity ? initialCapacity : DEFAULTCAPACITYINCREMENT);
349 ma->alignment = (alignment ? alignment : objectSize);
350
351 if (ma->alignment > 2*sizeof(REAL)) {
352 WARNING("large alignment %d requested.\n", ma->alignment);
353 }
354
355 /* Need at least space for one pointer ... better do not allocate
356 * integers with this structure ...
357 */
358 objectSize = MAX(sizeof(FREEMEMORY), objectSize);
359 #if !ALBERTA_DEBUG_MEMORY
360 ma->objectSize = ALIGNMEM(objectSize, ma->alignment);
361 #else
362 /* allocate enough space to include that magic suffix */
363 ma->objectSize = ALIGNMEM(objectSize+strlen(MEMMAGIC), ma->alignment);
364
365 /* data stored in front of each object */
366 ma->prefixSize = ALIGNMEM(sizeof(MEMORY_INFO) + strlen(MEMMAGIC),
367 ma->alignment);
368 ma->objectSize += ma->prefixSize;
369 ma->suffixOffset = ma->prefixSize + objectSize;
370 #endif
371 ma->blocks = NULL;
372 ma->freeMem = NULL;
373
374 if(initialCapacity) {
375 newBlock(ma, initialCapacity);
376 }
377
378 return ma;
379 }
380
381 /* Default case: the required alignment is the object size if <
382 * sizeof(void *), otherwise it is the sizeof(void *).
383 */
newObject(size_t objectSize,unsigned int initialCapacity,const char * name)384 static void *newObject(size_t objectSize, unsigned int initialCapacity,
385 const char *name)
386 {
387 return _newObject(objectSize,
388 objectSize < sizeof(void *) ? objectSize : sizeof(void *),
389 initialCapacity, name);
390 }
391
392 /* returns a pointer to an object of size objectsize given by creating the
393 * memoryadmin.
394 */
getMemory(void * memoryAdmin)395 static void* getMemory(void* memoryAdmin)
396 {
397 #if ALBERTA_DEBUG_MEMORY
398 FUNCNAME("getMemory");
399 MEMORY_INFO *memInfo;
400 char *thisChunk;
401 #endif
402 MEMORYADMIN *admin;
403 FREEMEMORY *tmp;
404
405 admin = (MEMORYADMIN *)memoryAdmin;
406
407 #if ALBERTA_USE_MALLOC
408 return malloc(admin->objectSize);
409 #endif
410
411 while(admin->capacity <= 0)
412 newBlock(admin, admin->capacityIncrement);
413
414 #if ALBERTA_DEBUG_MEMORY
415 /* Check the magic numbers, check the status, remove the block
416 * from the additional freeList and hook it in to the allocList,
417 * change its status.
418 */
419 thisChunk = (char *)admin->freeMem - admin->prefixSize;
420 memInfo = (MEMORY_INFO *)thisChunk;
421 TEST_EXIT(memInfo->next != memInfo->prev,
422 "Unconnected free block.\n");
423 TEST_EXIT(memInfo->admin == admin,
424 "memInfo->admin != admin.\n");
425 TEST_EXIT(memInfo->status == MemFree,
426 "Status %d != MemFree (= %d).\n", memInfo->status, MemFree);
427 TEST_EXIT(memcmp(thisChunk + admin->prefixSize - strlen(MEMMAGIC),
428 MEMMAGIC, strlen(MEMMAGIC)) == 0,
429 "Wrong magic in front of free block.\n");
430 TEST_EXIT(memcmp(thisChunk + admin->suffixOffset,
431 MEMMAGIC, strlen(MEMMAGIC)) == 0,
432 "Wrong magic after free block.\n");
433
434 /* remove block from freeList */
435 memListDel(memInfo);
436 memListAdd((void *)&allocList, memInfo);
437 memInfo->status = MemAllocated;
438 #endif
439
440 tmp = admin->freeMem;
441 admin->freeMem = admin->freeMem->next;
442
443 admin->capacity--;
444
445 #ifdef __GNUC__
446 /* gcc has a bug when run with -O3 and the i386 default target; so
447 * insert an optimization barrier here :(
448 */
449 __asm__ __volatile__("" ::: "memory");
450 #endif
451
452 return (void *)tmp;
453 }
454
455 /* returns memory at *object to memoryadmin. !!! object must be the
456 * result of an call of getMemory(...) with the same memoryadmin. !!!
457 */
freeMemory(void * object,void * memoryAdmin)458 static void freeMemory(void *object, void *memoryAdmin)
459 {
460 #if ALBERTA_DEBUG_MEMORY
461 FUNCNAME("freeMemory");
462 MEMORY_INFO *memInfo;
463 char *thisChunk;
464 #endif
465 MEMORYADMIN *memadm = (MEMORYADMIN *)memoryAdmin;
466 FREEMEMORY *freemem = (FREEMEMORY *)object;
467
468 #if ALBERTA_USE_MALLOC
469 free(object);
470 return;
471 #endif
472
473 #if ALBERTA_DEBUG_MEMORY
474 _AI_freeListCheck();
475 _AI_allocListCheck();
476
477 thisChunk = (char *)object - memadm->prefixSize;
478 memInfo = (MEMORY_INFO *)thisChunk;
479
480 TEST_EXIT(memInfo->next != memInfo->prev,
481 "Unconnected allocated block.\n");
482 TEST_EXIT(memInfo->admin == memadm,
483 "memInfo->admin != memadm.\n");
484 TEST_EXIT(memInfo->status == MemAllocated,
485 "Status %d != MemAllocated (= %d).\n",
486 memInfo->status, MemAllocated);
487 TEST_EXIT(memcmp(thisChunk + memadm->prefixSize - 4,
488 MEMMAGIC, strlen(MEMMAGIC)) == 0,
489 "Wrong magic in front of allocated block.\n");
490 TEST_EXIT(memcmp(thisChunk + memadm->suffixOffset,
491 MEMMAGIC, strlen(MEMMAGIC)) == 0,
492 "Wrong magic after allocated block.\n");
493
494 memListDel(memInfo);
495 memListAdd((void *)&freeList, memInfo);
496 memInfo->status = MemFree;
497 #endif
498
499 freemem->next = memadm->freeMem;
500 memadm->freeMem = freemem;
501 memadm->capacity ++;
502 }
503
504 /* the counterpart to newObject, frees all memory allocated with the given
505 * memoryadmin.
506 */
deleteObject(void * memoryAdmin)507 static void deleteObject(void* memoryAdmin)
508 {
509 #if ALBERTA_DEBUG_MEMORY
510 FUNCNAME("deleteObject");
511 #endif
512 BLOCK *b, *tmp;
513 MEMORYADMIN *adm = (MEMORYADMIN *)memoryAdmin;
514 size_t size;
515
516 DEBUG_TEST_EXIT(memoryAdmin, "memoryAdmin == NULL\n");
517
518 #if ALBERTA_DEBUG_MEMORY
519 /* inefficient like hell, but this is debug code anyway ... */
520 adminDelFromMemList((void *)&freeList, adm);
521 _AI_freeListCheck();
522 adminDelFromMemList((void *)&allocList, adm);
523 _AI_allocListCheck();
524 #endif
525
526 b = adm->blocks;
527 while(b != NULL) {
528 tmp = b->next;
529 size = b->size;
530 #if ALBERTA_DEBUG_MEMORY
531 memset(b, -1, size);
532 #endif
533 alberta_free(b, size);
534 b = tmp;
535 }
536 if (adm->name)
537 free((char *)adm->name);
538 alberta_free(adm, sizeof(MEMORYADMIN));
539 }
540
541 typedef struct dof_admin_mem_info DOF_ADMIN_MEM_INFO;
542 struct dof_admin_mem_info
543 {
544 void *dof_matrix;
545 void *real_matrix_row;
546 void *real_d_matrix_row;
547 void *real_dd_matrix_row;
548 void *dof_int_vec;
549 void *dof_dof_vec;
550 void *int_dof_vec;
551 void *dof_uchar_vec;
552 void *dof_schar_vec;
553 void *dof_real_vec;
554 void *dof_real_d_vec;
555 void *dof_real_dd_vec;
556 void *dof_ptr_vec;
557 };
558
559
560 /*--------------------------------------------------------------------------*/
561 /* memory management for DOF admin structures */
562 /*--------------------------------------------------------------------------*/
563
add_dof_admin_to_mesh(DOF_ADMIN * admin,MESH * mesh)564 static void add_dof_admin_to_mesh(DOF_ADMIN *admin, MESH *mesh)
565 {
566 FUNCNAME("add_dof_admin_to_mesh");
567 int i, n, dim = mesh->dim;
568
569 admin->mesh = mesh;
570 n = mesh->n_dof_admin;
571 if ((n > 0) && (mesh->dof_admin == NULL))
572 ERROR_EXIT("no mesh->dof_admin but n_dof_admin=%d\n", n);
573 if ((n <= 0) && (mesh->dof_admin != NULL))
574 ERROR_EXIT("found mesh->dof_admin but n_dof_admin=%d\n", n);
575
576 for (i = 0; i < n; i++)
577 if (mesh->dof_admin[i] == admin)
578 ERROR_EXIT("admin %s is already associated to mesh %s\n",
579 NAME(admin), NAME(mesh));
580
581 mesh->dof_admin = MEM_REALLOC(mesh->dof_admin, n, n+1, DOF_ADMIN *);
582 n++;
583
584 mesh->dof_admin[n-1] = admin;
585 mesh->n_dof_admin = n;
586
587 mesh->n_dof_el = 0;
588
589 admin->n0_dof[VERTEX] = mesh->n_dof[VERTEX];
590 mesh->n_dof[VERTEX] += admin->n_dof[VERTEX];
591 mesh->n_dof_el += N_VERTICES(dim) * mesh->n_dof[VERTEX];
592
593 admin->n0_dof[CENTER] = mesh->n_dof[CENTER];
594 mesh->n_dof[CENTER] += admin->n_dof[CENTER];
595 mesh->n_dof_el += mesh->n_dof[CENTER];
596
597 if(dim > 1) {
598 admin->n0_dof[EDGE] = mesh->n_dof[EDGE];
599 mesh->n_dof[EDGE] += admin->n_dof[EDGE];
600 mesh->n_dof_el += N_EDGES(dim) * mesh->n_dof[EDGE];
601 }
602
603 if(dim == 3) {
604 admin->n0_dof[FACE] = mesh->n_dof[FACE];
605 mesh->n_dof[FACE] += admin->n_dof[FACE];
606 mesh->n_dof_el += N_FACES_3D * mesh->n_dof[FACE];
607 }
608
609 #define REALLY_DONT_USE_THIS_CODE true
610 #if !REALLY_DONT_USE_THIS_CODE
611 /* It is not so easy: the ordering of the node-types in the
612 * DOF-pointers in the elements (i.e. el->dof[][][]) is still the
613 * same as before, so: first come the vertexc nodes, then come the
614 * edge nodes, then the face nodes and last the center node.
615 */
616 mesh->n_node_el = 0;
617 for (i = 0; i < N_NODE_TYPES; i++) {
618 mesh->node[i] = mesh->n_node_el;
619 if (mesh->n_dof[i] > 0) {
620 mesh->n_node_el += n_nodes[i];
621 }
622 }
623 #else
624 mesh->node[VERTEX] = 0;
625 if (mesh->n_dof[VERTEX] > 0) {
626 mesh->n_node_el = N_VERTICES(dim);
627 } else {
628 mesh->n_node_el = 0;
629 }
630
631 if(dim > 1) {
632 mesh->node[EDGE] = mesh->n_node_el;
633 if (mesh->n_dof[EDGE] > 0) mesh->n_node_el += N_EDGES(dim);
634 }
635
636 if (dim == 3) {
637 mesh->node[FACE] = mesh->n_node_el;
638 if (mesh->n_dof[FACE] > 0) mesh->n_node_el += N_FACES_3D;
639 }
640
641 mesh->node[CENTER] = mesh->n_node_el;
642 if (mesh->n_dof[CENTER] > 0) mesh->n_node_el += 1;
643 #endif
644 }
645
646 /*--------------------------------------------------------------------------*/
647
AI_get_dof_admin(MESH * mesh,const char * name,const int n_dof[N_NODE_TYPES])648 DOF_ADMIN *AI_get_dof_admin(MESH *mesh, const char *name,
649 const int n_dof[N_NODE_TYPES])
650 {
651 FUNCNAME("AI_get_dof_admin");
652 DOF_ADMIN *admin;
653 int i;
654 DOF_ADMIN_MEM_INFO *mem_info;
655
656 admin = MEM_CALLOC(1, DOF_ADMIN);
657 admin->mesh = mesh;
658 admin->name = name ? strdup(name) : NULL;
659
660 admin->dof_free = NULL;
661 admin->dof_free_size = admin->first_hole = 0;
662
663 TEST_EXIT((mesh->dim > 1) || (n_dof[EDGE]==0),
664 "EDGE DOFs only make sense for mesh->dim > 1!\n");
665
666 TEST_EXIT((mesh->dim == 3) || (n_dof[FACE]==0),
667 "FACE DOFs only make sense for mesh->dim == 3!\n");
668
669 for (i = 0; i < N_NODE_TYPES; i++) admin->n_dof[i] = n_dof[i];
670
671 mem_info = (DOF_ADMIN_MEM_INFO *)
672 (admin->mem_info = MEM_ALLOC(1, DOF_ADMIN_MEM_INFO));
673
674 mem_info->dof_matrix = newObject(sizeof(DOF_MATRIX), 10, "dof_matrix");
675 mem_info->real_matrix_row =
676 newObject(sizeof(MATRIX_ROW_REAL), 0, "real_matrix_row");
677 mem_info->real_d_matrix_row =
678 newObject(sizeof(MATRIX_ROW_REAL_D), 0, "real_d_matrix_row");
679 mem_info->real_dd_matrix_row =
680 newObject(sizeof(MATRIX_ROW_REAL_DD), 0, "real_dd_matrix_row");
681
682 mem_info->dof_int_vec = newObject(sizeof(DOF_INT_VEC), 10,
683 "dof_int_vec");
684 mem_info->dof_dof_vec = newObject(sizeof(DOF_DOF_VEC), 10,
685 "dof_dof_vec");
686 mem_info->int_dof_vec = newObject(sizeof(DOF_DOF_VEC), 10,
687 "int_dof_vec");
688 mem_info->dof_uchar_vec = newObject(sizeof(DOF_UCHAR_VEC), 10,
689 "dof_uchar_vec");
690 mem_info->dof_schar_vec = newObject(sizeof(DOF_SCHAR_VEC), 10,
691 "dof_schar_vec");
692 mem_info->dof_real_vec = newObject(sizeof(DOF_REAL_VEC), 10,
693 "dof_real_vec");
694 mem_info->dof_real_d_vec = newObject(sizeof(DOF_REAL_D_VEC), 10,
695 "dof_real_d_vec");
696 mem_info->dof_real_dd_vec = newObject(sizeof(DOF_REAL_DD_VEC), 10,
697 "dof_real_dd_vec");
698 mem_info->dof_ptr_vec = newObject(sizeof(DOF_PTR_VEC), 10,
699 "dof_ptr_vec");
700
701 admin->compress_hooks.next =
702 admin->compress_hooks.prev =
703 &admin->compress_hooks;
704
705 add_dof_admin_to_mesh(admin, mesh);
706
707 return admin;
708 }
709
710
711 /*--------------------------------------------------------------------------*/
712 /* Free all DOF_ADMINs in a rather brutal way, only makes sense when */
713 /* freeing an entire mesh. */
714 /*--------------------------------------------------------------------------*/
715
free_dof_admins(MESH * mesh)716 static void free_dof_admins(MESH *mesh)
717 {
718 FUNCNAME("free_dof_admins");
719
720 DOF_ADMIN **admin = mesh->dof_admin;
721 DOF_ADMIN_MEM_INFO *admmeminfo;
722 int i, n;
723 DOF_MATRIX *dof_matrix, *dof_matrix_next;
724 DOF_INT_VEC *dof_int_vec, *dof_int_vec_next;
725 DOF_DOF_VEC *dof_dof_vec, *dof_dof_vec_next;
726 DOF_DOF_VEC *int_dof_vec, *int_dof_vec_next;
727 DOF_UCHAR_VEC *dof_uchar_vec, *dof_uchar_vec_next;
728 DOF_SCHAR_VEC *dof_schar_vec, *dof_schar_vec_next;
729 DOF_REAL_VEC *dof_real_vec, *dof_real_vec_next;
730 DOF_REAL_D_VEC *dof_real_d_vec, *dof_real_d_vec_next;
731 DOF_REAL_DD_VEC *dof_real_dd_vec, *dof_real_dd_vec_next;
732 DOF_PTR_VEC *dof_ptr_vec, *dof_ptr_vec_next;
733
734 n = mesh->n_dof_admin;
735 if ((n > 0) && !admin)
736 ERROR_EXIT("no mesh->dof_admin but n_dof_admin=%d\n", n);
737 if ((n <= 0) && admin)
738 ERROR_EXIT("found mesh->dof_admin but n_dof_admin=%d\n", n);
739
740 for(i = 0; i < n; i++) {
741 #define FREE_DOF_OBJ(admin, obj) \
742 for(obj = admin[i]->obj; obj; obj = obj##_next) { \
743 obj##_next = obj->next; \
744 free_##obj(obj); \
745 }
746 FREE_DOF_OBJ(admin, dof_matrix);
747 FREE_DOF_OBJ(admin, dof_int_vec);
748 FREE_DOF_OBJ(admin, dof_dof_vec);
749 FREE_DOF_OBJ(admin, int_dof_vec);
750 FREE_DOF_OBJ(admin, dof_uchar_vec);
751 FREE_DOF_OBJ(admin, dof_schar_vec);
752 FREE_DOF_OBJ(admin, dof_real_vec);
753 FREE_DOF_OBJ(admin, dof_real_d_vec);
754 FREE_DOF_OBJ(admin, dof_real_dd_vec);
755 FREE_DOF_OBJ(admin, dof_ptr_vec);
756
757 admmeminfo = (DOF_ADMIN_MEM_INFO *)admin[i]->mem_info;
758
759 deleteObject(admmeminfo->dof_matrix);
760 deleteObject(admmeminfo->real_matrix_row);
761 deleteObject(admmeminfo->real_d_matrix_row);
762 deleteObject(admmeminfo->real_dd_matrix_row);
763 deleteObject(admmeminfo->dof_int_vec);
764 deleteObject(admmeminfo->dof_dof_vec);
765 deleteObject(admmeminfo->int_dof_vec);
766 deleteObject(admmeminfo->dof_uchar_vec);
767 deleteObject(admmeminfo->dof_schar_vec);
768 deleteObject(admmeminfo->dof_real_vec);
769 deleteObject(admmeminfo->dof_real_d_vec);
770 deleteObject(admmeminfo->dof_real_dd_vec);
771 deleteObject(admmeminfo->dof_ptr_vec);
772
773 MEM_FREE(admin[i]->mem_info, 1, DOF_ADMIN_MEM_INFO);
774 MEM_FREE(admin[i]->dof_free, admin[i]->dof_free_size, DOF_FREE_UNIT);
775 }
776
777 return;
778 }
779
780
781 /*--------------------------------------------------------------------------*/
782 /* allocate memory for a new mesh, and build a simply linked list of these */
783 /*--------------------------------------------------------------------------*/
784
AI_advance_cookies_rec(MESH * mesh)785 void AI_advance_cookies_rec(MESH *mesh)
786 {
787 FUNCNAME("AI_advance_cookies_rec");
788 int i;
789 MESH_MEM_INFO *mem_info;
790
791 TEST_EXIT(mesh,"Oops, did not get a mesh!\n");
792 mem_info = (MESH_MEM_INFO *)mesh->mem_info;
793
794 mesh->cookie += 1;
795
796 for(i = 0; i < mem_info->n_slaves; i++)
797 AI_advance_cookies_rec(mem_info->slaves[i]);
798
799 return;
800 }
801
802 MESH *
_AI_get_mesh(int dim,const char * name,const MACRO_DATA * macro_data,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),AFF_TRAFO * (* init_wall_trafos)(MESH *,MACRO_EL *,int wall),bool strict_periodic)803 _AI_get_mesh(int dim, const char *name,
804 const MACRO_DATA *macro_data,
805 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
806 AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall),
807 bool strict_periodic)
808 {
809 FUNCNAME("get_mesh");
810 MESH *mesh;
811 MESH_MEM_INFO *mem_info;
812
813 mesh = MEM_CALLOC(1,MESH);
814
815 mesh->dim = dim;
816 mesh->name = name ? strdup(name) : NULL;
817
818 mem_info = (MESH_MEM_INFO *)
819 (mesh->mem_info = MEM_CALLOC(1, MESH_MEM_INFO));
820
821 mem_info->element = newObject(sizeof(EL), 0, "element");
822
823 if(mesh->dim == 3)
824 mem_info->rc_list = NULL;
825
826 mem_info->real_d = _newObject(sizeof(REAL_D), sizeof(REAL), 0, "real_d");
827
828 mem_info->leaf_data = NULL;
829
830 mem_info->next_trace_id = 0;
831
832 /* Flags them as invalid, originally */
833 mesh->n_vertices =
834 mesh->n_edges =
835 mesh->n_faces =
836 mesh->per_n_vertices =
837 mesh->per_n_edges =
838 mesh->per_n_faces = -1;
839
840 if (macro_data) {
841 _AI_macro_data2mesh(mesh, macro_data,
842 init_node_proj, init_wall_trafos, strict_periodic);
843 }
844
845 #if ALBERTA_DEBUG
846 srand(13);
847 #else
848 srand((unsigned int) time(NULL));
849 #endif
850
851 mesh->cookie = rand();
852
853 mesh->trace_id = -1;
854
855 check_mesh(mesh);
856
857 return mesh;
858 }
859
free_mesh(MESH * mesh)860 void free_mesh(MESH *mesh)
861 {
862 FUNCNAME("free_mesh");
863 MESH_MEM_INFO *mem_info;
864 int i;
865
866 if (!mesh) {
867 ERROR("No mesh specified!\n");
868 return;
869 }
870
871 mem_info = (MESH_MEM_INFO *)mesh->mem_info;
872
873 /* If this mesh is itself a slave mesh, then unchain it first. */
874 if(mem_info->master) unchain_submesh(mesh);
875 /* Free all slave meshes... */
876
877 for (i = 0; i < mem_info->n_slaves; i++) {
878 unchain_submesh(mem_info->slaves[i]);
879 }
880
881 /* free all elements/dofs/... */
882 if (mem_info->dof_ptrs) {
883 deleteObject(mem_info->dof_ptrs);
884 }
885
886 for (i = 0; i < N_NODE_TYPES; i++) {
887 if (mem_info->dofs[i]) {
888 deleteObject(mem_info->dofs[i]);
889 }
890 }
891
892 deleteObject(mem_info->element);
893 if (mem_info->rc_list)
894 free_rc_list(mesh, (RC_LIST_EL *)mem_info->rc_list);
895
896 deleteObject(mem_info->real_d);
897
898 if (mem_info->leaf_data) {
899 deleteObject(mem_info->leaf_data);
900 }
901
902 AI_free_dof_vec_list(mesh);
903 if (mesh->is_periodic) {
904 AI_free_dof_vec_list_np(mesh);
905 }
906
907 MEM_FREE(mem_info->coords, mem_info->count, REAL_D);
908
909 MEM_FREE(mem_info, 1, MESH_MEM_INFO);
910
911 MEM_FREE(mesh->macro_els, mesh->n_macro_el, MACRO_EL);
912
913 free_dof_admins(mesh);
914
915 MEM_FREE(mesh->dof_admin, mesh->n_dof_admin, DOF_ADMIN);
916
917 if(mesh->name) free((char *)mesh->name);
918
919 if (mesh->is_periodic) {
920 if (mesh->n_wall_trafos) {
921 MEM_FREE(mesh->wall_trafos,
922 mesh->n_wall_trafos*(sizeof(AFF_TRAFO *)+sizeof(AFF_TRAFO)),
923 char);
924 }
925 }
926
927 MEM_FREE(mesh, 1, MESH);
928
929 return;
930 }
931
932
933 MESH *
check_and_get_mesh(int dim,int dow,int debug,const char * version,const char * name,const MACRO_DATA * macro_data,NODE_PROJECTION * (* init_node_proj)(MESH *,MACRO_EL *,int),AFF_TRAFO * (* init_wall_trafos)(MESH *,MACRO_EL *,int wall))934 check_and_get_mesh(int dim, int dow, int debug,
935 const char *version, const char *name,
936 const MACRO_DATA *macro_data,
937 NODE_PROJECTION *(*init_node_proj)(MESH *, MACRO_EL *, int),
938 AFF_TRAFO *(*init_wall_trafos)(MESH *, MACRO_EL *, int wall))
939 {
940 FUNCNAME("check_and_get_mesh");
941 int error = 0;
942
943 if (dow != DIM_OF_WORLD)
944 {
945 ERROR("%s = %d, but you are using a lib with %s = %d\n",
946 "DIM_OF_WORLD", dow, "DIM_OF_WORLD", DIM_OF_WORLD);
947 error++;
948 }
949 if (dim > DIM_MAX)
950 {
951 ERROR("dim == %d > %d == DIM_MAX!\n", dim, DIM_MAX);
952 error++;
953 }
954 if (debug != ALBERTA_DEBUG)
955 {
956 ERROR("%s = %d, but you are using a lib with %s = %d\n",
957 "DEBUG", debug, "DEBUG", ALBERTA_DEBUG);
958 error++;
959 }
960 if (strcmp(version,ALBERTA_VERSION))
961 {
962 ERROR("you are using %s but a lib with %s\n", version, ALBERTA_VERSION);
963 error++;
964 }
965 if (error) ERROR_EXIT("Bye!\n");
966
967 return _AI_get_mesh(dim, name, macro_data,
968 init_node_proj, init_wall_trafos, false);
969 }
970
971
972 /*--------------------------------------------------------------------------*/
973 /* memory management for DOF pointers */
974 /*--------------------------------------------------------------------------*/
975
976 #if ALBERTA_DEBUG
977 static const int max_dof_ptrs[4] = {1, 3, 7, 15};
978 #endif
979
get_dof_ptrs(MESH * mesh)980 static DOF **get_dof_ptrs(MESH *mesh)
981 {
982 FUNCNAME("get_dof_ptrs");
983 int i, n;
984 DOF **ptrs;
985 MESH_MEM_INFO *mem_info;
986
987 DEBUG_TEST_EXIT(mesh, "mesh=NULL\n");
988 DEBUG_TEST_EXIT(mesh->mem_info, "mesh \"%s\": mesh->mem_info=NULL\n",
989 mesh->name);
990 mem_info = (MESH_MEM_INFO*)(mesh->mem_info);
991 n = mesh->n_node_el;
992 if (n <= 0)
993 return(NULL);
994
995 DEBUG_TEST_EXIT(n <= max_dof_ptrs[mesh->dim],
996 "mesh \"%s\": too many nodes: %d > %d\n",
997 mesh->name, n, max_dof_ptrs[mesh->dim]);
998 DEBUG_TEST_EXIT(mem_info->dof_ptrs,
999 "mesh \"%s\": mesh->mem_info->dof_ptrs=NULL\n", mesh->name);
1000
1001 ptrs = (DOF **)getMemory(mem_info->dof_ptrs);
1002 for (i = 0; i < n; i++)
1003 ptrs[i] = NULL;
1004
1005 return(ptrs);
1006 }
1007
1008
free_dof_ptrs(DOF ** ptrs,MESH * mesh)1009 static void free_dof_ptrs(DOF **ptrs, MESH *mesh)
1010 {
1011 FUNCNAME("free_dof_ptrs");
1012 int n;
1013 MESH_MEM_INFO *mem_info;
1014
1015 DEBUG_TEST_EXIT(ptrs,"ptrs=NULL\n");
1016 DEBUG_TEST_EXIT(mesh, "mesh=NULL\n");
1017 DEBUG_TEST_EXIT(mesh->mem_info,
1018 "mesh \"%s\": mesh->mem_info=NULL\n", mesh->name);
1019 n = mesh->n_node_el;
1020 if (n <= 0)
1021 return;
1022
1023 DEBUG_TEST_EXIT(n <= max_dof_ptrs[mesh->dim],
1024 "mesh \"%s\": too many nodes: %d > %d\n",
1025 mesh->name, n, max_dof_ptrs[mesh->dim]);
1026
1027 mem_info = (MESH_MEM_INFO*)(mesh->mem_info);
1028
1029 DEBUG_TEST_EXIT(mem_info->dof_ptrs,
1030 "mesh \"%s\": mesh->mem_info->dof_ptrs=NULL\n", mesh->name);
1031
1032 freeMemory((void *)ptrs, mem_info->dof_ptrs);
1033 return;
1034 }
1035
AI_get_dof_ptr_list(MESH * mesh)1036 void AI_get_dof_ptr_list(MESH *mesh)
1037 {
1038 FUNCNAME("AI_get_dof_ptr_list");
1039 MESH_MEM_INFO *mem_info;
1040
1041 DEBUG_TEST_EXIT(mesh, "No mesh given!\n");
1042
1043 if (mesh->n_node_el == 0) {
1044 return;
1045 }
1046
1047 mem_info = (MESH_MEM_INFO*)mesh->mem_info;
1048 DEBUG_TEST_EXIT(mem_info, "No mesh memory info structure present!\n");
1049
1050 mem_info->dof_ptrs = newObject(mesh->n_node_el * sizeof(DOF*), 1000,
1051 "dof_ptrs");
1052 }
1053
1054 /*--------------------------------------------------------------------------*/
1055 /* memory management for DOFs */
1056 /*--------------------------------------------------------------------------*/
1057
1058 #define DOF_BLOCK 1000
1059
1060 /****************************************************************************/
1061 /* AI_reactivate_dof(mesh, el): */
1062 /* When coarsening the mesh, we must replace all -1 settings with new DOF */
1063 /* indices. We no longer allocate new memory for DOF pointers in the */
1064 /* el->dof[] vector - this memory was also freed in old versions depending */
1065 /* on the "mesh->preserve_coarse_dofs" setting. */
1066 /****************************************************************************/
1067
AI_reactivate_dof(MESH * mesh,const EL * el,DOF ** edge_twins,DOF ** face_twins)1068 void AI_reactivate_dof(MESH *mesh, const EL *el,
1069 DOF **edge_twins, DOF **face_twins)
1070 {
1071 FUNCNAME("AI_reactivate_dof");
1072 DOF_ADMIN *admin;
1073 int i, j, n, n0, node;
1074
1075 DEBUG_TEST_EXIT(mesh,"mesh=NULL\n");
1076 DEBUG_TEST_EXIT(el, "el=NULL\n");
1077
1078 for (i = 0; i < mesh->n_dof_admin; i++) {
1079 admin = mesh->dof_admin[i];
1080 DEBUG_TEST_EXIT(admin, "mesh \"%s\": no dof_admin[%d]\n", mesh->name, i);
1081
1082 if(mesh->n_dof[CENTER]) {
1083 node = mesh->node[CENTER];
1084 n = admin->n_dof[CENTER];
1085
1086 if(n) {
1087 n0 = admin->n0_dof[CENTER];
1088
1089 DEBUG_TEST_EXIT(n+n0 <= mesh->n_dof[CENTER],
1090 "dof_admin \"%s\": n=%d, n0=%d too large: ndof[CENTER]=%d\n",
1091 admin->name, n, n0, mesh->n_dof[CENTER]);
1092
1093 if(el->dof[node][n0] == -1)
1094 for (j = 0; j < n; j++)
1095 el->dof[node][n0+j] = get_dof_index(admin);
1096 }
1097 }
1098
1099 #if DIM_MAX > 1
1100 if(mesh->n_dof[EDGE]) {
1101 int dim = mesh->dim;
1102 int k;
1103
1104 for(k = 0; k < N_EDGES(dim); k++) {
1105 node = mesh->node[EDGE] + k;
1106 n = admin->n_dof[EDGE];
1107
1108 if(n) {
1109 n0 = admin->n0_dof[EDGE];
1110
1111 DEBUG_TEST_EXIT(n+n0 <= mesh->n_dof[EDGE],
1112 "dof_admin \"%s\": n=%d, n0=%d too large: ndof[EDGE]=%d\n",
1113 admin->name, n, n0, mesh->n_dof[EDGE]);
1114
1115 if(el->dof[node][n0] == -1) {
1116 if ((admin->flags & ADM_PERIODIC) &&
1117 edge_twins && edge_twins[k] && edge_twins[k][n0] != -1) {
1118 for (j = 0; j < n; j++)
1119 el->dof[node][n0+j] = edge_twins[k][n0+j];
1120 } else {
1121 for (j = 0; j < n; j++)
1122 el->dof[node][n0+j] = get_dof_index(admin);
1123 }
1124 }
1125 }
1126 }
1127 }
1128
1129 #if DIM_MAX > 2
1130 if(mesh->n_dof[FACE]) {
1131 int k;
1132 for(k = 0; k < N_FACES_3D; k++) {
1133 node = mesh->node[FACE] + k;
1134 n = admin->n_dof[FACE];
1135
1136 if(n) {
1137 n0 = admin->n0_dof[FACE];
1138
1139 DEBUG_TEST_EXIT(n+n0 <= mesh->n_dof[FACE],
1140 "dof_admin \"%s\": n=%d, n0=%d too large: ndof[FACE]=%d\n",
1141 admin->name, n, n0, mesh->n_dof[FACE]);
1142
1143 if(el->dof[node][n0] == -1) {
1144 if ((admin->flags & ADM_PERIODIC) &&
1145 face_twins && face_twins[k] && face_twins[k][n0] != -1) {
1146 for (j = 0; j < n; j++)
1147 el->dof[node][n0+j] = face_twins[k][n0+j];
1148 } else {
1149 for (j = 0; j < n; j++)
1150 el->dof[node][n0+j] = get_dof_index(admin);
1151 }
1152 }
1153 }
1154 }
1155 }
1156 #endif
1157 #endif
1158 }
1159
1160 return;
1161 }
1162
AI_get_dof_memory(MESH * mesh,int position)1163 DOF *AI_get_dof_memory(MESH *mesh, int position)
1164 {
1165 FUNCNAME("AI_get_dof_memory");
1166 MESH_MEM_INFO *mem_info;
1167
1168
1169 DEBUG_TEST_EXIT(mesh, "mesh=NULL\n");
1170 mem_info = (MESH_MEM_INFO *)mesh->mem_info;
1171 DEBUG_TEST_EXIT(mem_info, "mesh \"%s\": mesh->mem_info=NULL\n", mesh->name);
1172
1173 DEBUG_TEST_EXIT(position >= 0 && position < N_NODE_TYPES,
1174 "mesh \"%s\": unknown position %d\n", mesh->name, position);
1175
1176 DEBUG_TEST_EXIT(mesh->n_dof[position], "mesh->n_dof[%d] == 0!\n", position);
1177
1178 return (DOF *)getMemory(mem_info->dofs[position]);
1179 }
1180
AI_free_dof_memory(DOF * dof,MESH * mesh,int position)1181 void AI_free_dof_memory(DOF *dof, MESH *mesh, int position)
1182 {
1183 FUNCNAME("AI_get_dof_memory");
1184 MESH_MEM_INFO *mem_info;
1185
1186
1187 DEBUG_TEST_EXIT(mesh, "mesh=NULL\n");
1188 mem_info = (MESH_MEM_INFO *)mesh->mem_info;
1189 DEBUG_TEST_EXIT(mem_info, "mesh \"%s\": mesh->mem_info=NULL\n", mesh->name);
1190
1191 DEBUG_TEST_EXIT(position >= 0 && position < N_NODE_TYPES,
1192 "mesh \"%s\": unknown position %d\n", mesh->name, position);
1193
1194 DEBUG_TEST_EXIT(mesh->n_dof[position], "mesh->n_dof[%d] == 0!\n", position);
1195
1196 freeMemory((void *)dof, mem_info->dofs[position]);
1197 }
1198
_AI_get_dof(MESH * mesh,int position,bool alloc_index)1199 DOF *_AI_get_dof(MESH *mesh, int position, bool alloc_index)
1200 {
1201 FUNCNAME("get_dof");
1202 DOF_ADMIN *admin;
1203 DOF *dof;
1204 int i, j, n, n0, ndof;
1205
1206 ndof = mesh->n_dof[position];
1207 if (ndof <= 0) return(NULL);
1208
1209 dof = AI_get_dof_memory(mesh, position);
1210
1211 for (i = 0; i < mesh->n_dof_admin; i++) {
1212 admin = mesh->dof_admin[i];
1213 DEBUG_TEST_EXIT(admin, "mesh \"%s\": no dof_admin[%d]\n", mesh->name, i);
1214
1215 n = admin->n_dof[position];
1216 n0 = admin->n0_dof[position];
1217
1218 DEBUG_TEST_EXIT(n+n0 <= ndof,
1219 "dof_admin \"%s\": n=%d, n0=%d too large: ndof=%d\n",
1220 admin->name, n, n0, ndof);
1221
1222 if (alloc_index) {
1223 for (j = 0; j < n; j++) {
1224 dof[n0+j] = get_dof_index(admin);
1225 }
1226 }
1227 }
1228
1229 return dof;
1230 }
1231
get_dof(MESH * mesh,int position)1232 DOF *get_dof(MESH *mesh, int position)
1233 {
1234 return _AI_get_dof(mesh, position, true);
1235 }
1236
1237 /* Get a "periodic copy" of "twin" where some DOF_ADMINs are periodic,
1238 * others aren't.
1239 */
get_periodic_dof(MESH * mesh,int position,const DOF * twin)1240 DOF *get_periodic_dof(MESH *mesh, int position, const DOF *twin)
1241 {
1242 FUNCNAME("get_dof");
1243 DOF_ADMIN *admin;
1244 DOF *dof;
1245 int i, j, n, n0, ndof;
1246
1247 ndof = mesh->n_dof[position];
1248 if (ndof <= 0) return(NULL);
1249
1250 dof = AI_get_dof_memory(mesh, position);
1251
1252 for (i = 0; i < mesh->n_dof_admin; i++) {
1253 admin = mesh->dof_admin[i];
1254 DEBUG_TEST_EXIT(admin, "mesh \"%s\": no dof_admin[%d]\n", mesh->name, i);
1255
1256 n = admin->n_dof[position];
1257 n0 = admin->n0_dof[position];
1258
1259 DEBUG_TEST_EXIT(n+n0 <= ndof,
1260 "dof_admin \"%s\": n=%d, n0=%d too large: ndof=%d\n",
1261 admin->name, n, n0, ndof);
1262
1263 if (twin && (admin->flags & ADM_PERIODIC)) {
1264 for (j = 0; j < n; j++)
1265 dof[n0+j] = twin[n0+j];
1266 } else {
1267 for (j = 0; j < n; j++)
1268 dof[n0+j] = get_dof_index(admin);
1269 }
1270 }
1271
1272 return(dof);
1273 }
1274
1275 /* "flags" can take the same values as DOF_ADMIN->flags:
1276 * ADM_PRESERVE_COARSE_DOFS: dof index is not freed if this admin
1277 * preserves coarse DOFs. The DOF memory is also not freed.
1278 * ADM_PERIODIC: dof index is not freed if the admin is periodic.
1279 */
free_dof(DOF * dof,MESH * mesh,int position,FLAGS flags)1280 void free_dof(DOF *dof, MESH *mesh, int position, FLAGS flags)
1281 {
1282 FUNCNAME("free_dof");
1283 DOF_ADMIN *admin;
1284 int i, j, n, n0, ndof;
1285 MESH_MEM_INFO *mem_info;
1286 FLAGS admflags;
1287
1288 DEBUG_TEST_EXIT(mesh, "mesh=NULL\n");
1289 mem_info = (MESH_MEM_INFO*)mesh->mem_info;
1290 DEBUG_TEST_EXIT(mem_info,"mesh \"%s\": mesh->mem_info=NULL\n", mesh->name);
1291
1292 DEBUG_TEST_EXIT(position >= 0 && position < N_NODE_TYPES,
1293 "mesh \"%s\": unknown position %d\n",mesh->name, position);
1294
1295 ndof = mesh->n_dof[position];
1296 (void)ndof;
1297
1298 DEBUG_TEST_EXIT(!ndof || dof,
1299 "dof = NULL, but ndof=%d\n", ndof);
1300 DEBUG_TEST_EXIT(ndof || !dof,
1301 "dof != NULL, but ndof=0\n");
1302
1303 DEBUG_TEST_EXIT(mem_info->dofs[position],
1304 "mesh \"%s\": no memory management present for %d DOFs.",
1305 mesh->name, position);
1306
1307 for (i = 0; i < mesh->n_dof_admin; i++) {
1308 admin = mesh->dof_admin[i];
1309 DEBUG_TEST_EXIT(admin, "mesh \"%s\": no dof_admin[%d]\n", mesh->name, i);
1310
1311 admflags = flags & admin->flags;
1312
1313 n0 = admin->n0_dof[position];
1314
1315 #if 0
1316 if (admin->stride) {
1317 n = admin->n_dof_stride[position];
1318 if ((admflags & ADM_PRESERVE_COARSE_DOFS) == 0) {
1319 for (j = 0; j < n; j++) {
1320 if (!(admflags & ADM_PERIODIC)) {
1321 free_dof_index(admin, dof[n0+j], true);
1322 }
1323 /* Mark this DOF as unused! For stride DOFs only one index
1324 * is stored in the mesh.
1325 */
1326 dof[n0+j] = -1;
1327 }
1328 }
1329 n0 += n;
1330 n = admin->n_dof[position] - n;
1331 } else
1332 #endif
1333 n = admin->n_dof[position];
1334
1335 DEBUG_TEST_EXIT(n+n0 <= ndof,
1336 "dof_admin \"%s\": n=%d, n0=%d too large: ndof=%d\n",
1337 admin, n, n0, ndof);
1338
1339 if ((admflags & ADM_PRESERVE_COARSE_DOFS) == 0) {
1340 for (j = 0; j < n; j++) {
1341 if (!(admflags & ADM_PERIODIC)) {
1342 free_dof_index(admin, dof[n0+j]);
1343 }
1344 dof[n0+j] = DOF_UNUSED; /* Mark this DOF as unused! */
1345 }
1346 }
1347 }
1348
1349 if ((flags & ADM_PRESERVE_COARSE_DOFS) == 0) {
1350 /* Also free the DOF memory. */
1351 freeMemory((void *)dof, mem_info->dofs[position]);
1352 }
1353
1354 return;
1355 }
1356
1357 /****************************************************************************/
1358 /* AI_get_dof_list(mesh, position): */
1359 /* Allocate a memory management object to administrate DOFs. The el->dof[] */
1360 /* entries will point into this memory space. Please note that we will also */
1361 /* allocate space for freed coarse DOFs - the allocated DOFs will have value*/
1362 /* -1 to show that they are no longer being used. */
1363 /* */
1364 /* This routine is also called from read_mesh.c. */
1365 /****************************************************************************/
1366
AI_get_dof_list(MESH * mesh,int position)1367 void AI_get_dof_list(MESH *mesh, int position)
1368 {
1369 FUNCNAME("AI_get_dof_list");
1370 MESH_MEM_INFO *mem_info;
1371
1372 DEBUG_TEST_EXIT(mesh, "No mesh given!\n");
1373 DEBUG_TEST_EXIT(position >=0 && position < N_NODE_TYPES,
1374 "Illegal position %d!\n", position);
1375 DEBUG_TEST_EXIT(mesh->n_dof[position],
1376 "Mesh has no DOFs on this position!\n");
1377
1378 mem_info = (MESH_MEM_INFO*)mesh->mem_info;
1379
1380 DEBUG_TEST_EXIT(mem_info, "No mesh memory info structure found!\n");
1381
1382 mem_info->dofs[position] =
1383 newObject(sizeof(DOF)*mesh->n_dof[position], DOF_BLOCK, "dof[pos]");
1384
1385 return;
1386 }
1387
1388
1389 /*--------------------------------------------------------------------------*/
1390 /* memory management for FE_SPACE and DOF_ADMIN structures */
1391 /*--------------------------------------------------------------------------*/
1392
1393 /* transfer_dofs(mesh, new_admin, old_dof, position,
1394 * is_coarse_dof, periodic_twin):
1395 *
1396 * We allocate and return memory for a new dof pointer in an el->dof[]
1397 * entry. The field is filled either with new DOF indices for each
1398 * admin or with -1 to mark it as unused. If "periodic_twin" != NULL
1399 * and (new_admin->flags & ADM_PERIODIC), then dof indices will not be
1400 * allocated but copied over from periodic_twin.
1401 */
transfer_dofs(MESH * mesh,DOF_ADMIN * new_admin,DOF * old_dof,int position,int is_coarse_dof,const DOF * periodic_twin)1402 static DOF *transfer_dofs(MESH *mesh, DOF_ADMIN *new_admin,
1403 DOF *old_dof, int position,
1404 int is_coarse_dof, const DOF *periodic_twin)
1405 {
1406 /* FUNCNAME("transfer_dofs"); */
1407 DOF_ADMIN *admin;
1408 DOF *new_dof = NULL;
1409 int i, j, n, n0, ndof = mesh->n_dof[position];
1410
1411 if (ndof <= 0) return NULL;
1412
1413 new_dof = AI_get_dof_memory(mesh, position);
1414
1415 for (i = 0; i < mesh->n_dof_admin; i++) {
1416 admin = mesh->dof_admin[i];
1417
1418 n = admin->n_dof[position];
1419 n0 = admin->n0_dof[position];
1420
1421 for (j = 0; j < n; j++) {
1422 if (admin == new_admin) {
1423 if ((admin->flags & ADM_PERIODIC) && periodic_twin) {
1424 new_dof[n0+j] = periodic_twin[n0+j];
1425 } else if (!is_coarse_dof ||
1426 (admin->flags & ADM_PRESERVE_COARSE_DOFS)) {
1427 new_dof[n0+j] = get_dof_index(admin);
1428 } else {
1429 new_dof[n0+j] = -1;
1430 }
1431 } else {
1432 if(old_dof) {
1433 new_dof[n0+j] = old_dof[n0 + j];
1434 } else {
1435 new_dof[n0+j] = -1;
1436 }
1437 }
1438 }
1439 }
1440
1441 return new_dof;
1442 }
1443
1444
1445 #include "memory_0d.c"
1446 #include "memory_1d.c"
1447 #if DIM_MAX > 1
1448 #include "memory_2d.c"
1449 #if DIM_MAX > 2
1450 #include "memory_3d.c"
1451 #endif
1452 #endif
1453
fe_space_ref(const FE_SPACE * fe_space)1454 static inline int fe_space_ref(const FE_SPACE *fe_space)
1455 {
1456 FE_SPACE *mutable = (FE_SPACE *)fe_space;
1457
1458 return ++mutable->ref_cnt;
1459 }
1460
fe_space_unref(const FE_SPACE * fe_space)1461 static inline int fe_space_unref(const FE_SPACE *fe_space)
1462 {
1463 FE_SPACE *mutable = (FE_SPACE *)fe_space;
1464
1465 return --mutable->ref_cnt;
1466 }
1467
copy_fe_space(const FE_SPACE * fe_space)1468 const FE_SPACE *copy_fe_space(const FE_SPACE *fe_space)
1469 {
1470 #if 0
1471 FE_SPACE *clone;
1472
1473 if (fe_space->bas_fcts == NULL) {
1474 clone = (FE_SPACE *)get_dof_space(fe_space->mesh,
1475 fe_space->name,
1476 fe_space->admin->n_dof,
1477 fe_space->admin->flags);
1478 } else if (fe_space->admin) {
1479 clone = (FE_SPACE *)get_fe_space(fe_space->mesh,
1480 fe_space->name,
1481 fe_space->bas_fcts,
1482 fe_space->rdim,
1483 fe_space->admin->flags);
1484 } else {
1485 clone = MEM_CALLOC(1, FE_SPACE);
1486 *clone = *fe_space;
1487 clone->name = fe_space->name ? strdup(fe_space->name) : NULL;
1488 CHAIN_INIT(clone);
1489 }
1490
1491 return clone;
1492 #else
1493 if (fe_space) {
1494 CHAIN_DO(fe_space, const FE_SPACE) {
1495 fe_space_ref(fe_space);
1496 fe_space_ref(fe_space->unchained);
1497 } CHAIN_WHILE(fe_space, const FE_SPACE);
1498 }
1499
1500 return fe_space;
1501 #endif
1502 }
1503
1504 /* Same as copy_fe_space(), but with a potentially different range dimension. */
clone_fe_space(const FE_SPACE * fe_space,int rdim)1505 const FE_SPACE *clone_fe_space(const FE_SPACE *fe_space, int rdim)
1506 {
1507 if (fe_space->bas_fcts == NULL || fe_space->rdim == rdim) {
1508 return copy_fe_space(fe_space);
1509 } else {
1510 return get_fe_space(fe_space->mesh,
1511 fe_space->name,
1512 fe_space->bas_fcts,
1513 rdim,
1514 fe_space->admin->flags);
1515 }
1516 }
1517
1518
1519 /*
1520 * get_fe_space(mesh, name, n_dof, bas_fcts, flags):
1521 *
1522 * Unlike in older versions of ALBERTA, this routine can be called at
1523 * any time, even after mesh refinement. As this is implemented at the
1524 * price of temporary memory usage and CPU time, it should still be
1525 * avoided. DK.
1526 *
1527 * Actually, extra temporary memory usage and computing power is only
1528 * need when the underlying DOF_ADMIN does not exist yet. cH.
1529 */
get_fe_space(MESH * mesh,const char * name,const BAS_FCTS * bas_fcts,int rdim,FLAGS adm_flags)1530 const FE_SPACE *get_fe_space(MESH *mesh,
1531 const char *name,
1532 const BAS_FCTS *bas_fcts,
1533 int rdim,
1534 FLAGS adm_flags)
1535 {
1536 FUNCNAME("get_fe_space");
1537 FE_SPACE *fe_space, *fe_space_chain, *fesp;
1538 BAS_FCTS *bfcts_chain;
1539 char fe_name[1024];
1540 const char *oname = name;
1541 int maxdim;
1542
1543 TEST_EXIT(bas_fcts->dim == mesh->dim,
1544 "Dimension of basis functions %d "
1545 "does not match mesh dimension %d!\n",
1546 bas_fcts->dim, mesh->dim);
1547
1548 if (oname && oname != bas_fcts->name) {
1549 snprintf(fe_name, sizeof(fe_name), "%s (@%s)", oname, bas_fcts->name);
1550 name = fe_name;
1551 } else {
1552 name = bas_fcts->name;
1553 }
1554 if (bas_fcts->trace_admin >= 0) {
1555 MESH *trace_mesh = lookup_submesh_by_id(mesh, bas_fcts->trace_admin);
1556 TEST_EXIT(trace_mesh != NULL,
1557 "Required trace-mesh with id %d not found.\n",
1558 bas_fcts->trace_admin);
1559 fe_space =
1560 (FE_SPACE *)get_dof_space(trace_mesh, name, bas_fcts->n_dof, adm_flags);
1561 fe_space->mesh = mesh; /* but the admin keeps the trace-mesh */
1562 } else {
1563 fe_space =
1564 (FE_SPACE *)get_dof_space(mesh, name, bas_fcts->n_dof, adm_flags);
1565 }
1566 fe_space->bas_fcts = bas_fcts;
1567 fe_space->rdim = rdim;
1568
1569 if (CHAIN_SINGLE(bas_fcts)) {
1570 fe_space->unchained = fe_space;
1571 } else {
1572 /* make a flat copy */
1573 fe_space->unchained = (fesp = MEM_ALLOC(1, FE_SPACE));
1574 --fe_space->ref_cnt; /* unchained has been deleted here */
1575 *fesp = *fe_space;
1576 CHAIN_INIT(fesp);
1577 fesp->bas_fcts = bas_fcts->unchained;
1578 if (fesp->name) {
1579 fesp->name = strdup(fesp->name);
1580 }
1581 }
1582
1583 /* If bas_fcts is a direct sum, then aquire an fe-space for each
1584 * component and chain the fe-spaces together.
1585 */
1586 maxdim = bas_fcts->rdim;
1587 CHAIN_FOREACH(bfcts_chain, bas_fcts, BAS_FCTS) {
1588 if (oname && oname != bas_fcts->name) {
1589 snprintf(fe_name, sizeof(fe_name), "%s (@%s)", oname, bfcts_chain->name);
1590 } else {
1591 name = bfcts_chain->name;
1592 }
1593 if (bfcts_chain->trace_admin >= 0) {
1594 MESH *trace_mesh = lookup_submesh_by_id(mesh, bfcts_chain->trace_admin);
1595 TEST_EXIT(trace_mesh != NULL,
1596 "Required trace-mesh with id %d not found.\n",
1597 bfcts_chain->trace_admin);
1598 fe_space_chain = (FE_SPACE *)get_dof_space(
1599 trace_mesh, name, bfcts_chain->n_dof, adm_flags);
1600 fe_space_chain->mesh = mesh; /* but the admin keeps the trace-mesh */
1601 } else {
1602 fe_space_chain =
1603 (FE_SPACE *)get_dof_space(mesh, name, bfcts_chain->n_dof, adm_flags);
1604 }
1605 fe_space_chain->bas_fcts = bfcts_chain;
1606 fe_space_chain->rdim = rdim;
1607 maxdim = MAX(bfcts_chain->rdim, maxdim);
1608
1609 fe_space_chain->unchained = (fesp = MEM_ALLOC(1, FE_SPACE));
1610 --fe_space_chain->ref_cnt; /* unchained has been deleted here */
1611 *fesp = *fe_space_chain;
1612 CHAIN_INIT(fesp);
1613 fesp->bas_fcts = bfcts_chain->unchained;
1614 if (fesp->name) {
1615 fesp->name = strdup(fesp->name);
1616 }
1617 CHAIN_ADD_TAIL(fe_space, fe_space_chain);
1618 }
1619 if (rdim != -1 && maxdim > rdim) {
1620 WARNING("%d dimensional range requested < "
1621 "range dimension %d of basis functions",
1622 rdim, maxdim);
1623 }
1624
1625 TEST_EXIT(VECTOR_BASIS_FUNCTIONS || maxdim == 1,
1626 "DIM_OF_WORLD-valued basis functions not supported in this "
1627 "version of ALBERTA. re-configure and re-compile the library "
1628 "without \"--disable-vector-basis-functions\".\n");
1629
1630 return fe_space;
1631 }
1632
get_dof_space(MESH * mesh,const char * name,const int n_dof[N_NODE_TYPES],FLAGS flags)1633 const FE_SPACE *get_dof_space(MESH *mesh,
1634 const char *name,
1635 const int n_dof[N_NODE_TYPES],
1636 FLAGS flags)
1637 {
1638 FUNCNAME("get_dof_space");
1639 DOF_ADMIN *admin = NULL;
1640 FE_SPACE *fe_space;
1641 int i, j, good_admin;
1642
1643 if (!mesh->is_periodic) {
1644 flags &= ~ADM_PERIODIC;
1645 }
1646
1647 fe_space = MEM_CALLOC(1, FE_SPACE);
1648 fe_space->name = name ? strdup(name) : NULL;
1649
1650 /****************************************************************************/
1651 /* Search for a fitting DOF_ADMIN to serve the required n_dof field and */
1652 /* preserve_coarse_dofs status. Both have to match exactly! */
1653 /****************************************************************************/
1654 for (i = 0; i < mesh->n_dof_admin; i++) {
1655 admin = mesh->dof_admin[i];
1656 good_admin = true;
1657
1658 for (j = 0; j < N_NODE_TYPES; j++) {
1659 if (admin->n_dof[j] != n_dof[j]) {
1660 good_admin = false;
1661 break;
1662 }
1663 }
1664 if (admin->flags != flags)
1665 good_admin = false;
1666
1667 if(good_admin == true) break;
1668
1669 admin = NULL;
1670 }
1671
1672 if (!admin) {
1673 /****************************************************************************/
1674 /* We did not find a fitting admin. In this case, we must adjust the */
1675 /* mem_info->dof_ptrs and mem_info->dofs[i] fields. */
1676 /****************************************************************************/
1677 int old_n_node_el;
1678 int old_n_dof[N_NODE_TYPES];
1679 int old_node[N_NODE_TYPES];
1680 void *old_dof_ptrs;
1681 void *old_dofs[N_NODE_TYPES];
1682 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
1683
1684 if (!mesh->n_dof[VERTEX] &&
1685 (!n_dof[VERTEX] ||
1686 (!(flags & ADM_PERIODIC) && mesh->is_periodic))) {
1687 /* The first DOF_ADMIN with vertex DOFs is used to determine
1688 * things like the relative orientation of elements etc. On
1689 * periodic meshes this should be a periodic DOF_ADMIN, so we
1690 * make sure that this is the case.
1691 *
1692 * get_vertex_admin() will re-enter get_fe_space() (this
1693 * function), but this is ok, it will not enter this
1694 * if-clause.
1695 */
1696 (void)get_vertex_admin(mesh, ADM_PERIODIC);
1697 }
1698
1699 old_n_node_el = mesh->n_node_el;
1700 old_dof_ptrs = mem_info->dof_ptrs;
1701 for(i = 0; i < N_NODE_TYPES; i++) {
1702 old_n_dof[i] = mesh->n_dof[i];
1703 old_node[i] = mesh->node[i];
1704 old_dofs[i] = mem_info->dofs[i];
1705 }
1706
1707 admin = AI_get_dof_admin(mesh, name, n_dof);
1708 admin->flags = flags;
1709
1710 /* We must now adjust the mem_info->dofs[i] memory management object. */
1711 for(i = 0; i < N_NODE_TYPES; i++)
1712 if(n_dof[i])
1713 AI_get_dof_list(mesh, i);
1714
1715 /* Did mesh->n_node_el change while adding the new admin? If so, we have to */
1716 /* adjust the mem_info->dof_ptrs memory management object. */
1717 if(old_n_node_el < mesh->n_node_el)
1718 AI_get_dof_ptr_list(mesh);
1719
1720 /* Now we magically adjust all pointers by traversing the mesh. :-) */
1721
1722 switch(mesh->dim) {
1723 case 0:
1724 adjust_dofs_and_dof_ptrs_0d(mesh, admin,
1725 old_n_node_el, old_n_dof, old_node);
1726 break;
1727 case 1:
1728 adjust_dofs_and_dof_ptrs_1d(mesh, admin,
1729 old_n_node_el, old_n_dof, old_node);
1730 break;
1731 #if DIM_MAX > 1
1732 case 2:
1733 adjust_dofs_and_dof_ptrs_2d(mesh, admin,
1734 old_n_node_el, old_n_dof, old_node);
1735 break;
1736 #if DIM_MAX >= 3
1737 case 3:
1738 adjust_dofs_and_dof_ptrs_3d(mesh, admin,
1739 old_n_node_el, old_n_dof, old_node);
1740 break;
1741 #endif
1742 #endif
1743 default:
1744 ERROR_EXIT("Illegal mesh dimension!\n");
1745 }
1746 /* Release old memory management objects. */
1747
1748 if(old_n_node_el < mesh->n_node_el && old_dof_ptrs)
1749 deleteObject(old_dof_ptrs);
1750
1751 for(i = 0; i < N_NODE_TYPES; i++)
1752 if(n_dof[i] && old_dofs[i])
1753 deleteObject(old_dofs[i]);
1754 }
1755
1756 fe_space->admin = admin;
1757 fe_space->bas_fcts = NULL;
1758 fe_space->mesh = mesh;
1759 fe_space->unchained = fe_space;
1760 fe_space->rdim = -1;
1761 fe_space->ref_cnt = 2; /* unchained counts twice */
1762
1763 CHAIN_INIT(fe_space);
1764
1765 return fe_space;
1766 }
1767
free_fe_space(const FE_SPACE * fe_space)1768 void free_fe_space(const FE_SPACE *fe_space)
1769 {
1770 FUNCNAME("free_fe_space");
1771 DBL_LIST_NODE *next;
1772 const FE_SPACE *chain;
1773 bool destruction = false;
1774 bool direct_sum = false;
1775
1776 (void)destruction;
1777 (void)direct_sum;
1778
1779 if (fe_space == NULL) {
1780 ERROR("No fe_space specified!\n");
1781 return;
1782 }
1783
1784 CHAIN_FOREACH_SAFE(chain, next, fe_space, const FE_SPACE) {
1785 direct_sum = true;
1786
1787 fe_space_unref(chain);
1788 fe_space_unref(chain->unchained);
1789
1790 DEBUG_TEST_EXIT(chain->ref_cnt >= 0 && chain->unchained->ref_cnt >= 0,
1791 "Negative reference counts.\n");
1792
1793 if (chain->unchained != chain && chain->unchained->ref_cnt == 0) {
1794 if (chain->unchained->name != NULL) {
1795 free((void *)chain->unchained->name);
1796 }
1797 MEM_FREE(chain->unchained, 1, FE_SPACE);
1798 }
1799
1800 if (chain->ref_cnt == 0) {
1801 if(chain->name) {
1802 free((char *)chain->name);
1803 }
1804 MEM_FREE(chain, 1, FE_SPACE);
1805 destruction = true;
1806 }
1807 }
1808
1809 fe_space_unref(fe_space);
1810 fe_space_unref(fe_space->unchained);
1811
1812 DEBUG_TEST_EXIT(fe_space->ref_cnt >= 0 && fe_space->unchained->ref_cnt >= 0,
1813 "Negative reference counts.\n");
1814
1815 if (fe_space->unchained != fe_space && fe_space->unchained->ref_cnt == 0) {
1816 if (fe_space->unchained->name != NULL) {
1817 free((void *)fe_space->unchained->name);
1818 }
1819 MEM_FREE(fe_space->unchained, 1, FE_SPACE);
1820 }
1821
1822 DEBUG_TEST_EXIT(!direct_sum || destruction == (fe_space->ref_cnt == 0),
1823 "Reference counts are inconsistent within different "
1824 "members of a direct sum.\n");
1825
1826 if (fe_space->ref_cnt == 0) {
1827 if(fe_space->name) {
1828 free((char *)fe_space->name);
1829 }
1830 MEM_FREE(fe_space, 1, FE_SPACE);
1831 }
1832 }
1833
1834
1835 /****************************************************************************/
1836 /* AI_fill_missing_dofs(mesh): */
1837 /* This routine allocates currently unused element DOF pointers since this */
1838 /* is not done during read_mesh(). These corresponding values pointed to */
1839 /* will be -1. The reason: write_mesh() does not save these pointers in */
1840 /* the file, since this information is unnecessary to recreate DOF_ADMINs. */
1841 /* */
1842 /* called by read_mesh(). */
1843 /****************************************************************************/
1844
AI_fill_missing_dofs(MESH * mesh)1845 void AI_fill_missing_dofs(MESH *mesh)
1846 {
1847 FUNCNAME("AI_fill_missing_dofs");
1848
1849 DEBUG_TEST_EXIT(mesh, "Did not supply a mesh!\n");
1850
1851 switch(mesh->dim) {
1852 case 0:
1853 break;
1854 case 1:
1855 fill_missing_dofs_1d(mesh);
1856 break;
1857 #if DIM_MAX > 1
1858 case 2:
1859 fill_missing_dofs_2d(mesh);
1860 break;
1861 #if DIM_MAX >= 3
1862 case 3:
1863 fill_missing_dofs_3d(mesh);
1864 break;
1865 #endif
1866 #endif
1867 default:
1868 ERROR_EXIT("Illegal mesh dimension!\n");
1869 }
1870
1871 return;
1872 }
1873
1874
1875 /*--------------------------------------------------------------------------*/
1876 /* memory management for leaf_data */
1877 /*--------------------------------------------------------------------------*/
1878
1879 #define LEAF_DATA_BLOCK 1000
1880
AI_get_leaf_data(MESH * mesh)1881 void *AI_get_leaf_data(MESH *mesh)
1882 {
1883 FUNCNAME("AI_get_leaf_data");
1884 MEMORYADMIN *ldbi;
1885
1886 DEBUG_TEST_EXIT(mesh, "pointer to mesh = NULL\n");
1887
1888 ldbi = (MEMORYADMIN *)((MESH_MEM_INFO*)(mesh->mem_info))->leaf_data;
1889
1890 if (ldbi)
1891 return(getMemory(ldbi));
1892 else
1893 return(NULL);
1894 }
1895
AI_free_leaf_data(void * leaf_data,MESH * mesh)1896 void AI_free_leaf_data(void *leaf_data, MESH *mesh)
1897 {
1898 FUNCNAME("AI_free_leaf_data");
1899 MEMORYADMIN *ldbi;
1900
1901 if (!leaf_data) return;
1902
1903 DEBUG_TEST_EXIT(mesh, "pointer to mesh = NULL\n");
1904 ldbi = (MEMORYADMIN *) ((MESH_MEM_INFO*)(mesh->mem_info))->leaf_data;
1905
1906 if (!ldbi) return;
1907
1908 freeMemory((void *)leaf_data, ldbi);
1909 return;
1910 }
1911
1912
1913 /****************************************************************************/
1914 /* init_leaf_data(mesh, name, size, rld, cld): */
1915 /* Initialize leaf data on the mesh. */
1916 /****************************************************************************/
1917
init_leaf_data(MESH * mesh,size_t size,void (* refine_leaf_data)(EL * parent,EL * child[2]),void (* coarsen_leaf_data)(EL * parent,EL * child[2]))1918 size_t init_leaf_data(MESH *mesh, size_t size,
1919 void (*refine_leaf_data)(EL *parent, EL *child[2]),
1920 void (*coarsen_leaf_data)(EL *parent, EL *child[2]))
1921 {
1922 MESH_MEM_INFO *mem_info;
1923 size_t new_size;
1924 TRAVERSE_STACK *stack = get_traverse_stack();
1925 const EL_INFO *el_info = NULL;
1926
1927 TEST_EXIT(mesh, "No mesh specified!\n");
1928 TEST_EXIT(size, "size must be > 0!\n");
1929 TEST_EXIT(mesh->mem_info, "No memory management present for mesh!\n");
1930
1931 mem_info = (MESH_MEM_INFO *)mesh->mem_info;
1932
1933 TEST_EXIT(!mem_info->leaf_data, "Leaf data was already initialized!\n");
1934
1935 new_size = ALIGNMEM(size, sizeof(void *));
1936
1937 if (new_size != size)
1938 WARNING("installing leafdata of size %d with aligned size %d\n",
1939 size, new_size);
1940
1941 mem_info->leaf_data_info->leaf_data_size = new_size;
1942 mem_info->leaf_data_info->refine_leaf_data = refine_leaf_data;
1943 mem_info->leaf_data_info->coarsen_leaf_data = coarsen_leaf_data;
1944
1945 mem_info->leaf_data = newObject(new_size, 0, "leaf_data");
1946
1947 el_info = traverse_first(stack, mesh, -1, CALL_LEAF_EL);
1948
1949 while (el_info) {
1950 el_info->el->child[1] = (EL *) AI_get_leaf_data(mesh);
1951
1952 el_info = traverse_next(stack, el_info);
1953 }
1954
1955 free_traverse_stack(stack);
1956
1957 return new_size;
1958 }
1959
1960
1961 #if 0 /* not used at the moment */
1962 /*--------------------------------------------------------------------------*/
1963 /* memory management for quadratures */
1964 /*--------------------------------------------------------------------------*/
1965
1966 #define QUAD_VEC_BLOCK 100
1967 #define MAX_QUAD_INFO 3
1968
1969 static MEM_BLOCK_INFO *quad_vec_block_infos;
1970
1971 REAL *get_quad_vec(const QUAD *quad)
1972 {
1973 FUNCNAME("get_quad_vec");
1974 REAL *quad_vec;
1975 MEM_BLOCK_INFO *qvbi;
1976 int qvbi_idx;
1977
1978 DEBUG_TEST_EXIT(dim >= 1 && dim <= 3, "dim = %d not allowed\n");
1979
1980 qvbi_idx = quad->n_points / 16;
1981
1982 qvbi = quad_vec_block_infos[dim];
1983
1984 if (!gvbi)
1985 {
1986 gvbi = MEM_ALLOC(1, MEM_BLOCK_INFO);
1987 gvbi->free = NULL;
1988 gvbi->info = NULL;
1989 gvbi->free_count = 0;
1990 gvbi->used_count = 0;
1991 gvbi->unit_size = get_max_no_quad_point(dim)*sizeof(REAL);
1992 gvbi->block_size = QUAD_VEC_BLOCK;
1993
1994 quad_vec_block_infos[dim] = gvbi;
1995 }
1996
1997 quad_vec = (REAL *) get_mem(gvbi);
1998
1999 return(quad_vec);
2000 }
2001
2002 void free_quad_vec(REAL *quad_vec, int dim)
2003 {
2004 FUNCNAME("free_quad_vec");
2005 int n;
2006 MEM_BLOCK_INFO *gvbi;
2007
2008 if (!quad_vec) return;
2009
2010 DEBUG_TEST_EXIT(dim >= 1 && dim <= 3, "dim = %d not allowed\n");
2011
2012 gvbi = quad_vec_block_infos[dim];
2013
2014 if (!gvbi) return;
2015
2016 free_mem((void *) quad_vec, gvbi);
2017 return;
2018 }
2019
2020 static MEM_BLOCK_INFO *quad_vec_d_block_infos[MAX_QUAD_INFO+1];
2021
2022 REAL_D *get_quad_vec_d(int dim)
2023 {
2024 FUNCNAME("get_quad_vec_d");
2025 int i, n;
2026 REAL_D *quad_vec_d;
2027 MEM_BLOCK_INFO *gvdbi;
2028
2029 DEBUG_TEST_EXIT(dim >= 1 && dim <= 3, "dim = %d not allowed\n");
2030
2031 gvdbi = quad_vec_d_block_infos[dim];
2032
2033 if (!gvdbi)
2034 {
2035 gvdbi = MEM_ALLOC(1, MEM_BLOCK_INFO);
2036 gvdbi->free = NULL;
2037 gvdbi->info = NULL;
2038 gvdbi->free_count = 0;
2039 gvdbi->used_count = 0;
2040 gvdbi->unit_size = get_max_no_quad_point(dim)*sizeof(REAL_D);
2041 gvdbi->block_size = QUAD_VEC_BLOCK;
2042
2043 quad_vec_d_block_infos[dim] = gvdbi;
2044 }
2045
2046 quad_vec_d = (REAL_D *) get_mem(gvdbi);
2047
2048 return(quad_vec_d);
2049 }
2050
2051 void free_quad_vec_d(REAL_D *quad_vec_d, int dim)
2052 {
2053 FUNCNAME("free_quad_vec_d");
2054 int n;
2055 MEM_BLOCK_INFO *gvdbi;
2056
2057 if (!quad_vec_d) return;
2058
2059 DEBUG_TEST_EXIT(dim >= 1 && dim <= 3, "dim = %d not allowed\n");
2060
2061 gvdbi = quad_vec_d_block_infos[dim];
2062
2063 if (!gvdbi) return;
2064
2065 free_mem((void *) quad_vec_d, gvdbi);
2066 return;
2067 }
2068 #endif
2069
2070 /*--------------------------------------------------------------------------*/
2071 /* memory management for elements */
2072 /*--------------------------------------------------------------------------*/
2073
2074 #define EL_BLOCK 1000
2075
get_element(MESH * mesh)2076 EL *get_element(MESH *mesh)
2077 {
2078 FUNCNAME("get_element");
2079 EL *el;
2080 #if ALBERTA_DEBUG
2081 static int el_index = 0;
2082 #endif
2083
2084 DEBUG_TEST_EXIT(mesh, "mesh == NULL\n");
2085 DEBUG_TEST_EXIT(mesh->mem_info,
2086 "mesh \"%s\": no memory management present.\n", mesh->name);
2087
2088 el = (EL *)getMemory(((MESH_MEM_INFO*)(mesh->mem_info))->element);
2089 el->child[0] = NULL;
2090 el->child[1] = (EL *) AI_get_leaf_data(mesh);
2091 el->dof = get_dof_ptrs(mesh);
2092 #if ALBERTA_DEBUG
2093 el->index = el_index++;
2094 #endif
2095 el->mark = 0;
2096 el->new_coord = NULL;
2097
2098 return el;
2099 }
2100
free_element(EL * el,MESH * mesh)2101 void free_element(EL *el, MESH *mesh)
2102 {
2103 free_dof_ptrs(el->dof, mesh);
2104 if(mesh->dim > 1 && el->new_coord) {
2105 free_real_d(mesh, el->new_coord);
2106 el->new_coord = NULL;
2107 }
2108
2109 if (el->child[1]) AI_free_leaf_data((void *) el->child[1], mesh);
2110
2111 freeMemory(el, ((MESH_MEM_INFO*)(mesh->mem_info))->element);
2112 }
2113
2114
2115 /*--------------------------------------------------------------------------*/
2116 /* some routines for allocation and deallocation of list for collecting */
2117 /* neighbours at the refinement edge in 3 dimensions */
2118 /*--------------------------------------------------------------------------*/
2119
2120 #define LIST_BLOCK 20
2121
2122 /* static MEM_BLOCK_INFO rc_list_block_info = {NULL, NULL, 0, 0, 0, LIST_BLOCK}; */
2123
get_rc_list(MESH * mesh)2124 RC_LIST_EL *get_rc_list(MESH *mesh)
2125 {
2126 FUNCNAME("get_rc_list");
2127 MESH_MEM_INFO *mem_info = (MESH_MEM_INFO *)mesh->mem_info;
2128 int n_elem =
2129 mesh->is_periodic ? 2*mesh->max_edge_neigh : mesh->max_edge_neigh;
2130
2131 if (!mem_info->rc_list) {
2132 mem_info->rc_list =
2133 newObject(n_elem * sizeof(RC_LIST_EL), LIST_BLOCK, "rc_list");
2134 } else {
2135 DEBUG_TEST_EXIT(
2136 ((MEMORYADMIN *)mem_info->rc_list)->objectSize
2137 >=
2138 n_elem * sizeof(RC_LIST_EL),
2139 "mesh \"%s\": mesh->max_edge_neigh changed\n", mesh->name);
2140 }
2141
2142 return (RC_LIST_EL *)getMemory(mem_info->rc_list);
2143 }
2144
free_rc_list(MESH * mesh,RC_LIST_EL * list)2145 void free_rc_list(MESH *mesh, RC_LIST_EL *list)
2146 {
2147 freeMemory((void *)list, ((MESH_MEM_INFO *)(mesh->mem_info))->rc_list);
2148 }
2149
2150 /*--------------------------------------------------------------------------*/
2151 /*--------------------------------------------------------------------------*/
2152
get_real_d(MESH * mesh)2153 REAL *get_real_d(MESH *mesh)
2154 {
2155 FUNCNAME("get_real_d");
2156
2157 DEBUG_TEST_EXIT(mesh, "mesh==NULL\n");
2158 return (REAL *)getMemory(((MESH_MEM_INFO*)(mesh->mem_info))->real_d);
2159 }
2160
free_real_d(MESH * mesh,REAL * ptr)2161 void free_real_d(MESH *mesh, REAL *ptr)
2162 {
2163 FUNCNAME("free_real_d");
2164
2165 DEBUG_TEST_EXIT(mesh, "mesh==NULL\n");
2166 freeMemory((void *)ptr, ((MESH_MEM_INFO*)(mesh->mem_info))->real_d);
2167 return;
2168 }
2169
2170 /*--------------------------------------------------------------------------*/
2171 /* memory management for matrix rows */
2172 /* matrix rows not connected to a FE_SPACE, e.g. in multigrid.c are */
2173 /* treated specially. Others are fixed to a FE_SPACE and thus to a unique */
2174 /* DOF_ADMIN and MESH. */
2175 /*--------------------------------------------------------------------------*/
2176
2177 static void *unconnected_real_rows = NULL;
2178
get_matrix_row_real(const FE_SPACE * fe_space)2179 static MATRIX_ROW_REAL *get_matrix_row_real(const FE_SPACE *fe_space)
2180 {
2181 MATRIX_ROW_REAL *row;
2182 void *matrix_row_real_info;
2183 int j;
2184
2185 if (!fe_space || !fe_space->admin) {
2186 if (!unconnected_real_rows) {
2187 unconnected_real_rows = newObject(sizeof(MATRIX_ROW_REAL), 100,
2188 "unconnected rows");
2189 }
2190 matrix_row_real_info = unconnected_real_rows;
2191 } else {
2192 matrix_row_real_info =
2193 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_matrix_row;
2194 }
2195 row = (MATRIX_ROW_REAL *)getMemory(matrix_row_real_info);
2196 row->type = MATENT_REAL;
2197 row->next = NULL;
2198 for (j = 0; j < ROW_LENGTH; j++) {
2199 row->col[j] = NO_MORE_ENTRIES;
2200 }
2201
2202 return row;
2203 }
2204
free_matrix_row_real(const FE_SPACE * fe_space,MATRIX_ROW_REAL * row)2205 static void free_matrix_row_real(const FE_SPACE *fe_space, MATRIX_ROW_REAL *row)
2206 {
2207 void *matrix_row_info = NULL;
2208
2209 if (!fe_space || !fe_space->admin) {
2210 matrix_row_info = unconnected_real_rows;
2211 } else {
2212 matrix_row_info =
2213 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_matrix_row;
2214 }
2215 freeMemory((void *)row, matrix_row_info);
2216 }
2217
ALBERTA_DEFUNUSED(static void delete_unconnected_real_rows (void))2218 ALBERTA_DEFUNUSED(static void delete_unconnected_real_rows(void))
2219 {
2220 if (unconnected_real_rows) {
2221 deleteObject(unconnected_real_rows);
2222 }
2223 unconnected_real_rows = NULL;
2224 }
2225
2226 static void *unconnected_real_d_rows = NULL;
2227
get_matrix_row_real_d(const FE_SPACE * fe_space)2228 static MATRIX_ROW_REAL_D *get_matrix_row_real_d(const FE_SPACE *fe_space)
2229 {
2230 MATRIX_ROW_REAL_D *row;
2231 void *matrix_row_real_d_info;
2232 int j;
2233
2234 if (!fe_space || !fe_space->admin) {
2235 if (!unconnected_real_d_rows) {
2236 unconnected_real_d_rows = newObject(sizeof(MATRIX_ROW_REAL_D), 100,
2237 "unconnected rows");
2238 }
2239 matrix_row_real_d_info = unconnected_real_d_rows;
2240 } else {
2241 matrix_row_real_d_info =
2242 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_d_matrix_row;
2243 }
2244 row = (MATRIX_ROW_REAL_D *)getMemory(matrix_row_real_d_info);
2245 row->type = MATENT_REAL_D;
2246 row->next = NULL;
2247 for (j=0; j<ROW_LENGTH; j++) {
2248 row->col[j] = NO_MORE_ENTRIES;
2249 }
2250
2251 return row;
2252 }
2253
2254 static void
free_matrix_row_real_d(const FE_SPACE * fe_space,MATRIX_ROW_REAL_D * row)2255 free_matrix_row_real_d(const FE_SPACE *fe_space, MATRIX_ROW_REAL_D *row)
2256 {
2257 void *matrix_row_info = NULL;
2258
2259 if (!fe_space || !fe_space->admin) {
2260 matrix_row_info = unconnected_real_d_rows;
2261 } else {
2262 matrix_row_info =
2263 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_d_matrix_row;
2264 }
2265 freeMemory((void *)row, matrix_row_info);
2266 }
2267
ALBERTA_DEFUNUSED(static void delete_unconnected_real_d_rows (void))2268 ALBERTA_DEFUNUSED(static void delete_unconnected_real_d_rows(void))
2269 {
2270 if (unconnected_real_d_rows) {
2271 deleteObject(unconnected_real_d_rows);
2272 }
2273 unconnected_real_d_rows = NULL;
2274 }
2275
2276 static void *unconnected_real_dd_rows = NULL;
2277
get_matrix_row_real_dd(const FE_SPACE * fe_space)2278 static MATRIX_ROW_REAL_DD *get_matrix_row_real_dd(const FE_SPACE *fe_space)
2279 {
2280 MATRIX_ROW_REAL_DD *row;
2281 void *matrix_row_real_dd_info;
2282 int j;
2283
2284 if (!fe_space || !fe_space->admin) {
2285 if (!unconnected_real_dd_rows) {
2286 unconnected_real_dd_rows = newObject(sizeof(MATRIX_ROW_REAL_DD), 100,
2287 "unconnected rows");
2288 }
2289 matrix_row_real_dd_info = unconnected_real_dd_rows;
2290 } else {
2291 matrix_row_real_dd_info =
2292 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_dd_matrix_row;
2293 }
2294 row = (MATRIX_ROW_REAL_DD *)getMemory(matrix_row_real_dd_info);
2295 row->type = MATENT_REAL_DD;
2296 row->next = NULL;
2297 for (j=0; j<ROW_LENGTH; j++) {
2298 row->col[j] = NO_MORE_ENTRIES;
2299 }
2300
2301 return row;
2302 }
2303
2304 static void
free_matrix_row_real_dd(const FE_SPACE * fe_space,MATRIX_ROW_REAL_DD * row)2305 free_matrix_row_real_dd(const FE_SPACE *fe_space, MATRIX_ROW_REAL_DD *row)
2306 {
2307 void *matrix_row_info = NULL;
2308
2309 if (!fe_space || !fe_space->admin) {
2310 matrix_row_info = unconnected_real_dd_rows;
2311 } else {
2312 matrix_row_info =
2313 ((DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info)->real_dd_matrix_row;
2314 }
2315 freeMemory((void *)row, matrix_row_info);
2316 }
2317
ALBERTA_DEFUNUSED(static void delete_unconnected_real_dd_rows (void))2318 ALBERTA_DEFUNUSED(static void delete_unconnected_real_dd_rows(void))
2319 {
2320 if (unconnected_real_dd_rows) {
2321 deleteObject(unconnected_real_dd_rows);
2322 }
2323 unconnected_real_dd_rows = NULL;
2324 }
2325
get_matrix_row(const FE_SPACE * fe_space,MATENT_TYPE type)2326 MATRIX_ROW *get_matrix_row(const FE_SPACE *fe_space, MATENT_TYPE type)
2327 {
2328 FUNCNAME("get_matrix_row");
2329 switch (type) {
2330 case MATENT_REAL:
2331 return (MATRIX_ROW *)get_matrix_row_real(fe_space);
2332 case MATENT_REAL_D:
2333 return (MATRIX_ROW *)get_matrix_row_real_d(fe_space);
2334 case MATENT_REAL_DD:
2335 return (MATRIX_ROW *)get_matrix_row_real_dd(fe_space);
2336 default:
2337 ERROR_EXIT("Unsupported MATENT_TYPE: %d\n", type);
2338 }
2339 return NULL;
2340 }
2341
free_matrix_row(const FE_SPACE * fe_space,MATRIX_ROW * row)2342 void free_matrix_row(const FE_SPACE *fe_space, MATRIX_ROW *row)
2343 {
2344 FUNCNAME("free_matrix_row");
2345 switch (row->type) {
2346 case MATENT_REAL:
2347 free_matrix_row_real(fe_space, (MATRIX_ROW_REAL *)row);
2348 break;
2349 case MATENT_REAL_D:
2350 free_matrix_row_real_d(fe_space, (MATRIX_ROW_REAL_D *)row);
2351 break;
2352 case MATENT_REAL_DD:
2353 free_matrix_row_real_dd(fe_space, (MATRIX_ROW_REAL_DD *)row);
2354 break;
2355 default:
2356 ERROR_EXIT("Unsupported MATENT_TYPE: %d\n", row->type);
2357 break;
2358 }
2359 }
2360
2361 /*--------------------------------------------------------------------------*/
2362 /* memory management for matrices */
2363 /* unconnected_matrices: see unconnected rows above... */
2364 /*--------------------------------------------------------------------------*/
2365
2366 static void *unconnected_matrices = NULL;
2367
_AI_get_dof_matrix(const char * name,const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space)2368 static inline DOF_MATRIX *_AI_get_dof_matrix(const char *name,
2369 const FE_SPACE *row_fe_space,
2370 const FE_SPACE *col_fe_space)
2371 {
2372 DOF_MATRIX *mat;
2373 void *matrix_info = NULL;
2374
2375 if (!row_fe_space || !row_fe_space->admin) {
2376 if (!unconnected_matrices) {
2377 unconnected_matrices = newObject(sizeof(DOF_MATRIX), 10,
2378 "unconnected matrices");
2379 }
2380 matrix_info = unconnected_matrices;
2381 } else {
2382 matrix_info =
2383 ((DOF_ADMIN_MEM_INFO *)row_fe_space->admin->mem_info)->dof_matrix;
2384 }
2385
2386 mat = (DOF_MATRIX *)getMemory(matrix_info);
2387
2388 /* This is not timing critical in general, so we waste some CPU time
2389 * here and simply zap the entire matrix structure.
2390 */
2391 memset(mat, 0, sizeof(*mat));
2392
2393 mat->next = NULL;
2394 mat->row_fe_space = row_fe_space;
2395 mat->col_fe_space = col_fe_space;
2396 mat->name = name ? strdup(name) : NULL;
2397 mat->matrix_row = NULL;
2398 mat->size = 0;
2399 mat->type = MATENT_NONE;
2400 mat->n_entries = 0;
2401 mat->is_diagonal = false;
2402 mat->diagonal.real = NULL;
2403 mat->diag_cols = NULL;
2404 mat->inv_diag.real = NULL;
2405 mat->unchained = NULL;
2406 ROW_CHAIN_INIT(mat);
2407 COL_CHAIN_INIT(mat);
2408
2409 mat->refine_interpol = mat->coarse_restrict = NULL;
2410
2411 mat->mem_info = matrix_info;
2412
2413 if (row_fe_space && row_fe_space->admin) {
2414 #if 0
2415 if (row_fe_space->admin->n_dof[CENTER] == 1 &&
2416 row_fe_space->admin->n_dof[VERTEX] == 0 &&
2417 row_fe_space->admin->n_dof[EDGE] == 0 &&
2418 row_fe_space->admin->n_dof[FACE] == 0 &&
2419 (col_fe_space == NULL ||
2420 (col_fe_space->admin->n_dof[CENTER] == 1 &&
2421 col_fe_space->admin->n_dof[VERTEX] == 0 &&
2422 col_fe_space->admin->n_dof[EDGE] == 0 &&
2423 col_fe_space->admin->n_dof[FACE] == 0))) {
2424 mat->is_diagonal = true; /* e.g. element bubbles */
2425 }
2426 #endif
2427 add_dof_matrix_to_admin(mat, (DOF_ADMIN *)row_fe_space->admin);
2428 }
2429
2430 return mat;
2431 }
2432
_AI_free_dof_matrix(DOF_MATRIX * mat)2433 static inline void _AI_free_dof_matrix(DOF_MATRIX *mat)
2434 {
2435 if (mat->row_fe_space && mat->row_fe_space->admin) {
2436 remove_dof_matrix_from_admin(mat);
2437 }
2438
2439 clear_dof_matrix(mat); /* free all matrix_rows */
2440
2441 if (mat->matrix_row) {
2442 MEM_FREE(mat->matrix_row, mat->size, MATRIX_ROW *);
2443 mat->matrix_row = NULL;
2444 }
2445 if (mat->diag_cols) {
2446 free_dof_int_vec(mat->diag_cols);
2447 }
2448 mat->size = 0;
2449
2450 if (mat->name) {
2451 free((void *)mat->name);
2452 }
2453
2454 if (mat->mem_info) {
2455 freeMemory((void *)mat, mat->mem_info);
2456 } else {
2457 memset(mat, 0, sizeof(*mat));
2458 }
2459 }
2460
2461 /* Allocate a (possibly: bock-matrix), according to the structure of
2462 * row_fe_space and col_fe_space
2463 */
get_dof_matrix(const char * name,const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space)2464 DOF_MATRIX *get_dof_matrix(const char *name,
2465 const FE_SPACE *row_fe_space,
2466 const FE_SPACE *col_fe_space)
2467 {
2468 DOF_MATRIX *mat, *row_chain, *col_chain;
2469 const FE_SPACE *row_fesp, *col_fesp;
2470
2471 if (col_fe_space == NULL) {
2472 col_fe_space = row_fe_space;
2473 }
2474 row_fe_space = copy_fe_space(row_fe_space);
2475 col_fe_space = copy_fe_space(col_fe_space);
2476
2477 mat = _AI_get_dof_matrix(name, row_fe_space, col_fe_space);
2478
2479 if (row_fe_space) {
2480 /* Must be a call from MG-code if row_fe_space == NULL */
2481
2482 /* generate the first row */
2483 col_fe_space = mat->col_fe_space;
2484 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
2485 row_chain = _AI_get_dof_matrix(name, row_fe_space, col_fesp);
2486 ROW_CHAIN_ADD_TAIL(mat, row_chain);
2487 }
2488 CHAIN_FOREACH(row_fesp, row_fe_space, const FE_SPACE) {
2489 col_chain = _AI_get_dof_matrix(name, row_fesp, col_fe_space);
2490 COL_CHAIN_ADD_TAIL(mat, col_chain);
2491 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
2492 row_chain = _AI_get_dof_matrix(name, row_fesp, col_fesp);
2493 ROW_CHAIN_ADD_TAIL(col_chain, row_chain);
2494 mat = ROW_CHAIN_NEXT(mat, DOF_MATRIX);
2495 COL_CHAIN_ADD_TAIL(mat, row_chain);
2496 }
2497 mat = ROW_CHAIN_NEXT(mat, DOF_MATRIX); /* roll-over to the list head. */
2498 }
2499 }
2500
2501 return mat;
2502 }
2503
free_dof_matrix(DOF_MATRIX * mat)2504 void free_dof_matrix(DOF_MATRIX *mat)
2505 {
2506 DBL_LIST_NODE *row_next, *col_next;
2507 DOF_MATRIX *row_chain, *col_chain;
2508
2509 if (mat->row_fe_space) {
2510 free_fe_space(mat->row_fe_space);
2511 free_fe_space(mat->col_fe_space);
2512 }
2513
2514 /* delete all columns, safe the first one */
2515 ROW_CHAIN_FOREACH_SAFE(row_chain, row_next, mat, DOF_MATRIX) {
2516 /* delete the column down to the pre-last entry */
2517 COL_CHAIN_FOREACH_SAFE(col_chain, col_next, row_chain, DOF_MATRIX) {
2518 ROW_CHAIN_DEL(col_chain);
2519 COL_CHAIN_DEL(col_chain);
2520 _AI_free_dof_matrix(col_chain);
2521 }
2522 /* delete the first entry in the column */
2523 ROW_CHAIN_DEL(row_chain); /* col_chain is empty */
2524 _AI_free_dof_matrix(row_chain);
2525 }
2526 /* delete the first column */
2527 COL_CHAIN_FOREACH_SAFE(col_chain, col_next, mat, DOF_MATRIX) {
2528 COL_CHAIN_DEL(col_chain);
2529 _AI_free_dof_matrix(col_chain);
2530 }
2531 /* finally, delete the (1,1)-entry */
2532 _AI_free_dof_matrix(mat);
2533 }
2534
2535
ALBERTA_DEFUNUSED(static void delete_unconnected_matrices (void))2536 ALBERTA_DEFUNUSED(static void delete_unconnected_matrices(void))
2537 {
2538 if(unconnected_matrices) deleteObject(unconnected_matrices);
2539 unconnected_matrices = NULL;
2540 }
2541
2542 /*--------------------------------------------------------------------------*/
2543 /* memory management for dof_vectors */
2544 /* unconnected dof_vectors are handled as above */
2545 /*--------------------------------------------------------------------------*/
2546
2547 #define DEFUN_ALLOC_DOF_VEC(TYPE, typename) \
2548 static void *unconnected_##typename##_vecs = NULL; \
2549 \
2550 static inline DOF_##TYPE##_VEC * \
2551 _AI_get_##typename##_vec(const char *name, const FE_SPACE *fe_space) \
2552 { \
2553 DOF_##TYPE##_VEC *vec; \
2554 static void *vec_info = NULL; \
2555 \
2556 if(!fe_space || !fe_space->admin) { \
2557 if(!unconnected_##typename##_vecs) { \
2558 unconnected_##typename##_vecs = \
2559 newObject(sizeof(DOF_##TYPE##_VEC), 10, \
2560 "unconnected "#typename" vecs"); \
2561 vec_info = unconnected_##typename##_vecs; \
2562 } \
2563 } else { \
2564 DOF_ADMIN_MEM_INFO *mem_info = \
2565 (DOF_ADMIN_MEM_INFO *)fe_space->admin->mem_info; \
2566 \
2567 vec_info = mem_info->typename##_vec; \
2568 } \
2569 \
2570 vec = (DOF_##TYPE##_VEC *)getMemory(vec_info); \
2571 \
2572 vec->next = NULL; \
2573 vec->fe_space = fe_space; \
2574 vec->name = name ? strdup(name): NULL; \
2575 vec->size = 0; \
2576 vec->reserved = \
2577 ((sizeof(TYPE##_VEC_TYPE) + sizeof(REAL) - 1) / sizeof(REAL)); \
2578 vec->vec = NULL; \
2579 vec->refine_interpol = \
2580 vec->coarse_restrict = NULL; \
2581 vec->user_data = NULL; \
2582 vec->vec_loc = NULL; \
2583 vec->mem_info = vec_info; \
2584 CHAIN_INIT(vec); \
2585 vec->unchained = NULL; /* for now */ \
2586 \
2587 if (fe_space && fe_space->admin) { \
2588 add_##typename##_vec_to_admin(vec, (DOF_ADMIN *)fe_space->admin); \
2589 } \
2590 \
2591 return vec; \
2592 } \
2593 struct _AI_semicolon_dummy
2594
2595 #define DEFUN_FREE_DOF_VEC(TYPE, type, typename) \
2596 static inline void _AI_free_##typename##_vec(DOF_##TYPE##_VEC *vec) \
2597 { \
2598 if (vec->fe_space && vec->fe_space->admin) { \
2599 remove_##typename##_vec_from_admin(vec); \
2600 } \
2601 \
2602 MEM_FREE(vec->vec, vec->size, type); \
2603 if (vec->name) { \
2604 free((char *)vec->name); \
2605 } \
2606 \
2607 if (vec->mem_info) { \
2608 freeMemory((void *)vec, vec->mem_info); \
2609 } else { \
2610 memset(vec, 0, sizeof(*vec)); \
2611 } \
2612 } \
2613 \
2614 ALBERTA_DEFUNUSED(static void delete_unconnected_##typename##_vecs(void)) \
2615 { \
2616 if(unconnected_##typename##_vecs) { \
2617 deleteObject(unconnected_##typename##_vecs); \
2618 } \
2619 unconnected_##typename##_vecs = NULL; \
2620 } \
2621 struct _AI_semicolon_dummy
2622
2623
2624 #define DEFUN_ALLOC_DOF_VEC_STD(TYPE, type, typename) \
2625 DOF_##TYPE##_VEC * \
2626 get_##typename##_vec(const char *name, const FE_SPACE *fe_space) \
2627 { \
2628 DOF_##TYPE##_VEC *vec, *vec_chain; \
2629 const FE_SPACE *fe_chain; \
2630 EL_##TYPE##_VEC *vec_loc = NULL; \
2631 \
2632 vec = _AI_get_##typename##_vec(name, fe_space); \
2633 if (fe_space != NULL) { \
2634 fe_space = copy_fe_space(fe_space); \
2635 vec->fe_space = fe_space; \
2636 if (fe_space->bas_fcts) { \
2637 vec->vec_loc = vec_loc = get_el_##type##_vec(fe_space->bas_fcts); \
2638 } \
2639 CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) { \
2640 vec_chain = _AI_get_##typename##_vec(name, fe_chain); \
2641 CHAIN_ADD_TAIL(vec, vec_chain); \
2642 if (vec_loc) { \
2643 vec_loc = CHAIN_NEXT(vec_loc, EL_##TYPE##_VEC); \
2644 vec_chain->vec_loc = vec_loc; \
2645 } \
2646 } \
2647 } \
2648 \
2649 return vec; \
2650 } \
2651 struct _AI_semicolon_dummy
2652
2653
2654 #define DEFUN_FREE_DOF_VEC_STD(TYPE, type, typename) \
2655 void free_##typename##_vec(DOF_##TYPE##_VEC *vec) \
2656 { \
2657 DOF_##TYPE##_VEC *vec_chain; \
2658 DBL_LIST_NODE *next; \
2659 const FE_SPACE *fe_space = vec->fe_space; \
2660 \
2661 if (vec->vec_loc) { \
2662 free_el_##type##_vec(vec->vec_loc); \
2663 } \
2664 CHAIN_FOREACH_SAFE(vec_chain, next, vec, DOF_##TYPE##_VEC) { \
2665 _AI_free_##typename##_vec(vec_chain); \
2666 } \
2667 _AI_free_##typename##_vec(vec); \
2668 if (fe_space) { \
2669 free_fe_space(fe_space); \
2670 } \
2671 } \
2672 struct _AI_semicolon_dummy
2673
2674 DEFUN_ALLOC_DOF_VEC(INT, dof_int);
2675 DEFUN_FREE_DOF_VEC(INT, int, dof_int);
2676 DEFUN_ALLOC_DOF_VEC_STD(INT, int, dof_int);
2677 DEFUN_FREE_DOF_VEC_STD(INT, int, dof_int);
2678
2679 DEFUN_ALLOC_DOF_VEC(DOF, dof_dof);
2680 DEFUN_FREE_DOF_VEC(DOF, DOF, dof_dof);
2681 DEFUN_ALLOC_DOF_VEC_STD(DOF, dof, dof_dof);
2682 DEFUN_FREE_DOF_VEC_STD(DOF, dof, dof_dof);
2683
2684 DEFUN_ALLOC_DOF_VEC(DOF, int_dof);
2685 DEFUN_FREE_DOF_VEC(DOF, DOF, int_dof);
2686 DEFUN_ALLOC_DOF_VEC_STD(DOF, dof, int_dof);
2687 DEFUN_FREE_DOF_VEC_STD(DOF, dof, int_dof);
2688
2689 DEFUN_ALLOC_DOF_VEC(UCHAR, dof_uchar);
2690 DEFUN_FREE_DOF_VEC(UCHAR, U_CHAR, dof_uchar);
2691 DEFUN_ALLOC_DOF_VEC_STD(UCHAR, uchar, dof_uchar);
2692 DEFUN_FREE_DOF_VEC_STD(UCHAR, uchar, dof_uchar);
2693
2694 DEFUN_ALLOC_DOF_VEC(SCHAR, dof_schar);
2695 DEFUN_FREE_DOF_VEC(SCHAR, S_CHAR, dof_schar);
2696 DEFUN_ALLOC_DOF_VEC_STD(SCHAR, schar, dof_schar);
2697 DEFUN_FREE_DOF_VEC_STD(SCHAR, schar, dof_schar);
2698
2699 DEFUN_ALLOC_DOF_VEC(PTR, dof_ptr);
2700 DEFUN_FREE_DOF_VEC(PTR, void *, dof_ptr);
2701 DEFUN_ALLOC_DOF_VEC_STD(PTR, ptr, dof_ptr);
2702 DEFUN_FREE_DOF_VEC_STD(PTR, ptr, dof_ptr);
2703
2704 DEFUN_ALLOC_DOF_VEC(REAL, dof_real);
2705 DEFUN_FREE_DOF_VEC(REAL, REAL, dof_real);
2706 DEFUN_ALLOC_DOF_VEC_STD(REAL, real, dof_real);
2707 DEFUN_FREE_DOF_VEC_STD(REAL, real, dof_real);
2708
2709 DEFUN_ALLOC_DOF_VEC(REAL_D, dof_real_d);
2710 DEFUN_FREE_DOF_VEC(REAL_D, REAL_D, dof_real_d);
2711 DEFUN_ALLOC_DOF_VEC_STD(REAL_D, real_d, dof_real_d);
2712 DEFUN_FREE_DOF_VEC_STD(REAL_D, real_d, dof_real_d);
2713
2714 DEFUN_ALLOC_DOF_VEC(REAL_DD, dof_real_dd);
2715 DEFUN_FREE_DOF_VEC(REAL_DD, REAL_DD, dof_real_dd);
2716 DEFUN_ALLOC_DOF_VEC_STD(REAL_DD, real_dd, dof_real_dd);
2717 DEFUN_FREE_DOF_VEC_STD(REAL_DD, real_dd, dof_real_dd);
2718
2719 /* The DOF_REAL_VEC_D is special, it serves to encapsulate discrete
2720 * functions -- possibly living on a chain of FE-spaces -- which may
2721 * or may not be vector-valued. We look at BAS_FCTS::rdim and
2722 * FE_SPACE::rdim to decide this:
2723 *
2724 * - BAS_FCTS::rdim == DIM_OF_WORLD (== FE_SPACE::rdim)
2725 * allocate an ordinary DOF_REAL_VEC
2726 * - BAS_FCTS::rdim == 1 == FE_SPACE::rdim
2727 * allocate an ordinary DOF_REAL_VEC
2728 * - BAS_FCTS::rdim == 1, FE_SPACE::rdim == DIM_OF_WORLD
2729 * Allocate a DOF_REAL_D_VEC. This happens when chaining a Cartesian
2730 * product space with "really" vector-valued basis functions.
2731 */
2732 DOF_REAL_VEC_D *
get_dof_real_vec_d(const char * name,const FE_SPACE * fe_space)2733 get_dof_real_vec_d(const char *name, const FE_SPACE *fe_space)
2734 {
2735 FUNCNAME("get_dof_real_vec_d");
2736 DOF_REAL_VEC_D *vec, *vec_chain;
2737 const FE_SPACE *fe_chain;
2738 EL_REAL_VEC_D *vec_loc;
2739
2740 fe_space = copy_fe_space(fe_space);
2741 if (fe_space->rdim == DIM_OF_WORLD &&
2742 fe_space->bas_fcts->rdim == DIM_OF_WORLD) {
2743 vec = (DOF_REAL_VEC_D *)_AI_get_dof_real_vec(name, fe_space);
2744 } else if (fe_space->bas_fcts->rdim == 1 && fe_space->rdim == DIM_OF_WORLD) {
2745 vec = (DOF_REAL_VEC_D *)_AI_get_dof_real_d_vec(name, fe_space);
2746 } else {
2747 ERROR_EXIT("The combination FE_SPACE::rdim == %d and "
2748 "FE_SPACE::BAS_FCTS::rdim == %d does not make sense\n",
2749 fe_space->rdim , fe_space->bas_fcts->rdim);
2750 vec_chain = vec = NULL; /* not reached */
2751 }
2752 vec->vec_loc = vec_loc = get_el_real_vec_d(fe_space->bas_fcts);
2753 CHAIN_FOREACH(fe_chain, fe_space, const FE_SPACE) {
2754 if (fe_chain->rdim == fe_chain->bas_fcts->rdim) {
2755 vec_chain = (DOF_REAL_VEC_D *)_AI_get_dof_real_vec(name, fe_chain);
2756 } else if (fe_chain->bas_fcts->rdim == 1 &&
2757 fe_chain->rdim == DIM_OF_WORLD) {
2758 vec_chain = (DOF_REAL_VEC_D *)_AI_get_dof_real_d_vec(name, fe_chain);
2759 } else {
2760 ERROR_EXIT("The combination FE_SPACE::rdim == %d and "
2761 "FE_SPACE::BAS_FCTS::rdim == %d does not make sense\n",
2762 fe_chain->rdim , fe_chain->bas_fcts->rdim);
2763 }
2764 CHAIN_ADD_TAIL(vec, vec_chain);
2765 if (vec_loc) {
2766 vec_loc = CHAIN_NEXT(vec_loc, EL_REAL_VEC_D);
2767 vec_chain->vec_loc = vec_loc;
2768 }
2769 }
2770 return vec;
2771 }
2772
free_dof_real_vec_d(DOF_REAL_VEC_D * vec)2773 void free_dof_real_vec_d(DOF_REAL_VEC_D *vec)
2774 {
2775 FUNCNAME("free_dof_real_vec_d");
2776 DOF_REAL_VEC_D *vec_chain;
2777 DBL_LIST_NODE *next;
2778 const FE_SPACE *fe_space;
2779 const BAS_FCTS *bas_fcts;
2780
2781 if (vec->vec_loc) {
2782 free_el_real_vec_d(vec->vec_loc);
2783 }
2784 CHAIN_FOREACH_SAFE(vec_chain, next, vec, DOF_REAL_VEC_D) {
2785 fe_space = vec_chain->fe_space;
2786 bas_fcts = fe_space->bas_fcts;
2787 if (vec_chain->stride == 1) {
2788 _AI_free_dof_real_vec((DOF_REAL_VEC *)vec_chain);
2789 } else if (vec_chain->stride == DIM_OF_WORLD) {
2790 _AI_free_dof_real_d_vec((DOF_REAL_D_VEC *)vec_chain);
2791 } else {
2792 ERROR_EXIT("The combination FE_SPACE::rdim == %d and "
2793 "FE_SPACE::BAS_FCTS::rdim == %d and "
2794 "EL_REAL_VEC::stride == %d does not make sense\n",
2795 fe_space->rdim , bas_fcts->rdim, vec_chain->stride);
2796 }
2797 }
2798 fe_space = vec->fe_space;
2799 bas_fcts = fe_space->bas_fcts;
2800 if (vec->stride == 1) {
2801 _AI_free_dof_real_vec((DOF_REAL_VEC *)vec);
2802 } else if (vec->stride == DIM_OF_WORLD) {
2803 _AI_free_dof_real_d_vec((DOF_REAL_D_VEC *)vec);
2804 } else {
2805 ERROR_EXIT("The combination FE_SPACE::rdim == %d and "
2806 "FE_SPACE::BAS_FCTS::rdim == %d and "
2807 "EL_REAL_VEC::stride == %d does not make sense\n",
2808 fe_space->rdim , bas_fcts->rdim, vec->stride);
2809 }
2810 free_fe_space(fe_space);
2811 }
2812
2813 #define DEFUN_ALLOC_EL_VEC(VECNAME, vecname, vectype) \
2814 static inline EL_##VECNAME##_VEC * \
2815 _AI_get_el_##vecname##_vec(const BAS_FCTS *bas_fcts) \
2816 { \
2817 EL_##VECNAME##_VEC *el_vec; \
2818 \
2819 el_vec = (EL_##VECNAME##_VEC *) \
2820 MEM_CALLOC(sizeof(EL_##VECNAME##_VEC) \
2821 + (bas_fcts->n_bas_fcts_max - 1) * sizeof(vectype), \
2822 char); \
2823 el_vec->n_components = bas_fcts->n_bas_fcts; \
2824 el_vec->n_components_max = bas_fcts->n_bas_fcts_max; \
2825 el_vec->reserved = \
2826 ((sizeof(EL_##VECNAME##_VEC_TYPE) + sizeof(REAL) - 1) / sizeof(REAL)); \
2827 CHAIN_INIT(el_vec); \
2828 \
2829 return el_vec; \
2830 } \
2831 struct _AI_semicolon_dummy
2832
2833 #define DEFUN_ALLOC_EL_VEC_STD(VECNAME, vecname, vectype) \
2834 EL_##VECNAME##_VEC *get_el_##vecname##_vec(const BAS_FCTS *bas_fcts) \
2835 { \
2836 const BAS_FCTS *bf_chain; \
2837 EL_##VECNAME##_VEC *el_vec, *el_vec_chain; \
2838 \
2839 el_vec = _AI_get_el_##vecname##_vec(bas_fcts); \
2840 CHAIN_FOREACH(bf_chain, bas_fcts, const BAS_FCTS) { \
2841 el_vec_chain = _AI_get_el_##vecname##_vec(bf_chain); \
2842 CHAIN_ADD_TAIL(el_vec, el_vec_chain); \
2843 } \
2844 \
2845 return el_vec; \
2846 } \
2847 struct _AI_semicolon_dummy
2848
2849 #define DEFUN_FREE_EL_VEC(VECNAME, vecname, vectype) \
2850 void free_el_##vecname##_vec(EL_##VECNAME##_VEC *el_vec) \
2851 { \
2852 EL_##VECNAME##_VEC *el_vec_chain; \
2853 DBL_LIST_NODE *next; \
2854 \
2855 if (el_vec == NULL) { \
2856 return; \
2857 } \
2858 CHAIN_FOREACH_SAFE(el_vec_chain, next, el_vec, EL_##VECNAME##_VEC) { \
2859 CHAIN_DEL(el_vec_chain); \
2860 MEM_FREE(el_vec_chain, \
2861 (el_vec_chain->n_components_max - 1) * sizeof(vectype) \
2862 + \
2863 sizeof(EL_##VECNAME##_VEC), char); \
2864 } \
2865 MEM_FREE(el_vec, \
2866 (el_vec->n_components_max - 1) * sizeof(vectype) \
2867 + \
2868 sizeof(EL_##VECNAME##_VEC), char); \
2869 } \
2870 struct _AI_semicolon_dummy
2871
2872 DEFUN_ALLOC_EL_VEC(INT, int, int);
2873 DEFUN_FREE_EL_VEC(INT, int, int);
2874 DEFUN_ALLOC_EL_VEC_STD(INT, int, int);
2875
2876 DEFUN_ALLOC_EL_VEC(DOF, dof, DOF);
2877 DEFUN_FREE_EL_VEC(DOF, dof, DOF);
2878 DEFUN_ALLOC_EL_VEC_STD(DOF, dof, DOF);
2879
2880 DEFUN_ALLOC_EL_VEC(UCHAR, uchar, U_CHAR);
2881 DEFUN_FREE_EL_VEC(UCHAR, uchar, U_CHAR);
2882 DEFUN_ALLOC_EL_VEC_STD(UCHAR, uchar, U_CHAR);
2883
2884 DEFUN_ALLOC_EL_VEC(SCHAR, schar, S_CHAR);
2885 DEFUN_FREE_EL_VEC(SCHAR, schar, S_CHAR);
2886 DEFUN_ALLOC_EL_VEC_STD(SCHAR, schar, S_CHAR);
2887
2888 DEFUN_ALLOC_EL_VEC(PTR, ptr, void *);
2889 DEFUN_FREE_EL_VEC(PTR, ptr, void *);
2890 DEFUN_ALLOC_EL_VEC_STD(PTR, ptr, void *);
2891
2892 DEFUN_ALLOC_EL_VEC(REAL, real, REAL);
2893 DEFUN_FREE_EL_VEC(REAL, real, REAL);
2894 DEFUN_ALLOC_EL_VEC_STD(REAL, real, REAL);
2895
2896 DEFUN_ALLOC_EL_VEC(REAL_D, real_d, REAL_D);
2897 DEFUN_FREE_EL_VEC(REAL_D, real_d, REAL_D);
2898 DEFUN_ALLOC_EL_VEC_STD(REAL_D, real_d, REAL_D);
2899
2900 DEFUN_ALLOC_EL_VEC(REAL_DD, real_dd, REAL_DD);
2901 DEFUN_FREE_EL_VEC(REAL_DD, real_dd, REAL_DD);
2902 DEFUN_ALLOC_EL_VEC_STD(REAL_DD, real_dd, REAL_DD);
2903
2904 DEFUN_ALLOC_EL_VEC(BNDRY, bndry, BNDRY_FLAGS);
2905 DEFUN_FREE_EL_VEC(BNDRY, bndry, BNDRY_FLAGS);
2906 DEFUN_ALLOC_EL_VEC_STD(BNDRY, bndry, BNDRY_FLAGS);
2907
get_el_real_vec_d(const BAS_FCTS * bas_fcts)2908 EL_REAL_VEC_D *get_el_real_vec_d(const BAS_FCTS *bas_fcts)
2909 {
2910 FUNCNAME("get_el_real_vec_d");
2911 EL_REAL_VEC_D *el_vec, *el_vec_chain;
2912 const BAS_FCTS *bfcts_chain;
2913
2914 if (bas_fcts->rdim == DIM_OF_WORLD) {
2915 el_vec = (EL_REAL_VEC_D *)_AI_get_el_real_vec(bas_fcts);
2916 } else if (bas_fcts->rdim == 1) {
2917 el_vec = (EL_REAL_VEC_D *)_AI_get_el_real_d_vec(bas_fcts);
2918 el_vec->stride = DIM_OF_WORLD;
2919 } else {
2920 ERROR_EXIT("BAS_FCTS::rdim %d not in { 1, DIM_OF_WORLD = %d }.\n",
2921 bas_fcts->rdim, DIM_OF_WORLD);
2922 el_vec = NULL; /* not reached */
2923 }
2924 CHAIN_FOREACH(bfcts_chain, bas_fcts, const BAS_FCTS) {
2925 if (bfcts_chain->rdim == DIM_OF_WORLD) {
2926 el_vec_chain = (EL_REAL_VEC_D *)_AI_get_el_real_vec(bfcts_chain);
2927 } else if (bfcts_chain->rdim == 1) {
2928 el_vec_chain = (EL_REAL_VEC_D *)_AI_get_el_real_d_vec(bfcts_chain);
2929 el_vec_chain->stride = DIM_OF_WORLD;
2930 } else {
2931 ERROR_EXIT("BAS_FCTS::rdim %d not in { 1, DIM_OF_WORLD = %d }.\n",
2932 bfcts_chain->rdim, DIM_OF_WORLD);
2933 }
2934 CHAIN_ADD_TAIL(el_vec, el_vec_chain);
2935 }
2936
2937 return el_vec;
2938 }
2939
free_el_real_vec_d(EL_REAL_VEC_D * el_vec)2940 void free_el_real_vec_d(EL_REAL_VEC_D *el_vec)
2941 {
2942 EL_REAL_VEC_D *el_vec_chain;
2943 DBL_LIST_NODE *next;
2944
2945 if (el_vec == NULL) {
2946 return;
2947 }
2948
2949 CHAIN_FOREACH_SAFE(el_vec_chain, next, el_vec, EL_REAL_VEC_D) {
2950 CHAIN_DEL(el_vec_chain);
2951 MEM_FREE(el_vec_chain,
2952 (el_vec_chain->n_components_max - 1)
2953 *
2954 el_vec_chain->stride * sizeof(REAL)
2955 +
2956 sizeof(EL_REAL_VEC), char);
2957 }
2958 MEM_FREE(el_vec,
2959 (el_vec->n_components_max - 1)
2960 *
2961 el_vec->stride * sizeof(REAL)
2962 +
2963 sizeof(EL_REAL_VEC), char);
2964 }
2965
2966 /******************************************************************************/
2967
2968 /* allocation of element matrices */
2969 static inline EL_MATRIX *
_AI_get_el_matrix_single(const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space,MATENT_TYPE operator_type)2970 _AI_get_el_matrix_single(const FE_SPACE *row_fe_space,
2971 const FE_SPACE *col_fe_space,
2972 MATENT_TYPE operator_type)
2973 {
2974 EL_MATRIX *el_mat;
2975 const BAS_FCTS *row_bfcts, *col_bfcts;
2976
2977 row_bfcts = row_fe_space->bas_fcts;
2978 col_bfcts = col_fe_space->bas_fcts;
2979
2980 el_mat = MEM_ALLOC(1, EL_MATRIX);
2981 el_mat->type = matent_type(row_fe_space, col_fe_space, operator_type);
2982
2983 el_mat->n_row = row_bfcts->n_bas_fcts;
2984 el_mat->n_col = col_bfcts->n_bas_fcts;
2985 el_mat->n_row_max = row_bfcts->n_bas_fcts_max;
2986 el_mat->n_col_max = col_bfcts->n_bas_fcts_max;
2987
2988 ROW_CHAIN_INIT(el_mat);
2989 COL_CHAIN_INIT(el_mat);
2990
2991 #undef MAT_BODY
2992 #define MAT_BODY(F, CC, C, S, TYPE) \
2993 el_mat->data.S = (TYPE *const*) \
2994 MAT_ALLOC(el_mat->n_row_max, el_mat->n_col_max, TYPE)
2995 MAT_EMIT_BODY_SWITCH(el_mat->type);
2996
2997 return el_mat;
2998 }
2999
get_el_matrix(const FE_SPACE * row_fe_space,const FE_SPACE * col_fe_space,MATENT_TYPE op_type)3000 EL_MATRIX *get_el_matrix(const FE_SPACE *row_fe_space,
3001 const FE_SPACE *col_fe_space,
3002 MATENT_TYPE op_type)
3003 {
3004 EL_MATRIX *el_mat, *row_chain, *col_chain;
3005 const FE_SPACE *row_fesp, *col_fesp;
3006
3007 if (col_fe_space == NULL) {
3008 col_fe_space = row_fe_space;
3009 }
3010 /* generate the first row */
3011 el_mat = _AI_get_el_matrix_single(row_fe_space, col_fe_space, op_type);
3012 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
3013 row_chain = _AI_get_el_matrix_single(row_fe_space, col_fesp, op_type);
3014 ROW_CHAIN_ADD_TAIL(el_mat, row_chain);
3015 }
3016 /* Then generate all following rows */
3017 CHAIN_FOREACH(row_fesp, row_fe_space, const FE_SPACE) {
3018 col_chain = _AI_get_el_matrix_single(row_fesp, col_fe_space, op_type);
3019 COL_CHAIN_ADD_TAIL(el_mat, col_chain);
3020 CHAIN_FOREACH(col_fesp, col_fe_space, const FE_SPACE) {
3021 row_chain = _AI_get_el_matrix_single(row_fesp, col_fesp, op_type);
3022 ROW_CHAIN_ADD_TAIL(col_chain, row_chain);
3023 el_mat = ROW_CHAIN_NEXT(el_mat, EL_MATRIX);
3024 COL_CHAIN_ADD_TAIL(el_mat, row_chain);
3025 }
3026 el_mat = ROW_CHAIN_NEXT(el_mat, EL_MATRIX); /* roll over to first element */
3027 }
3028 return el_mat;
3029 }
3030
_AI_free_el_matrix_single(EL_MATRIX * el_mat)3031 static inline void _AI_free_el_matrix_single(EL_MATRIX *el_mat)
3032 {
3033 #undef MAT_BODY
3034 #define MAT_BODY(F, CC, C, S, TYPE) \
3035 MAT_FREE(el_mat->data.S, \
3036 el_mat->n_row_max, el_mat->n_col_max, TYPE)
3037 MAT_EMIT_BODY_SWITCH(el_mat->type);
3038 MEM_FREE(el_mat, 1, EL_MATRIX);
3039 }
3040
free_el_matrix(EL_MATRIX * el_mat)3041 void free_el_matrix(EL_MATRIX *el_mat)
3042 {
3043 DBL_LIST_NODE *row_next, *col_next;
3044 EL_MATRIX *row_chain, *col_chain;
3045
3046 /* delete all columns, safe the first one */
3047 ROW_CHAIN_FOREACH_SAFE(row_chain, row_next, el_mat, EL_MATRIX) {
3048 /* delete the column down to the pre-last entry */
3049 COL_CHAIN_FOREACH_SAFE(col_chain, col_next, row_chain, EL_MATRIX) {
3050 ROW_CHAIN_DEL(col_chain);
3051 COL_CHAIN_DEL(col_chain);
3052 _AI_free_el_matrix_single(col_chain);
3053 }
3054 /* delete the first entry in the column */
3055 ROW_CHAIN_DEL(row_chain); /* col_chain is empty */
3056 _AI_free_el_matrix_single(row_chain);
3057 }
3058 /* delete the first column */
3059 COL_CHAIN_FOREACH_SAFE(col_chain, col_next, el_mat, EL_MATRIX) {
3060 COL_CHAIN_DEL(col_chain);
3061 _AI_free_el_matrix_single(col_chain);
3062 }
3063 /* finally, delete the (1,1)-entry */
3064 _AI_free_el_matrix_single(el_mat);
3065 }
3066
3067 static inline
__print_el_matrix(const EL_MATRIX * el_mat)3068 void __print_el_matrix(const EL_MATRIX *el_mat)
3069 {
3070 FUNCNAME("print_el_matrix");
3071 int i, j, k;
3072
3073 switch (el_mat->type) {
3074 case MATENT_REAL:
3075 for (i = 0; i < el_mat->n_row; i++) {
3076 MSG("%2d: ", i);
3077 for (j = 0; j < el_mat->n_col; j++) {
3078 print_msg(" %.8e", el_mat->data.real[i][j]);
3079 }
3080 print_msg("\n");
3081 }
3082 break;
3083 case MATENT_REAL_D:
3084 for (i = 0; i < el_mat->n_row; i++) {
3085 MSG("%2d: ", i);
3086 for (j = 0; j < el_mat->n_col; j++) {
3087 print_msg(" "FORMAT_DOW, EXPAND_DOW(el_mat->data.real_d[i][j]));
3088 }
3089 print_msg("\n");
3090 }
3091 break;
3092 case MATENT_REAL_DD:
3093 for (i = 0; i < el_mat->n_row; i++) {
3094 for (k = 0; k < DIM_OF_WORLD; k++) {
3095 if (k == 0) {
3096 MSG("%2d: ", i);
3097 } else {
3098 MSG(" ");
3099 }
3100 for (j = 0; j < el_mat->n_col; j++) {
3101 print_msg(" "FORMAT_DOW, EXPAND_DOW(el_mat->data.real_dd[i][j][k]));
3102 /* break; this "break" caused errors */
3103 }
3104 print_msg("\n");
3105 }
3106 print_msg("\n");
3107 }
3108 break;
3109 case MATENT_NONE:
3110 default:
3111 ERROR_EXIT("Unknown or invalid block-matrix type: %d\n", el_mat->type);
3112 break;
3113 }
3114 }
3115
print_el_matrix(const EL_MATRIX * el_mat)3116 void print_el_matrix(const EL_MATRIX *el_mat)
3117 {
3118 FUNCNAME("print_el_matrix");
3119 int i, j;
3120
3121 i = 0;
3122 COL_CHAIN_DO(el_mat, const EL_MATRIX) {
3123 j = 0;
3124 ROW_CHAIN_DO(el_mat, const EL_MATRIX) {
3125 if (!COL_CHAIN_SINGLE(el_mat) || !ROW_CHAIN_SINGLE(el_mat)) {
3126 MSG("BLOCK(%d,%d):\n", i, j);
3127 }
3128 __print_el_matrix(el_mat);
3129 ++j;
3130 } ROW_CHAIN_WHILE(el_mat, const EL_MATRIX);
3131 ++i;
3132 } COL_CHAIN_WHILE(el_mat, const EL_MATRIX);
3133 }
3134
3135
3136 static inline
__print_el_real_vec(const EL_REAL_VEC * el_vec)3137 void __print_el_real_vec(const EL_REAL_VEC *el_vec)
3138 {
3139 int i;
3140 for(i = 0; i < el_vec->n_components; i++)
3141 {
3142 print_msg(" %.8e", el_vec->vec[i]);
3143 }
3144 print_msg("\n");
3145 }
3146
3147 /* copy of print_dof_real_vec */
print_el_real_vec(const EL_REAL_VEC * vec)3148 void print_el_real_vec(const EL_REAL_VEC *vec)
3149 {
3150 FUNCNAME("print_el_real_vec");
3151 int i;
3152 i = 0;
3153
3154 CHAIN_DO(vec, const EL_REAL_VEC) {
3155 if (!CHAIN_SINGLE(vec)) {
3156 MSG("BLOCK(%d):\n", i);
3157 }
3158 __print_el_real_vec(vec);
3159 ++i;
3160 } CHAIN_WHILE(vec, const EL_REAL_VEC);
3161 }
3162
3163 static inline
__print_el_real_d_vec(const EL_REAL_D_VEC * el_vec)3164 void __print_el_real_d_vec(const EL_REAL_D_VEC *el_vec)
3165 {
3166 int i;
3167
3168 for(i = 0; i < el_vec->n_components; i++)
3169 {
3170 print_msg(" "FORMAT_DOW, EXPAND_DOW(el_vec->vec[i]));
3171 }
3172 print_msg("\n");
3173 }
3174
print_el_real_d_vec(const EL_REAL_D_VEC * vec)3175 void print_el_real_d_vec(const EL_REAL_D_VEC *vec)
3176 {
3177 FUNCNAME("print_el_real_d_vec");
3178 int i;
3179
3180 i = 0;
3181 CHAIN_DO(vec, const EL_REAL_D_VEC) {
3182 if (!CHAIN_SINGLE(vec)) {
3183 MSG("BLOCK(%d):\n", i);
3184 }
3185 __print_el_real_d_vec(vec);
3186 ++i;
3187 } CHAIN_WHILE(vec, const EL_REAL_D_VEC);
3188 }
3189
print_el_real_vec_d(const EL_REAL_VEC_D * vec)3190 void print_el_real_vec_d(const EL_REAL_VEC_D *vec)
3191 {
3192 FUNCNAME("print_el_real_d_vec");
3193 int i;
3194 i = 0;
3195
3196 CHAIN_DO(vec, const EL_REAL_VEC_D) {
3197 if (!CHAIN_SINGLE(vec)) {
3198 MSG("BLOCK(%d):\n", i);
3199 }
3200 if (vec->stride != 1){
3201 __print_el_real_d_vec((const EL_REAL_D_VEC *) vec);
3202 } else {
3203 __print_el_real_vec((const EL_REAL_VEC *) vec);
3204 }
3205 ++i;
3206 } CHAIN_WHILE(vec, const EL_REAL_VEC_D);
3207 }
3208
3209 static inline
__print_el_dof_vec(const EL_DOF_VEC * el_vec)3210 void __print_el_dof_vec(const EL_DOF_VEC *el_vec)
3211 {
3212 int i;
3213 for(i = 0; i < el_vec->n_components; i++) {
3214 print_msg(" %d", el_vec->vec[i]);
3215 }
3216 print_msg("\n");
3217 }
3218
3219 /* copy of print_dof_real_vec */
print_el_dof_vec(const EL_DOF_VEC * vec)3220 void print_el_dof_vec(const EL_DOF_VEC *vec)
3221 {
3222 FUNCNAME("print_el_dof_vec");
3223 int i;
3224 i = 0;
3225
3226 CHAIN_DO(vec, const EL_DOF_VEC) {
3227 if (!CHAIN_SINGLE(vec)) {
3228 MSG("BLOCK(%d): ", i);
3229 }
3230 __print_el_dof_vec(vec);
3231 ++i;
3232 } CHAIN_WHILE(vec, const EL_DOF_VEC);
3233 }
3234
3235 static inline
__print_el_schar_vec(const EL_SCHAR_VEC * el_vec)3236 void __print_el_schar_vec(const EL_SCHAR_VEC *el_vec)
3237 {
3238 int i;
3239 for(i = 0; i < el_vec->n_components; i++) {
3240 print_msg(" %02x", (int)el_vec->vec[i]);
3241 }
3242 print_msg("\n");
3243 }
3244
3245 /* copy of print_schar_real_vec */
print_el_schar_vec(const EL_SCHAR_VEC * vec)3246 void print_el_schar_vec(const EL_SCHAR_VEC *vec)
3247 {
3248 FUNCNAME("print_el_schar_vec");
3249 int i;
3250 i = 0;
3251
3252 CHAIN_DO(vec, const EL_SCHAR_VEC) {
3253 if (!CHAIN_SINGLE(vec)) {
3254 MSG("BLOCK(%d): ", i);
3255 }
3256 __print_el_schar_vec(vec);
3257 ++i;
3258 } CHAIN_WHILE(vec, const EL_SCHAR_VEC);
3259 }
3260
3261 static inline
__print_el_bndry_vec(const EL_BNDRY_VEC * el_vec)3262 void __print_el_bndry_vec(const EL_BNDRY_VEC *el_vec)
3263 {
3264 int i, j;
3265 for(i = 0; i < el_vec->n_components; i++) {
3266 for (j = 0; j < 4; j++) {
3267 print_msg("%lx", (long)el_vec->vec[i][j]);
3268 }
3269 print_msg(" ");
3270 }
3271 print_msg("\n");
3272 }
3273
3274 /* copy of print_bndry_real_vec */
print_el_bndry_vec(const EL_BNDRY_VEC * vec)3275 void print_el_bndry_vec(const EL_BNDRY_VEC *vec)
3276 {
3277 FUNCNAME("print_el_bndry_vec");
3278 int i;
3279 i = 0;
3280
3281 CHAIN_DO(vec, const EL_BNDRY_VEC) {
3282 if (!CHAIN_SINGLE(vec)) {
3283 MSG("BLOCK(%d): ", i);
3284 }
3285 __print_el_bndry_vec(vec);
3286 ++i;
3287 } CHAIN_WHILE(vec, const EL_BNDRY_VEC);
3288 }
3289