/*--------------------------------------------------------------------------*/ /* ALBERTA: an Adaptive multi Level finite element toolbox using */ /* Bisectioning refinement and Error control by Residual */ /* Techniques for scientific Applications */ /* */ /* file: read_mesh.c */ /* */ /* description: reading data of mesh and vectors in machine-dependent */ /* binary format */ /* */ /*--------------------------------------------------------------------------*/ /* */ /* authors: Alfred Schmidt */ /* Zentrum fuer Technomathematik */ /* Fachbereich 3 Mathematik/Informatik */ /* Universitaet Bremen */ /* Bibliothekstr. 2 */ /* D-28359 Bremen, Germany */ /* */ /* Kunibert G. Siebert */ /* Institut fuer Mathematik */ /* Universitaet Augsburg */ /* Universitaetsstr. 14 */ /* D-86159 Augsburg, Germany */ /* */ /* http://www.alberta-fem.de/ */ /* */ /*--------------------------------------------------------------------------*/ /* (c) by A. Schmidt and K.G. Siebert (1996-2004) */ /* */ /* This program is free software; you can redistribute it and/or modify */ /* it under the terms of the GNU General Public License as published by */ /* the Free Software Foundation; either version 2 of the License, or */ /* any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /*--------------------------------------------------------------------------*/ /*--------------------------------------------------------------------------*/ /* to_do : number of dofs depending on node */ /* nach read_macro, damit dort pointer auf mel, v vorhanden sind!!! */ /* (fuers naechste Schreiben) */ /* error handling is terrible */ /*--------------------------------------------------------------------------*/ #include #include "alberta.h" static FILE *file = nil; static DOF_ADMIN *admin = nil; static MESH *mesh = nil; static U_CHAR preserve_coarse_dofs; static int n_vert_dofs = 0; static DOF **vert_dofs = nil; #if DIM > 1 static int n_edge_dofs = 0; static DOF **edge_dofs = nil; #endif #if DIM == 3 static int n_face_dofs = 0; static DOF **face_dofs = nil; #endif static EL *read_el_recursive(EL *parent); extern DOF_ADMIN *get_dof_admin(MESH *mesh, const char *name, const int n_dof[DIM+1]); /*--------------------------------------------------------------------------*/ static void read_dof_admins(MESH *mesh) { FUNCNAME("read_dof_admins"); int i, n_dof_admin, iadmin, used_count; int n_dof_el, n_dof[DIM+1], n_node_el, node[DIM+1]; int a_n_dof[DIM+1]; char *name; fread(&n_dof_el, sizeof(int), 1, file); fread(n_dof, sizeof(int), DIM+1, file); fread(&n_node_el, sizeof(int), 1, file); fread(node, sizeof(int), DIM+1, file); /* use data later for check */ fread(&n_dof_admin, sizeof(int), 1, file); for (iadmin = 0; iadmin < n_dof_admin; iadmin++) { fread(a_n_dof, sizeof(int), DIM+1, file); fread(&used_count, sizeof(int), 1, file); fread(&i, sizeof(int), 1, file); /* length without terminating \0 */ name = MEM_ALLOC(i+1, char); fread(name, sizeof(char), i, file); /* admin->name */ name[i] = 0; admin = get_dof_admin(mesh, name, a_n_dof); MEM_FREE(name, i+1, char); if (used_count > 0) enlarge_dof_lists(admin, used_count); } /* end for (iadmin) */ TEST(mesh->n_dof_el == n_dof_el)("wrong n_dof_el: %d %d\n", mesh->n_dof_el, n_dof_el); for (i=0; i<=DIM; i++) TEST(mesh->n_dof[i] == n_dof[i])("wrong n_dof[%d]: %d %d\n", i, mesh->n_dof[i], n_dof[i]); TEST(mesh->n_node_el == n_node_el)("wrong n_node_el: %d %d\n", mesh->n_node_el, n_node_el); for (i=0; i<=DIM; i++) TEST(mesh->node[i] == node[i])("wrong node[%d]: %d %d\n", i, mesh->node[i], node[i]); return; } /*--------------------------------------------------------------------------*/ MESH *read_mesh(const char *filename, REAL *timeptr, void (*init_leaf_data)(LEAF_DATA_INFO *), const BOUNDARY *(*init_boundary)(MESH *mesh, int ibdry)) { FUNCNAME("read_mesh"); MACRO_EL *mel; int i, j, n; REAL_D *v, x_min, x_max; int neigh_i[N_NEIGH]; char c99[99], *name; int size_REAL, iDIM, iDIM_OF_WORLD, ne, nv; REAL time, diam[DIM_OF_WORLD]; int vert_i[N_VERTICES]; static int funccount=0; #if DIM==3 S_CHAR edge_sc[N_EDGES]; #endif #if DIM > 1 S_CHAR bound_sc[N_NEIGH]; const BOUNDARY *(*init_bdry)(MESH *, int); init_bdry = (init_boundary) ? init_boundary : default_boundary; #endif TEST_EXIT(file=fopen(filename,"rb"))("cannot open file %s\n",filename); for (i=0; i < 98; i++) { fread(c99+i, sizeof(char), 1, file); if (c99[i] == 0) break; } if (i >= 98) { ERROR("problem reading id. abort.\n"); goto error_exit; } if (strncmp(c99, "ALBERTA", 6)) { c99[98] = 0; ERROR("unknown file id: \"%s\"\n",c99); goto error_exit; } fread(&size_REAL, sizeof(int), 1, file); if (size_REAL != sizeof(REAL)) { ERROR("wrong sizeof(REAL) %d. abort.\n", size_REAL); goto error_exit; } fread(&iDIM, sizeof(int), 1, file); if (iDIM != DIM) { ERROR("wrong DIM %d. abort.\n", iDIM); goto error_exit; } fread(&iDIM_OF_WORLD, sizeof(int), 1, file); if (iDIM_OF_WORLD != DIM_OF_WORLD) { ERROR("wrong DIM_OF_WORLD %d. abort.\n", iDIM_OF_WORLD); goto error_exit; } fread(&time, sizeof(REAL), 1, file); if (timeptr) *timeptr = time; fread(&i, sizeof(int), 1, file); /* length without terminating \0 */ if(i) { name = MEM_ALLOC(i+1, char); fread(name, sizeof(char), i, file); name[i] = 0; } else { funccount++; i=100; name = MEM_ALLOC(i+1, char); sprintf(name, "READ_MESH%d", funccount); } fread(diam, sizeof(REAL), DIM_OF_WORLD, file); fread(&preserve_coarse_dofs, sizeof(U_CHAR), 1, file); mesh = GET_MESH(name, read_dof_admins, init_leaf_data); MEM_FREE(name, i+1, char); for (i=0; idiam[i] = diam[i]; mesh->preserve_coarse_dofs = preserve_coarse_dofs; /* mesh->parametric = nil; */ fread(&n_vert_dofs, sizeof(int), 1, file); if (n_vert_dofs > 0) { vert_dofs = MEM_ALLOC(n_vert_dofs, DOF *); n = mesh->n_dof[VERTEX]; for (i = 0; i < n_vert_dofs; i++) { vert_dofs[i] = get_dof(mesh, VERTEX); fread(vert_dofs[i], sizeof(DOF), n, file); } } #if DIM > 1 fread(&n_edge_dofs, sizeof(int), 1, file); if (n_edge_dofs > 0) { edge_dofs = MEM_ALLOC(n_edge_dofs, DOF *); n = mesh->n_dof[EDGE]; for (i = 0; i < n_edge_dofs; i++) { edge_dofs[i] = get_dof(mesh, EDGE); fread(edge_dofs[i], sizeof(DOF), n, file); } } #endif #if DIM==3 fread(&n_face_dofs, sizeof(int), 1, file); if (n_face_dofs > 0) { face_dofs = MEM_ALLOC(n_face_dofs, DOF *); n = mesh->n_dof[FACE]; for (i = 0; i < n_face_dofs; i++) { face_dofs[i] = get_dof(mesh, FACE); fread(face_dofs[i], sizeof(DOF), n, file); } } #endif fread(&ne, sizeof(int), 1, file); fread(&nv, sizeof(int), 1, file); ((MESH_MEM_INFO *)mesh->mem_info)->count = nv; v = ((MESH_MEM_INFO *)mesh->mem_info)->coords = MEM_ALLOC(nv, REAL_D); fread(v, sizeof(REAL_D), nv, file); for (j = 0; j < DIM_OF_WORLD; j++) { x_min[j] = 1.E30; x_max[j] = -1.E30; } for (i = 0; i < nv; i++) for (j = 0; j < DIM_OF_WORLD; j++) { x_min[j] = MIN(x_min[j], v[i][j]); x_max[j] = MAX(x_max[j], v[i][j]); } for (j = 0; j < DIM_OF_WORLD; j++) mesh->diam[j] = x_max[j] - x_min[j]; mel = MEM_ALLOC(ne, MACRO_EL); for (i = 0; i < ne-1; i++) mel[i].next = mel+(i+1); for (i = 1; i < ne; i++) mel[i].last = mel+(i-1); mel[0].last = mel[ne-1].next = nil; mesh->n_macro_el = ne; mesh->first_macro_el = mel; for (n = 0; n < ne; n++) { mel[n].index = n; fread(vert_i, sizeof(int), N_VERTICES, file); for (i = 0; i < N_VERTICES; i++) { if ((vert_i[i] >= 0) && (vert_i[i] < nv)) mel[n].coord[i] = (REAL *)(v + (vert_i[i])); else mel[n].coord[i] = nil; } fread(mel[n].bound, sizeof(S_CHAR), N_VERTICES, file); #if DIM > 1 fread(bound_sc, sizeof(S_CHAR), N_NEIGH, file); for (i = 0; i < N_NEIGH; i++) { if (bound_sc[i]) mel[n].boundary[i] = init_bdry(mesh, bound_sc[i]); else mel[n].boundary[i] = nil; } #endif fread(neigh_i, sizeof(int), N_NEIGH, file); for (i = 0; i < N_NEIGH; i++) { if ((neigh_i[i] >= 0) && (neigh_i[i] < ne)) mel[n].neigh[i] = mel + (neigh_i[i]); else mel[n].neigh[i] = nil; } fread(mel[n].opp_vertex, sizeof(U_CHAR), N_NEIGH, file); #if DIM==3 fread(edge_sc, sizeof(S_CHAR), N_EDGES, file); for (i = 0; i < N_EDGES; i++) { #if 1 if (edge_sc[i]) mel[n].boundary[N_FACES+i] = init_bdry(mesh, edge_sc[i]); else mel[n].boundary[N_FACES+i] = nil; #else if (edge_sc[i] > 0) mel[n].boundary[N_FACES+i] = bound+1; else if (edge_sc[i] < 0) mel[n].boundary[N_FACES+i] = bound+2; else mel[n].boundary[N_FACES+i] = nil; #endif } fread(&(mel[n].el_type), sizeof(U_CHAR), 1, file); /* read el_type */ { extern S_CHAR get_orientation(MACRO_EL *mel); mel[n].orientation = get_orientation(&mel[n]); } #endif mel[n].el = read_el_recursive(nil); } fread(&mesh->n_vertices, sizeof(int), 1, file); if (!strncmp((char *)(&mesh->n_vertices), "EOF.", 4)) /* file end marker ? */ { MSG("Warning: old mesh file format.\n"); mesh->n_vertices = n_vert_dofs; goto error_exit; } if (mesh->n_vertices != n_vert_dofs) { ERROR("n_vertices != n_vert_dofs.\n"); mesh->n_vertices = n_vert_dofs; goto error_exit; } #if DIM > 1 fread(&mesh->n_edges, sizeof(int), 1, file); #endif fread(&mesh->n_elements, sizeof(int), 1, file); fread(&mesh->n_hier_elements, sizeof(int), 1, file); #if DIM == 3 fread(&mesh->n_faces, sizeof(int), 1, file); fread(&mesh->max_edge_neigh, sizeof(int), 1, file); #endif if (fread(c99, sizeof(char), 4, file) != 4) { ERROR("problem reading FILE END MARK.\n"); goto error_exit; } if (strncmp(c99, "EOF.", 4)) /* file end marker */ { ERROR("no FILE END MARK.\n"); goto error_exit; } #if NEIGH_IN_EL ERROR_EXIT("read_mesh for NEIGH_IN_EL=1 not implemented yet!!!\n"); #endif error_exit: fclose(file); file = nil; #if DIM == 3 if(n_face_dofs) MEM_FREE(face_dofs, n_face_dofs, DOF *); #endif #if DIM > 1 if(n_edge_dofs) MEM_FREE(edge_dofs, n_edge_dofs, DOF *); #endif if(n_vert_dofs) MEM_FREE(vert_dofs, n_vert_dofs, DOF *); return(mesh); } /*--------------------------------------------------------------------------*/ static EL *read_el_recursive(EL *parent) { FUNCNAME("read_el_recursive"); int i, j, n, node0; EL *el; U_CHAR uc; #if DIM > 1 U_CHAR nc; #endif el = get_element(mesh); mesh->n_hier_elements++; #if EL_INDEX el->index = mesh->n_hier_elements; #endif fread(&uc, sizeof(U_CHAR), 1, file); #if DIM > 1 fread(&nc, sizeof(U_CHAR), 1, file); if (nc) { el->new_coord = get_real_d(mesh); fread(el->new_coord, sizeof(REAL_D), 1, file); } else { el->new_coord = nil; } #endif if (mesh->n_dof[VERTEX] > 0) { node0 = mesh->node[VERTEX]; for (i = 0; i < N_VERTICES; i++) { fread(&j, sizeof(int), 1, file); TEST_EXIT(j < n_vert_dofs) ("vert_dofs index too large: %d >= %d\n", j, n_vert_dofs); el->dof[node0 + i] = vert_dofs[j]; } } if ((!uc) || preserve_coarse_dofs) { #if DIM > 1 if (mesh->n_dof[EDGE] > 0) { node0 = mesh->node[EDGE]; for (i = 0; i < N_EDGES; i++) { fread(&j, sizeof(int), 1, file); TEST_EXIT(j < n_edge_dofs) ("edge_dofs index too large: %d >= %d\n", j, n_edge_dofs); el->dof[node0 + i] = edge_dofs[j]; } } #endif #if DIM==3 if ((n = mesh->n_dof[FACE]) > 0) { node0 = mesh->node[FACE]; for (i = 0; i < N_FACES; i++) { fread(&j, sizeof(int), 1, file); TEST_EXIT(j < n_face_dofs) ("face_dofs index too large: %d >= %d\n", j, n_face_dofs); el->dof[node0 + i] = face_dofs[j]; } } #endif if ((n = mesh->n_dof[CENTER]) > 0) { node0 = mesh->node[CENTER]; el->dof[node0] = get_dof(mesh, CENTER); fread(el->dof[node0], sizeof(DOF), n, file); } } #if NEIGH_IN_EL for (i = 0; i < N_NEIGH; i++) { el->neigh[i] = nil; el->opp_vertex[i] = 0; } #endif if (uc) { el->child[0] = read_el_recursive(el); el->child[1] = read_el_recursive(el); } else mesh->n_elements++; return(el); } #if 0 /*--------------------------------------------------------------------------*/ /* read DOF vectors of various types */ /*--------------------------------------------------------------------------*/ typedef DOF_REAL_VEC DOF_VEC; static const DOF_ADMIN *read_dof_vec(const char *filename, DOF_VEC *dv, int unit_size, const char dofvectype[16], MESH *mesh, FE_SPACE *fe_space) { FUNCNAME("read_dof_vec"); int i, last; int n_dof[DIM+1]; const DOF_ADMIN *admin = nil; const BAS_FCTS *bas_fcts = nil; char c20[20], *name; TEST_EXIT(mesh)("no mesh given\n"); TEST_EXIT(file=fopen(filename,"rb"))("cannot open file %s\n",filename); if (fread(c20, sizeof(char), 16, file) != 16) { ERROR("cannot read file id\n"); goto error_exit; } if (strncmp(c20, dofvectype, 12)) { c20[16] = 0; ERROR("invalid file id; %s\n", c20); goto error_exit; } fread(&last, sizeof(int), 1, file); /* length of DOF_VEC name */ if (last) { name = MEM_ALLOC(last+1, char); fread(name, sizeof(char), last, file); name[last] = 0; } else name = strdup(filename); dv->name = name; fread(n_dof, sizeof(int), DIM+1, file); /* for DOF_ADMIN selection */ fread(&last, sizeof(int), 1, file); /* length of BAS_FCTS name */ if (last) { name = MEM_ALLOC(last+1, char); fread(name, sizeof(char), last, file); name[last] = 0; if (fe_space && (bas_fcts = fe_space->bas_fcts)) { if (strcmp(bas_fcts->name, name)) { ERROR("invalid name %s is not given fe_space->bas_fcts->name %s\n", name, bas_fcts->name); } } else { /* no given fe_space or no bas_fcts in given fe_space */ TEST_EXIT(bas_fcts = get_bas_fcts(name)) ("cannot get bas_fcts <%s>\n", name); if (fe_space) { /* use given fe_space */ fe_space->bas_fcts = bas_fcts; } else { /* create new fe_space */ TEST_EXIT(fe_space = (FE_SPACE *)get_fe_space(mesh, name, n_dof, bas_fcts)) ("cannot get fe_space for bas_fcts <%s>\n", name); } } for (i=0; i<=DIM; i++) { TEST_EXIT(n_dof[i] == bas_fcts->n_dof[i]) ("wrong n_dof in bas_fcts <%s>", name); } } else { /* no bas_fcts.name in file */ if (fe_space) { /* use given fe_space */ TEST_EXIT(admin = fe_space->admin)("no fe_space->admin"); for (i=0; i<=DIM; i++) { TEST_EXIT(n_dof[i] == admin->n_dof[i]) ("wrong n_dof in admin <%s>", NAME(admin)); } } else { /* create new fe_space */ TEST_EXIT(fe_space = (FE_SPACE *) get_fe_space(mesh, nil, n_dof, nil)) ("cannot get fe_space for given n_dof\n"); TEST_EXIT(admin = fe_space->admin)("no admin in new fe_space\n"); for (i=0; i<=DIM; i++) { TEST_EXIT(n_dof[i] == admin->n_dof[i]) ("wrong n_dof in admin <%s>", NAME(admin)); } } MEM_FREE(name, last+1, char); } TEST_EXIT(fe_space)("still no fe_space\n"); dv->fe_space = fe_space; TEST_EXIT(admin = fe_space->admin)("still no admin\n"); dof_compress(mesh); fread(&last, sizeof(int), 1, file); if (last != admin->size_used) { ERROR("size of dof vector `%s' does not fit to size_used in admin `%s'\n", dv->name, admin->name); ERROR_EXIT("can not read incompatible data\n"); } if (last) { dv->size = last; dv->vec = ((REAL *) alberta_alloc((size_t)(last*unit_size), funcName, __FILE__, __LINE__)); fread(dv->vec, unit_size, last, file); } else { ERROR("empty dof vector\n"); dv->size = 0; dv->vec = nil; } if (fread(c20, sizeof(char), 4, file) != 4) { print_error_msg("ERROR: problem reading file end mark.\n"); goto error_exit; } if (strncmp(c20, "EOF.", 4)) /* file end marker */ { print_error_msg("ERROR: no file end mark.\n"); goto error_exit; } error_exit: fclose(file); file = nil; return(admin); } /*--------------------------------------------------------------------------*/ DOF_REAL_VEC *read_dof_real_vec(const char *filename, MESH *mesh, FE_SPACE *fe_space) { DOF_REAL_VEC *dv; const DOF_ADMIN *admin; void add_dof_real_vec_to_admin(DOF_REAL_VEC *dv, DOF_ADMIN *admin); dv = get_dof_real_vec(filename, nil); admin = read_dof_vec(filename, (DOF_VEC *)dv, sizeof(REAL), "DOF_REAL_VEC ", mesh, fe_space); if (admin) add_dof_real_vec_to_admin(dv, (DOF_ADMIN *)admin); return(dv); } /*--------------------------------------------------------------------------*/ DOF_REAL_D_VEC *read_dof_real_d_vec(const char *filename, MESH *mesh, FE_SPACE *fe_space) { DOF_REAL_D_VEC *dv; const DOF_ADMIN *admin; void add_dof_real_d_vec_to_admin(DOF_REAL_D_VEC *dv, DOF_ADMIN *admin); dv = get_dof_real_d_vec(filename, nil); admin = read_dof_vec(filename, (DOF_VEC *)dv, sizeof(REAL_D), "DOF_REAL_D_VEC ", mesh, fe_space); if (admin) add_dof_real_d_vec_to_admin(dv, (DOF_ADMIN *)admin); return(dv); } /*--------------------------------------------------------------------------*/ DOF_INT_VEC *read_dof_int_vec(const char *filename, MESH *mesh, FE_SPACE *fe_space) { DOF_INT_VEC *dv; const DOF_ADMIN *admin; void add_dof_int_vec_to_admin(DOF_INT_VEC *dv, DOF_ADMIN *admin); dv = get_dof_int_vec(filename, nil); admin = read_dof_vec(filename, (DOF_VEC *)dv, sizeof(int), "DOF_INT_VEC ", mesh, fe_space); if (admin) add_dof_int_vec_to_admin(dv, (DOF_ADMIN *)admin); return(dv); } /*--------------------------------------------------------------------------*/ DOF_SCHAR_VEC *read_dof_schar_vec(const char *filename, MESH *mesh, FE_SPACE *fe_space) { DOF_SCHAR_VEC *dv; const DOF_ADMIN *admin; void add_dof_schar_vec_to_admin(DOF_SCHAR_VEC *dv, DOF_ADMIN *admin); dv = get_dof_schar_vec(filename, nil); admin = read_dof_vec(filename, (DOF_VEC *)dv, sizeof(S_CHAR), "DOF_SCHAR_VEC ", mesh, fe_space); if (admin) add_dof_schar_vec_to_admin(dv, (DOF_ADMIN *)admin); return(dv); } /*--------------------------------------------------------------------------*/ DOF_UCHAR_VEC *read_dof_uchar_vec(const char *filename, MESH *mesh, FE_SPACE *fe_space) { DOF_UCHAR_VEC *dv; const DOF_ADMIN *admin; void add_dof_uchar_vec_to_admin(DOF_UCHAR_VEC *dv, DOF_ADMIN *admin); dv = get_dof_uchar_vec(filename, nil); admin = read_dof_vec(filename, (DOF_VEC *)dv, sizeof(U_CHAR), "DOF_UCHAR_VEC ", mesh, fe_space); if (admin) add_dof_uchar_vec_to_admin(dv, (DOF_ADMIN *)admin); return(dv); } /*--------------------------------------------------------------------------*/ #endif