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: mesh-interface-common.c */
7 /* */
8 /*--------------------------------------------------------------------------*/
9 /* */
10 /* authors: Alfred Schmidt */
11 /* Zentrum fuer Technomathematik */
12 /* Fachbereich 3 Mathematik/Informatik */
13 /* Universitaet Bremen */
14 /* Bibliothekstr. 2 */
15 /* D-28359 Bremen, Germany */
16 /* */
17 /* Kunibert G. Siebert */
18 /* Institut fuer Mathematik */
19 /* Universitaet Augsburg */
20 /* Universitaetsstr. 14 */
21 /* D-86159 Augsburg, Germany */
22 /* */
23 /* Claus-Justus Heine */
24 /* Abteilung fuer Angewandte Mathematik */
25 /* Albert-Ludwigs-Universitaet Freiburg */
26 /* Hermann-Herder-Str. 10 */
27 /* D-79104 Freiburg im Breisgau, Germany */
28 /* */
29 /* Robert Kloefkorn */
30 /* Abteilung fuer Angewandte Mathematik */
31 /* Albert-Ludwigs-Universitaet Freiburg */
32 /* Hermann-Herder-Str. 10 */
33 /* D-79104 Freiburg im Breisgau, Germany */
34 /* */
35 /* http://www.alberta-fem.de/ */
36 /* */
37 /*--------------------------------------------------------------------------*/
38 /* (c) by A. Schmidt and K.G. Siebert (1996-2004) */
39 /* C.-J. Heine and R. Kloefkorn (1998-2007) */
40 /* */
41 /* This program is free software; you can redistribute it and/or modify */
42 /* it under the terms of the GNU General Public License as published by */
43 /* the Free Software Foundation; either version 2 of the License, or */
44 /* any later version. */
45 /* */
46 /* This program is distributed in the hope that it will be useful, */
47 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
48 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
49 /* GNU General Public License for more details. */
50 /*--------------------------------------------------------------------------*/
51
52 /****************************************************************************/
53 /* Mesh hierarchy traversal: first(), next(), ... */
54 /****************************************************************************/
55
56 #if FE_DIM == 3
57 typedef HELEMENT3D GRAPE_ELEMENT;
58 typedef HMESH3D GRAPE_MESH;
59 typedef GENMESH3D GRAPE_GENMESH;
60 typedef F_HEL_INFO3D F_EL_INFO;
61 typedef F_HDATA3D F_DATA;
62 #define Grape_Mesh HMesh3d
63 #else
64 typedef HELEMENT2D GRAPE_ELEMENT;
65 typedef GENMESH2D GRAPE_GENMESH;
66 typedef F_HEL_INFO2D F_EL_INFO;
67
68 /* higher order visual */
69 typedef F_HPDATA2D F_DATA;
70 typedef HPMESH2D GRAPE_MESH;
71 #define Grape_Mesh HPMesh2d
72 /*
73 typedef F_HDATA2D F_DATA;
74 #define Grape_Mesh HMesh2d
75 */
76 #endif
77
78 static int traverse_level = -1;
79 static FLAGS traverse_fill_flag = FILL_COORDS | FILL_BOUND | FILL_NEIGH
80 #if FE_DIM==3
81 | FILL_ORIENTATION
82 #endif
83 | CALL_LEAF_EL;
84 static const EL_INFO *el_info;
85
86 /* coord has to be double [n][3], because hmesh in grape is always
87 * of dimension 3 and for FE_DIM == 2 we set the third component to zero */
88 static double coord[N_VERTICES_LIMIT][3];
89 static int vindex[N_VERTICES_LIMIT];
90 static GRAPE_ELEMENT element;
91
92 static double *vertex[N_VERTICES_LIMIT] = {coord[0], coord[1],
93 coord[2], coord[3] };
94
95 /* implementation of functions */
96
setup_traverse(int level,FLAGS fill_flag)97 void setup_traverse(int level, FLAGS fill_flag)
98 {
99 traverse_level = level;
100 traverse_fill_flag = fill_flag;
101 return;
102 }
103
104 /****************************************************************************/
105
fill_grape_element(void)106 static void fill_grape_element(void)
107 {
108 int i, j;
109
110 for (i = 0; i < N_VERTICES(FE_DIM); i++)
111 {
112 vindex[i] = el_info->el->dof[i][0];
113 for (j = 0; j < DIM_OF_WORLD; j++) {
114 coord[i][j] = el_info->coord[i][j]; /* double to float conv */
115 }
116 }
117 element.vertex = (double G_CONST*G_CONST*) vertex;
118 element.vindex = vindex;
119
120 #if FE_DIM == 3
121 if (el_info->orientation > 0)
122 element.descr = &tetra_description_even;
123 else
124 element.descr = &tetra_description_odd;
125 #else
126 element.descr = &tria_description;
127 #endif
128
129 element.eindex = INDEX(el_info->el);
130 element.user_data = (void *)el_info;
131 }
132
133 /****************************************************************************/
134 static TRAVERSE_STACK *stack = NULL;
135
136 /***************************************************************************/
137
138 /* returns first leaf element, which is treated like a macro element */
139 /* because traverse_fill_flag contains the CALL_LEAF_EL flag */
first_grape_element(GRAPE_GENMESH * grape_mesh,MESH_ELEMENT_FLAGS flags)140 static GRAPE_ELEMENT *first_grape_element(GRAPE_GENMESH *grape_mesh,
141 MESH_ELEMENT_FLAGS flags)
142 {
143 if (!stack) stack = get_traverse_stack();
144 el_info = traverse_first(stack, (MESH *)(grape_mesh->user_data),
145 traverse_level, traverse_fill_flag);
146 if (el_info)
147 {
148 fill_grape_element();
149 element.mesh = (GRAPE_GENMESH *) grape_mesh;
150 return(&element);
151 }
152 else
153 return(NULL);
154 }
155
156 /****************************************************************************/
157
158 /* returns next leaf element, which is treated like a macro element */
159 /* because traverse_fill_flag contains the CALL_LEAF_EL flag */
next_grape_element(GRAPE_ELEMENT * el,MESH_ELEMENT_FLAGS flags)160 static GRAPE_ELEMENT *next_grape_element(GRAPE_ELEMENT *el,
161 MESH_ELEMENT_FLAGS flags)
162 {
163 if ((el_info = traverse_next(stack, (EL_INFO *)el->user_data)))
164 {
165 fill_grape_element();
166 return(&element);
167 }
168 else
169 return(NULL);
170 }
171
172 /* we have no children */
fake_child(GRAPE_ELEMENT * el,MESH_ELEMENT_FLAGS flags)173 static GRAPE_ELEMENT *fake_child(GRAPE_ELEMENT *el,MESH_ELEMENT_FLAGS flags)
174 {
175 return (NULL);
176 }
177 /* we dont want to select no child */
fake_select(GRAPE_ELEMENT * el,double * parent_coord,double * child_coord,MESH_ELEMENT_FLAGS flags)178 static GRAPE_ELEMENT *fake_select(GRAPE_ELEMENT *el,double *parent_coord,
179 double * child_coord, MESH_ELEMENT_FLAGS flags)
180 {
181 return (NULL);
182 }
183
184 /****************************************************************************/
complete_grape_element(GRAPE_ELEMENT * el,MESH_ELEMENT_FLAGS flags)185 static GRAPE_ELEMENT *complete_grape_element(GRAPE_ELEMENT *el,
186 MESH_ELEMENT_FLAGS flags)
187 {
188 if (el_info)
189 return(&element);
190 else
191 return(NULL);
192 }
193
194 /* mkae copy of element is not used, because we are acting on leaf level */
copy_grape_element(GRAPE_ELEMENT * el,MESH_ELEMENT_FLAGS flags)195 static GRAPE_ELEMENT *copy_grape_element(GRAPE_ELEMENT *el,
196 MESH_ELEMENT_FLAGS flags)
197 {
198 FUNCNAME("copy_grape_element");
199
200 MSG("not implemented yet!\n");
201 return ( NULL );
202 }
203
204 /****************************************************************************/
205 /* not used, because copy_element is not implemented */
free_grape_element(GRAPE_ELEMENT * el)206 static void free_grape_element(GRAPE_ELEMENT *el)
207 {
208 return;
209 }
210
211 /****************************************************************************/
212
f_bounds(GRAPE_ELEMENT * el,double * min,double * max,void * function_data)213 static void f_bounds(GRAPE_ELEMENT *el, double* min, double* max,
214 void *function_data)
215 {
216 (*min) = 1.0E+308;
217 (*max) = -1.0E+308;
218 }
219
220 /****************************************************************************/
221
grape_get_vertex_estimate(GRAPE_ELEMENT * el,double * value,void * function_data)222 static void grape_get_vertex_estimate(GRAPE_ELEMENT *el, double *value,
223 void *function_data)
224 {
225 *value = 0.0;
226 return;
227 }
228
229 /****************************************************************************/
230
grape_get_element_estimate(GRAPE_ELEMENT * el,void * function_data)231 static double grape_get_element_estimate(GRAPE_ELEMENT *el, void *function_data)
232 {
233 return 0.0;
234 }
235
236
237 /****************************************************************************/
f_real(GRAPE_ELEMENT * el,int ind,double G_CONST * coord,double * val,void * function_data)238 static void f_real(GRAPE_ELEMENT *el, int ind, double G_CONST *coord,
239 double *val, void *function_data)
240 {
241 FUNCNAME("f_real");
242 EL_INFO *el_info = (EL_INFO *)(el->user_data);
243 DOF_REAL_VEC *u_h;
244 const EL_REAL_VEC *uh_loc;
245 const BAS_FCTS *bas_fcts;
246
247 static REAL_B vert_coords[N_VERTICES_LIMIT] = {
248 INIT_BARY_0D(1.0),
249 INIT_BARY_1D(0.0, 1.0),
250 INIT_BARY_2D(0.0, 0.0, 1.0),
251 INIT_BARY_3D(0.0,0.0,0.0,1.0)
252 };
253
254 u_h = (DOF_REAL_VEC *)(function_data);
255 if (u_h) {
256 if (u_h->fe_space && (bas_fcts = u_h->fe_space->bas_fcts)) {
257 uh_loc = fill_el_real_vec(NULL, el_info->el, u_h);
258 if (coord) {
259 val[0] = eval_uh(coord, uh_loc, bas_fcts);
260 } else {
261 if ((ind >= 0) && (ind < N_VERTICES(FE_DIM))) {
262 val[0] = eval_uh(vert_coords[ind], uh_loc, bas_fcts);
263 } else {
264 MSG("invalid vertex number %d\n",ind);
265 val[0] = 0.0;
266 }
267 }
268 } else {
269 MSG("no fe_space or bas_fcts in dof_real_vec <%s>\n", u_h->name);
270 val[0] = 0.0;
271 }
272 } else {
273 MSG("function_data == NULL\n");
274 val[0] = 0.0;
275 }
276 }
277
278
f_real_el_info(GRAPE_ELEMENT * el,F_EL_INFO * f_el_info,void * function_data)279 static void f_real_el_info(GRAPE_ELEMENT *el, F_EL_INFO *f_el_info,
280 void *function_data)
281 {
282 FUNCNAME("f_real_el_info");
283 DOF_REAL_VEC *u_h = (DOF_REAL_VEC *)(function_data);
284
285 if (u_h && u_h->fe_space && u_h->fe_space->bas_fcts) {
286 f_el_info->polynomial_degree = u_h->fe_space->bas_fcts->degree;
287 } else {
288 ERROR("no uh or fe_space or bas_fcts\n");
289 f_el_info->polynomial_degree = 1;
290 }
291
292 return;
293 }
294
295 /****************************************************************************/
296
f_real_d(GRAPE_ELEMENT * el,int ind,double G_CONST * coord,double * val,void * function_data)297 static void f_real_d(GRAPE_ELEMENT *el, int ind, double G_CONST *coord,
298 double *val, void *function_data)
299 {
300 FUNCNAME("f_real_d");
301 EL_INFO *el_info = (EL_INFO *)(el->user_data);
302 DOF_REAL_VEC_D *u_h;
303 const EL_REAL_VEC_D *uh_loc;
304 const BAS_FCTS *bas_fcts;
305 static REAL_B vert_coords[N_VERTICES_LIMIT] = {
306 INIT_BARY_0D(1.0),
307 INIT_BARY_1D(0.0, 1.0),
308 INIT_BARY_2D(0.0, 0.0, 1.0),
309 INIT_BARY_3D(0.0,0.0,0.0,1.0)
310 };
311
312 u_h = (DOF_REAL_VEC_D *)function_data;
313 if (u_h) {
314 if (u_h->fe_space && (bas_fcts = u_h->fe_space->bas_fcts)) {
315 uh_loc = fill_el_real_vec_d(NULL, el_info->el, u_h);
316 if (coord) {
317 eval_uh_dow(val, coord, uh_loc, bas_fcts);
318 } else {
319 if ((ind >= 0) && (ind < N_VERTICES(FE_DIM))) {
320 eval_uh_dow(val, vert_coords[ind], uh_loc, bas_fcts);
321 } else {
322 MSG("invalid vertex number %d\n",ind);
323 val[0] = 0.0;
324 }
325 }
326 } else {
327 MSG("no fe_space or bas_fcts in dof_real_vec <%s>\n", u_h->name);
328 val[0] = 0.0;
329 }
330 } else {
331 MSG("no dof_real_vec as function_data\n");
332 val[0] = 0.0;
333 }
334 }
335
336
f_real_d_el_info(GRAPE_ELEMENT * el,F_EL_INFO * f_el_info,void * function_data)337 static void f_real_d_el_info(GRAPE_ELEMENT *el, F_EL_INFO *f_el_info,
338 void *function_data)
339 {
340 FUNCNAME("f_real_el_info");
341 DOF_REAL_VEC_D *u_h = (DOF_REAL_VEC_D *)function_data;
342
343 if (u_h && u_h->fe_space && u_h->fe_space->bas_fcts) {
344 f_el_info->polynomial_degree = u_h->fe_space->bas_fcts->degree;
345 } else {
346 ERROR("no uh or fe_space or bas_fcts\n");
347 f_el_info->polynomial_degree = 1;
348 }
349
350 return;
351 }
352
353 /*************************************************************************/
354
355
grape_ini_f_data(GRAPE_MESH * grape_mesh,DOF_REAL_VEC * dof_vec)356 void grape_ini_f_data(GRAPE_MESH *grape_mesh, DOF_REAL_VEC *dof_vec)
357 {
358 FUNCNAME("grape_ini_f_data");
359 F_DATA *f_data = NULL;
360
361 if (!grape_mesh)
362 {
363 MSG("no grape_mesh\n");
364 return;
365 }
366
367 if (dof_vec && dof_vec->vec)
368 {
369 for (f_data = (void *)grape_mesh->f_data; f_data;
370 f_data = (void *)f_data->next)
371 {
372 if (f_data->function_data == dof_vec)
373 break;
374 }
375
376 if (!f_data)
377 for (f_data = (void *)grape_mesh->f_data; f_data;
378 f_data = (void *)f_data->last)
379 {
380 if (f_data->function_data == dof_vec)
381 break;
382 }
383
384
385 if (!f_data)
386 {
387 MSG("generate scalar f_data for `%s'\n", dof_vec->name);
388 f_data = (F_DATA *)g_mem_alloc(sizeof(F_DATA));
389 f_data->name = (char *) dof_vec->name;
390 f_data->dimension_of_value = 1;
391 if (strstr(dof_vec->fe_space->bas_fcts->name, "disc"))
392 f_data->continuous_data = 0;
393 else
394 f_data->continuous_data = 1;
395
396 f_data->f = f_real;
397 f_data->f_el_info = f_real_el_info;
398
399 f_data->next = NULL;
400 f_data->last = NULL;
401
402 /* this is done by add-funtion
403 f_data->next = grape_mesh->f_data;
404 f_data->last = grape_mesh->f_data ? grape_mesh->f_data->last : NULL;
405 */
406
407 f_data->function_data = (void *) dof_vec;
408
409 f_data->get_bounds = f_bounds;
410 f_data->get_vertex_estimate = grape_get_vertex_estimate;
411 f_data->get_element_estimate = grape_get_element_estimate;
412 f_data->threshold = 0.0;
413 #if FE_DIM == 3
414 f_data->geometry_threshold = 0.0;
415 #else
416 f_data->get_element_p_estimates = NULL;
417 f_data->get_edge_p_estimates = NULL;
418 #endif
419 f_data->hp_threshold = 0.0;
420 f_data->hp_maxlevel = 0;
421
422 /* grape_mesh->f_data = (GENMESH_FDATA *) f_data; */
423 /* use the hamehs add-function method to add data to mesh */
424 grape_mesh = (GRAPE_MESH *) GRAPE(grape_mesh,"add-function")(f_data);
425
426 } else if (grape_mesh->f_data != (GENMESH_FDATA *)f_data) {
427 MSG("select f_data for `%s'\n", dof_vec->name);
428 grape_mesh = (GRAPE_MESH *) GRAPE(grape_mesh,"add-function")(f_data);
429 /* use the hamehs add-function method to add data to mesh */
430 /* grape_mesh->f_data = (GENMESH_FDATA *)f_data; */
431 }
432 } else {
433 MSG("no dof_vec, or no vec\n");
434 }
435
436 return;
437 }
438
grape_ini_f_data_d(GRAPE_MESH * grape_mesh,DOF_REAL_VEC_D * dof_vec)439 void grape_ini_f_data_d(GRAPE_MESH *grape_mesh, DOF_REAL_VEC_D *dof_vec)
440 {
441 FUNCNAME("grape_ini_f_data_d");
442 F_DATA *f_data=NULL;
443
444 if (!grape_mesh)
445 {
446 MSG("no grape_mesh\n");
447 return;
448 }
449
450 if (dof_vec && dof_vec->vec)
451 {
452 for (f_data = (void *)grape_mesh->f_data; f_data;
453 f_data = (void *)f_data->next)
454 {
455 if (f_data->function_data == dof_vec)
456 break;
457 }
458
459 if (!f_data)
460 for (f_data = (void *)grape_mesh->f_data; f_data;
461 f_data = (void *)f_data->last)
462 {
463 if (f_data->function_data == dof_vec)
464 break;
465 }
466
467
468 if (!f_data) {
469 MSG("generate vector-valued f_data for `%s'\n", dof_vec->name);
470 f_data = (F_DATA *)g_mem_alloc(sizeof(F_DATA));
471 f_data->name = (char *) dof_vec->name;
472 f_data->dimension_of_value = DIM_OF_WORLD;
473 if (strstr(dof_vec->fe_space->bas_fcts->name, "disc"))
474 f_data->continuous_data = 0;
475 else
476 f_data->continuous_data = 1;
477
478 f_data->f = f_real_d;
479 f_data->f_el_info = f_real_d_el_info;
480
481 f_data->next = NULL;
482 f_data->last = NULL;
483
484 /* this is done by add-function
485 f_data->next = grape_mesh->f_data;
486 f_data->last = grape_mesh->f_data ? grape_mesh->f_data->last : NULL;
487 */
488
489 f_data->function_data = (void *) dof_vec;
490
491 f_data->get_bounds = f_bounds;
492 f_data->get_vertex_estimate = grape_get_vertex_estimate;
493 f_data->get_element_estimate = grape_get_element_estimate;
494 f_data->threshold = 0.0;
495 #if FE_DIM == 3
496 f_data->geometry_threshold = 0.0;
497 #endif
498
499 f_data->hp_threshold = 0.0;
500 f_data->hp_maxlevel = 0;
501
502 /*grape_mesh->f_data = (GENMESH_FDATA *)f_data;*/
503 grape_mesh = (GRAPE_MESH *) GRAPE(grape_mesh,"add-function")(f_data);
504 } else if (grape_mesh->f_data != (GENMESH_FDATA *)f_data) {
505 MSG("select f_data for `%s'\n", dof_vec->name);
506 grape_mesh = (GRAPE_MESH *) GRAPE(grape_mesh,"add-function")(f_data);
507 /*grape_mesh->f_data = (GENMESH_FDATA *)f_data;*/
508 }
509 } else {
510 MSG("no dof_vec, or no vec\n");
511 }
512
513 return;
514 }
515
516 /****************************************************************************/
517 /* handling of multiple functions (selection by next/last) */
518 /****************************************************************************/
519
next_f_data_send(void)520 static GRAPE_MESH *next_f_data_send(void)
521 {
522 GRAPE_MESH *self;
523
524 self = (GRAPE_MESH *)START_METHOD(G_INSTANCE);
525 ASSURE(self, "", END_METHOD(NULL));
526
527 if (self->f_data && self->f_data->next)
528 {
529 self->f_data->next->last = self->f_data; /*only to be sure...*/
530 self->f_data = self->f_data->next;
531 }
532 if (self->f_data)
533 printf("new f_data is: %s\n", self->f_data->name);
534
535 END_METHOD(self);
536 }
537
prev_f_data_send(void)538 static GRAPE_MESH *prev_f_data_send(void)
539 {
540 GRAPE_MESH *self;
541
542 self = (GRAPE_MESH *)START_METHOD(G_INSTANCE);
543 ASSURE(self, "", END_METHOD(NULL));
544
545 if (self->f_data && self->f_data->last)
546 {
547 self->f_data->last->next = self->f_data; /*only to be sure...*/
548 self->f_data = self->f_data->last;
549 }
550 if (self->f_data)
551 printf("new f_data is: %s\n", self->f_data->name);
552
553 END_METHOD(self);
554 }
555
556 /* copy function data */
copyFdata(F_DATA * copy,F_DATA * org)557 static inline void copyFdata(F_DATA *copy, F_DATA *org)
558 {
559 copy->name = org->name;
560 copy->last = org->last;
561 copy->next = org->next;
562
563 copy->dimension_of_value = org->dimension_of_value;
564 copy->continuous_data = org->continuous_data;
565 copy->function_data = org->function_data;
566
567 copy->f = org->f;
568 copy->f_el_info = org->f_el_info;
569
570 copy->get_bounds = org->get_bounds;
571 copy->get_vertex_estimate = org->get_vertex_estimate;
572 copy->get_element_estimate = org->get_element_estimate;
573
574 copy->threshold = org->threshold;
575 #if DIM_OF_WORLD == 3
576 copy->geometry_threshold = org->geometry_threshold;
577 #else
578 copy->get_element_p_estimates = org->get_element_p_estimates;
579 copy->get_edge_p_estimates = org->get_edge_p_estimates;
580 #endif
581 copy->hp_threshold = org->hp_threshold;
582 copy->hp_maxlevel = org->hp_maxlevel;
583 }
584
585 /* also used in albert_movi */
copyHmeshes(GRAPE_MESH * orgMesh,GRAPE_MESH * self)586 void copyHmeshes(GRAPE_MESH *orgMesh, GRAPE_MESH *self)
587 {
588 FUNCNAME("copyHmeshes");
589
590 DEBUG_TEST_EXIT(self != NULL, "self == NULL");
591 DEBUG_TEST_EXIT(orgMesh != NULL, "orgMesh == NULL");
592
593 if (!self->f_data && orgMesh->f_data) {
594 self->level_of_interest = orgMesh->level_of_interest;
595 F_DATA *next_data = (F_DATA *) orgMesh->f_data;
596 /* to keep the same order we have to go backward */
597 while (next_data) {
598 if (next_data->next) {
599 next_data = (F_DATA *) next_data->next;
600 } else {
601 break;
602 }
603 }
604 while (next_data) {
605 F_DATA * f_data = (F_DATA *) malloc(sizeof(F_DATA));
606
607 DEBUG_TEST_EXIT(f_data != NULL, "f_data == NULL");
608 copyFdata(f_data,next_data);
609
610 self = (GRAPE_MESH *)GRAPE(self,"add-function")(f_data);
611 next_data = (F_DATA *)next_data->last;
612 }
613 }
614
615 self->max_dimension_of_coord = orgMesh->max_dimension_of_coord;
616 self->max_eindex = orgMesh->max_eindex;
617 self->max_vindex = orgMesh->max_vindex;
618 self->max_dindex = orgMesh->max_dindex;
619 self->max_number_of_vertices = orgMesh->max_number_of_vertices;
620
621 self->access_mode = orgMesh->access_mode;
622 self->access_capability = orgMesh->access_capability;
623
624 /* we have to do that, GRAPE sucks */
625 /* set other function_data pointers */
626 if (orgMesh->f_data) {
627 GENMESH_FDATA * sf = self->f_data;
628 while(sf != NULL) {
629 const char * sfname = sf->name;
630 GENMESH_FDATA * nf = orgMesh->f_data;
631 int length = strlen(sfname);
632 while( (nf != NULL) ) {
633 /* compare the real function name, ha, not with me */
634 const char * nfname = nf->name;
635 if (strncmp(sfname,nfname,length) == 0) {
636 sf->function_data = nf->function_data;
637 break;
638 }
639
640 /* GRAPE sucks, sucks, sucks
641 * Robert after debugin' for this shit more than one day */
642 if (nf != nf->last ) {
643 nf = nf->next;
644 } else {
645 break;
646 }
647 }
648
649 /* go next f_data */
650 if (sf != sf->last ) {
651 sf = sf->next;
652 }
653 }
654 }
655
656 /* copy current function selections to orgMesh */
657 self = (GRAPE_MESH *) GRAPE(self, "copy-function-selector")(orgMesh);
658
659 self->user_data = orgMesh->user_data;
660
661 self->copy_element = orgMesh->copy_element;
662 self->free_element = orgMesh->free_element;
663
664 self->complete_element = orgMesh->complete_element;
665 self->set_time = orgMesh->set_time;
666 self->get_time = orgMesh->get_time;
667
668 self->first_macro = orgMesh->first_macro;
669 self->next_macro = orgMesh->next_macro;
670 self->first_child = orgMesh->first_child;
671 self->next_child = orgMesh->next_child;
672 self->select_child = orgMesh->select_child;
673
674 self->max_level = orgMesh->max_level;
675 self->level_of_interest = orgMesh->level_of_interest;
676 /* do not set level_of_interest, because is set by user during run time */
677
678 self->get_geometry_vertex_estimate = orgMesh->get_geometry_vertex_estimate;
679 self->get_geometry_element_estimate = orgMesh->get_geometry_element_estimate;
680 self->get_lens_element_estimate = orgMesh->get_lens_element_estimate;
681 self->threshold = orgMesh->threshold;
682
683 #if DIM_OF_WORLD == 2
684 self->dimension_of_world = orgMesh->dimension_of_world;
685 #endif
686 }
687
688 /* make hard copy of hmesh, i.e. create new object and copy data */
newHmeshHardCopy(void)689 static GRAPE_MESH* newHmeshHardCopy(void)
690 {
691 GRAPE_MESH * hmesh = (GRAPE_MESH*) START_METHOD (G_INSTANCE);
692 ALERT (hmesh, "hmesh-hardcopy: No hmesh!", END_METHOD(NULL));
693 char *newName = malloc(strlen(hmesh->name)+sizeof("H: "));
694 // get new hmesh
695 GRAPE_MESH* copy;
696
697 sprintf(newName, "H: %s", hmesh->name);
698
699 copy = (GRAPE_MESH *) GRAPE(Grape_Mesh,"new-instance")(newName);
700 ALERT (copy, "hmesh_hardcopy: No new instance!", END_METHOD(NULL));
701
702 // copy mesh
703 copyHmeshes(hmesh, copy);
704
705 END_METHOD (copy);
706 }
707
grape_add_methods(void)708 static void grape_add_methods(void)
709 {
710 static bool init_done;
711
712 if (!init_done) {
713 init_done = true;
714 GRAPE(Grape_Mesh, "add-method")("next-f-data-send",next_f_data_send);
715 GRAPE(Grape_Mesh, "add-method")("prev-f-data-send",prev_f_data_send);
716
717 if (GRAPE(Grape_Mesh, "find-method")("=hardcopy")) {
718 GRAPE(Grape_Mesh, "delete-method")("hardcopy");
719 }
720 GRAPE(Grape_Mesh, "add-method")("hardcopy", &newHmeshHardCopy);
721 }
722 }
723
724 /****************************************************************************/
725 /* convert ALBERTA mesh into GRAPE mesh */
726 /****************************************************************************/
727
setup_grape_mesh(MESH * mesh,char * name)728 GRAPE_MESH *setup_grape_mesh(MESH *mesh, char *name)
729 {
730 FUNCNAME("setup_grape_mesh");
731 GRAPE_MESH *grape_mesh;
732 int i, j;
733
734 grape_add_methods();
735
736 ASSURE((grape_mesh = (GRAPE_MESH *)GRAPE(Grape_Mesh,"new-instance")(name)),
737 "can't get Grape Mesh instance",
738 return(NULL));
739
740 grape_mesh->max_level = 0;
741
742 for (i=0; i < N_VERTICES(FE_DIM); i++)
743 for (j=0; j<3; j++)
744 coord[i][j] = 0.0;
745
746 /* Funktioniert halt nur fuer alle leaf Element, gell! */
747 grape_mesh->first_child = fake_child;
748 grape_mesh->next_child = fake_child;
749 grape_mesh->first_macro = first_grape_element;
750 grape_mesh->next_macro = next_grape_element;
751 grape_mesh->copy_element = copy_grape_element;
752 grape_mesh->free_element = free_grape_element;
753
754 grape_mesh->select_child = fake_select;
755 grape_mesh->complete_element = complete_grape_element;
756
757 /* robert: hier lieber nitch NULL setzen, wer weiss?
758 grape_mesh->get_lens_element_estimate = NULL;
759 grape_mesh->get_geometry_element_estimate = NULL;
760 grape_mesh->get_geometry_vertex_estimate = NULL;
761 grape_mesh->set_time = NULL;
762 grape_mesh->get_time = NULL;
763 */
764
765
766 grape_mesh->max_number_of_vertices = N_VERTICES(FE_DIM);
767 grape_mesh->max_eindex = mesh->n_hier_elements;
768 grape_mesh->max_vindex = mesh->n_vertices;
769 grape_mesh->max_dimension_of_coord = N_VERTICES(FE_DIM);
770
771 grape_mesh->max_dindex = 1;
772
773 grape_mesh->level_of_interest = grape_mesh->max_level;
774
775 #if FE_DIM==2
776 #if (DIM_OF_WORLD > 2) || ADD_Z
777 grape_mesh->dimension_of_world = 3; /* !!!! */
778 #else
779 grape_mesh->dimension_of_world = 2;
780 #endif
781 #endif
782
783 grape_mesh->user_data = (void *)mesh;
784
785 element.mesh = (GRAPE_GENMESH *) grape_mesh;
786 #if FE_DIM==2
787 element.descr = &tria_description;
788 #else
789 element.descr = &tetra_description_even;
790 #endif
791
792 /* adjust max_vindex such that simple algorithm is ok. (dof[v][0]) */
793 for (i = 0; i < mesh->n_dof_admin; i++) {
794 if ((mesh->dof_admin[i]->n_dof[VERTEX] > 0) &&
795 (mesh->dof_admin[i]->n0_dof[VERTEX] == 0)) {
796 grape_mesh->max_vindex = mesh->dof_admin[i]->size_used;
797 MSG("use admin <%s> for vertex numbering, max=%d\n",
798 NAME(mesh->dof_admin[i]), mesh->dof_admin[i]->size_used);
799 break;
800 }
801 }
802
803 return grape_mesh;
804 }
805
806 /****************************************************************************/
807