source: src/par/MESH_LabelPType.c @ 689:4235fdcc79c8

Revision 689:4235fdcc79c8, 7.1 KB checked in by duowang, 2 years ago (diff)

fixed the bug: two partitions touch at a corner but has no adjancy information

Line 
1#include <stdio.h>
2#include <stdlib.h>
3#include <string.h>
4#include <time.h>
5#include "MSTK.h"
6#include "MSTK_private.h"
7#ifdef __cplusplus
8extern "C" {
9#endif
10
11
12
13  /*
14     This function is a collective call
15     It assign vertex PType similar as MESH_AssignGlobalIDs(),
16     except that here global IDs are given
17
18     It assigns all the elements with a POVERLAP vertex as POVERLAP element
19
20     It also assign master partition ID
21
22     Author(s): Duo Wang, Rao Garimella
23  */
24
25  int MESH_LabelPType_Vertex(Mesh_ptr submesh, int rank, int num, MPI_Comm comm);
26  int MESH_LabelPType_Face(Mesh_ptr submesh, int rank, int num, MPI_Comm comm);
27  int MESH_LabelPType_Region(Mesh_ptr submesh, int rank, int num, MPI_Comm comm);
28
29int MESH_LabelPType(Mesh_ptr submesh, int rank, int num,  MPI_Comm comm) {
30  int nf, nr;
31  RepType rtype;
32
33  /* basic mesh information */
34  rtype = MESH_RepType(submesh);
35  nf = MESH_Num_Faces(submesh);
36  nr = MESH_Num_Regions(submesh);
37  /* first label vertex */
38  MESH_LabelPType_Vertex(submesh, rank, num, comm);
39  if (nr)
40    MESH_LabelPType_Region(submesh, rank, num, comm);
41  else if(nf) 
42    MESH_LabelPType_Face(submesh, rank, num, comm);
43  else {
44    MSTK_Report("MESH_LabelPType()","only send volume or surface mesh",MSTK_ERROR);
45    exit(-1);
46  }
47  return 1;
48}
49
50int MESH_LabelPType_Vertex(Mesh_ptr submesh, int rank, int num, MPI_Comm comm) {
51  int i, j, k, nv, nbv, ne, nf, mesh_info[10], nevs, nfes, nfv, natt, nset, ncomp, dir;
52  int nfe;
53  MVertex_ptr mv;
54  MEdge_ptr me;
55  MFace_ptr mf;
56  List_ptr boundary_verts, fverts, mfedges;
57  RepType rtype;
58  char attname[256], msetname[256];
59  MType mtype;
60  MAttType att_type;
61  MAttrib_ptr attrib;
62  MSet_ptr mset;
63  int *list_attr, *list_mset;
64  char *list_attr_names, *list_mset_names;
65  int *loc;
66  int iloc, num_ghost_verts, global_id;
67  double coor[3];
68  for (i = 0; i < 10; i++) mesh_info[i] = 0;
69
70  /* mesh_info store the mesh reptype, nv, ne, nf, nbv */
71  rtype = MESH_RepType(submesh);
72  nv = MESH_Num_Vertices(submesh);
73  ne = MESH_Num_Edges(submesh);
74  nf = MESH_Num_Faces(submesh);
75
76  mesh_info[0] = rtype;
77  mesh_info[1] = nv;
78  mesh_info[2] = ne;
79  mesh_info[3] = nf;
80
81  /* calculate number of boundary vertices */ 
82  nbv = 0;
83  boundary_verts = List_New(10);
84  for(i = 0; i < nv; i++) {
85    mv = MESH_Vertex(submesh,i);
86    if (MV_PType(mv)) {
87      List_Add(boundary_verts,mv);
88      nbv++;
89    }
90  }
91  mesh_info[4] = nbv;
92 
93  /* sort boundary vertices based on global ID, for binary search */
94  List_Sort(boundary_verts,nbv,sizeof(MVertex_ptr),compareGlobalID);
95
96  /*
97     gather submeshes information
98     right now we only need nv and nbv, and later num_ghost_verts, but we gather all mesh_info
99  */
100  int *global_mesh_info = (int *)MSTK_malloc(10*num*sizeof(int));
101  MPI_Allgather(mesh_info,10,MPI_INT,global_mesh_info,10,MPI_INT,comm);
102
103  /* get largest number of boundary vertices of all the processors */
104  int max_nbv = 0;
105  for(i = 0; i < num; i++)
106    if(max_nbv < global_mesh_info[10*i+4])
107      max_nbv = global_mesh_info[10*i+4];
108
109  int *list_boundary_vertex = (int *)MSTK_malloc(max_nbv*sizeof(int));
110  int *recv_list_vertex = (int *)MSTK_malloc(num*max_nbv*sizeof(int));
111
112  /* only global ID are sent */
113  int index_nbv = 0;
114  for(i = 0; i < nbv; i++) {
115    mv = List_Entry(boundary_verts,i);
116    list_boundary_vertex[index_nbv] = MV_GlobalID(mv);
117    index_nbv++;
118  }
119
120  /* gather boundary vertices */
121  MPI_Allgather(list_boundary_vertex,max_nbv,MPI_INT,recv_list_vertex,max_nbv,MPI_INT,comm);
122
123  /* indicate if a vertex is overlapped */
124  int *vertex_ov_label = (int *)MSTK_malloc(num*max_nbv*sizeof(int));
125
126  /*
127     store the local boundary id on ov processor
128     it is used to assign global id of local ghost vertices
129     do not need to store master partition id, MV_MasterParID(mv) is assigned
130  */
131  for (i = 0; i < num*max_nbv; i++)
132    vertex_ov_label[i] = 0;
133  num_ghost_verts = 0;
134  /* for processor other than 0 */
135  for(i = 0; i < nbv; i++) {
136    mv = List_Entry(boundary_verts,i);
137    global_id = MV_GlobalID(mv);
138    /* check which processor has the same coordinate vertex
139     Different from assigning global id, we check all the other processors, from large to small
140     by rank, whenever a vertex is found, mark rank and j as neighbors, and the masterparid is still
141     the smallest rank processor
142    */
143    for(j = num-1; j >= 0; j--) {
144      if(j == rank) continue;
145      /* since each processor has sorted the boundary vertices, use binary search */
146      loc = (int *)bsearch(&global_id,
147                           &recv_list_vertex[max_nbv*j],
148                           global_mesh_info[10*j+4],
149                           sizeof(int),
150                           compareINT);
151      /* if found the vertex on previous processors */
152      if(loc) {
153        /* here the location iloc is relative to the beginning of the jth processor */
154        iloc = (int)(loc - &recv_list_vertex[max_nbv*j]);
155        MV_Set_PType(mv,PGHOST);
156        MV_Set_MasterParID(mv,j);
157        MESH_Flag_Has_Ghosts_From_Prtn(submesh,j,MVERTEX);
158        num_ghost_verts++;
159        /* label the original vertex as overlapped */
160        vertex_ov_label[max_nbv*j+iloc] |= 1;
161        /* if found on processor j, no need to test other processors*/
162        break;
163      }
164    }
165  }
166
167  mesh_info[5] = num_ghost_verts;
168
169  /* update ghost verts number */
170  MPI_Allgather(mesh_info,10,MPI_INT,global_mesh_info,10,MPI_INT,comm);
171  /* since this is a or reduction, we can use MPI_IN_PLACE, send buffer same as recv buffer */
172  MPI_Allreduce(MPI_IN_PLACE,vertex_ov_label,num*max_nbv,MPI_INT,MPI_LOR,comm);   
173
174  for(i = 0; i < nv; i++) {
175    mv = MESH_Vertex(submesh,i);
176    if (MV_PType(mv) == PGHOST)
177      continue;
178    MV_Set_MasterParID(mv,rank);
179  }
180
181  /* label OVERLAP vertices */
182  for(i = 0; i < nbv; i++) {
183    if(vertex_ov_label[rank*max_nbv+i]) {
184      mv = List_Entry(boundary_verts,i);
185      MV_Set_PType(mv,POVERLAP);
186    }
187  }
188
189  List_Delete(boundary_verts);
190  MSTK_free(global_mesh_info);
191  MSTK_free(vertex_ov_label);
192
193  MSTK_free(list_boundary_vertex);
194  MSTK_free(recv_list_vertex);
195  return 1;
196}
197 
198  /*
199     Right now assume faces are not overlapped across processors
200     Label the face that has OVERLAP or Ghost vertex as OVERLAP
201     
202  */
203int MESH_LabelPType_Face(Mesh_ptr submesh, int rank, int num, MPI_Comm comm) {
204  int i, idx;
205  MVertex_ptr mv;
206  MFace_ptr mf;
207  List_ptr fverts;
208  idx = 0;
209  while( (mf = MESH_Next_Face(submesh,&idx)) ) {
210      fverts = MF_Vertices(mf,1,0);
211      for(i = 0; i < List_Num_Entries(fverts); i++) {
212        mv = List_Entry(fverts,i);
213        if( MV_PType(mv) == PGHOST || MV_PType(mv) == POVERLAP ) {
214          MF_Set_PType(mf,POVERLAP);
215          break;
216        }
217      }
218      List_Delete(fverts);
219    }
220   
221  return 1;
222}
223
224  /*
225     right now assume regions are not overlapped across processors
226     Label the region that has OVERLAP vertex as OVERLAP
227  */
228
229
230int MESH_LabelPType_Region(Mesh_ptr submesh, int rank, int num, MPI_Comm comm) {
231  int i, idx;
232  MRegion_ptr mr;
233  MVertex_ptr mv;
234  List_ptr rverts;
235  idx = 0;
236  while( (mr=MESH_Next_Region(submesh,&idx)) ) {
237      rverts = MR_Vertices(mr);
238      for(i = 0; i < List_Num_Entries(rverts); i++) {
239        mv = List_Entry(rverts,i);
240        if(MV_PType(mv) == PGHOST || MV_PType(mv) == POVERLAP) {
241          MR_Set_PType(mr,POVERLAP);
242          break;
243        }
244      }
245      List_Delete(rverts);
246    }
247   
248  return 1;
249}
250
251
252 
253#ifdef __cplusplus
254}
255#endif
256
Note: See TracBrowser for help on using the repository browser.