source: issm/trunk-jpl/src/c/modules/ModelProcessorx/ElementsAndVerticesPartitioning.cpp@ 23505

Last change on this file since 23505 was 23505, checked in by Mathieu Morlighem, 6 years ago

CHG: converting my_vertices back to bool
File size: 4.9 KB
Line 
1/*!\file: ElementsAndVerticesPartitioning.cpp
2 * \brief: partition elements and nodes and vertices
3 */
4
5#ifdef HAVE_CONFIG_H
6 #include <config.h>
7#else
8#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
9#endif
10
11#include <string.h>
12#include "../../classes/classes.h"
13#include "../../shared/shared.h"
14#include "../MeshPartitionx/MeshPartitionx.h"
15#include "../ModelProcessorx/ModelProcessorx.h"
16
17void ElementsAndVerticesPartitioning(bool** pmy_elements,bool** pmy_vertices,IoModel* iomodel){
18
19 int numberofelements2d;
20 int numberofvertices2d;
21 int numlayers;
22 int numvertex_pairing;
23
24 /*intermediary: */
25 int *epart = NULL; //element partitioning.
26 int *npart = NULL; //node partitioning.
27 int elements_width; //number of columns in elements (2d->3, 3d->6)
28 int *elements2d = NULL;
29
30 /*Get my_rank:*/
31 int my_rank = IssmComm::GetRank();
32 int num_procs = IssmComm::GetSize();
33
34 /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
35 if(*pmy_elements) return;
36
37 /*Number of vertices per elements, needed to correctly retrieve data: */
38 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
39 switch(iomodel->meshelementtype){
40 case TriaEnum:
41 elements_width=3;
42 numberofelements2d = 0;
43 numberofvertices2d = 0;
44 numlayers = 0;
45 break;
46 case TetraEnum:
47 elements_width=4;
48 numberofelements2d = 0;
49 numberofvertices2d = 0;
50 numlayers = 0;
51 break;
52 case PentaEnum:
53 elements_width=6;
54 iomodel->FetchData(&elements2d,NULL,NULL,"md.mesh.elements2d");
55 iomodel->FindConstant(&numberofelements2d,"md.mesh.numberofelements2d");
56 iomodel->FindConstant(&numberofvertices2d,"md.mesh.numberofvertices2d");
57 iomodel->FindConstant(&numlayers,"md.mesh.numberoflayers");
58 break;
59 default:
60 _error_("mesh elements "<< EnumToStringx(iomodel->meshelementtype) <<" not supported yet");
61 }
62
63 /*Partition and free resouces*/
64 MeshPartitionx(&epart,&npart,iomodel->numberofelements,iomodel->numberofvertices,iomodel->elements,numberofelements2d,numberofvertices2d,elements2d,numlayers,elements_width,iomodel->meshelementtype,num_procs);
65 xDelete<int>(elements2d);
66 xDelete<int>(npart);
67
68 /*Deal with rifts, they have to be included into one partition only, not several: */
69 int numrifts;
70 iomodel->FindConstant(&numrifts,"md.rifts.numrifts");
71 if(numrifts){
72 IssmDouble *riftinfo = NULL;
73 iomodel->FetchData(&riftinfo,&numrifts,NULL,"md.rifts.riftstruct");
74 for(int i=0;i<numrifts;i++){
75 const int RIFTINFOSIZE = 12;
76 int el1=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+2))-1; //matlab indexing to c indexing
77 int el2=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+3))-1; //matlab indexing to c indexing
78 epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding vertices;
79 }
80 iomodel->DeleteData(riftinfo,"md.rifts.riftstruct");
81 }
82
83 /*Create my_vertices and my_elements, used by each partition */
84 bool *my_elements = xNewZeroInit<bool>(iomodel->numberofelements);
85 bool *my_vertices = xNewZeroInit<bool>(iomodel->numberofvertices);
86
87 /*Start figuring out, out of the partition, which elements belong to this cpu: */
88 for(int i=0;i<iomodel->numberofelements;i++){
89
90 /*!All elements have been partitioned above, only deal with elements for this cpu: */
91 if(my_rank==epart[i]){
92 my_elements[i]=true;
93 /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
94 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
95 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
96 will hold which vertices belong to this partition*/
97 for(int j=0;j<elements_width;j++){
98 my_vertices[iomodel->elements[elements_width*i+j]-1]=true;
99 }
100 }
101 }
102
103 /*We might have vertex_pairing in which case, some vertices have to be cloned:
104 * penpair has 2 nodes that are poointing toward 2 vertices.
105 * The 2 vertices must be in the same cpu as the penpair*/
106 int *vertex_pairing = NULL;
107 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.stressbalance.vertex_pairing");
108 for(int i=0;i<numvertex_pairing;i++){
109 if(my_vertices[vertex_pairing[2*i+0]-1] && !my_vertices[vertex_pairing[2*i+1]-1]){
110 my_vertices[vertex_pairing[2*i+1]-1]=true;
111 }
112 }
113 xDelete<int>(vertex_pairing);
114 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,"md.masstransport.vertex_pairing");
115 for(int i=0;i<numvertex_pairing;i++){
116 if(my_vertices[vertex_pairing[2*i+0]-1] && !my_vertices[vertex_pairing[2*i+1]-1]){
117 my_vertices[vertex_pairing[2*i+1]-1]=true;
118 }
119 }
120 xDelete<int>(vertex_pairing);
121
122 /*cleanup and assign output pointer*/
123 xDelete<int>(epart);
124 *pmy_elements=my_elements;
125 *pmy_vertices=my_vertices;
126}
Note: See TracBrowser for help on using the repository browser.