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

Last change on this file since 17472 was 17472, checked in by Mathieu Morlighem, 12 years ago

NEW: working on Tetra elements
File size: 5.4 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, int** pmy_vertices, IoModel* iomodel){
18
19 int i;
20
21 const int RIFTINFOSIZE = 12;
22 int numberofelements2d;
23 int numberofvertices2d;
24 int numlayers;
25 int numrifts;
26 int numvertex_pairing;
27
28 /*output: */
29 bool *my_elements = NULL;
30 int *my_vertices = NULL;
31
32 /*intermediary: */
33 int *epart = NULL; //element partitioning.
34 int *npart = NULL; //node partitioning.
35 int elements_width; //number of columns in elements (2d->3, 3d->6)
36 int el1,el2;
37 int *elements2d = NULL;
38 int *vertex_pairing = NULL;
39 IssmDouble *riftinfo = NULL;
40
41 /*Get my_rank:*/
42 int my_rank = IssmComm::GetRank();
43 int num_procs = IssmComm::GetSize();
44
45 /*Fetch parameters: */
46 iomodel->Constant(&numrifts,RiftsNumriftsEnum);
47
48 /*First, check that partitioning has not yet been carryed out. Just check whether my_elements pointers is not already assigned a value: */
49 if(*pmy_elements)return;
50
51 /*Number of vertices per elements, needed to correctly retrieve data: */
52 /*Determine parallel partitioning of elements: we use Metis for now. First load the data, then partition*/
53 switch(iomodel->meshtype){
54 case Mesh2DhorizontalEnum:
55 elements_width=3;
56 numberofelements2d = 0;
57 numberofvertices2d = 0;
58 numlayers = 0;
59 break;
60 case Mesh2DverticalEnum:
61 elements_width=3;
62 numberofelements2d = 0;
63 numberofvertices2d = 0;
64 numlayers = 0;
65 break;
66 case Mesh3DtetrasEnum:
67 elements_width=4; //penta elements
68 numberofelements2d = 0;
69 numberofvertices2d = 0;
70 numlayers = 0;
71 break;
72 case Mesh3DEnum:
73 elements_width=6; //penta elements
74 iomodel->FetchData(&elements2d,NULL,NULL,MeshElements2dEnum);
75 iomodel->Constant(&numberofelements2d,MeshNumberofelements2dEnum);
76 iomodel->Constant(&numberofvertices2d,MeshNumberofvertices2dEnum);
77 iomodel->Constant(&numlayers,MeshNumberoflayersEnum);
78 break;
79 default:
80 _error_("mesh not supported yet");
81 }
82
83 MeshPartitionx(&epart,&npart,iomodel->numberofelements,iomodel->numberofvertices,iomodel->elements,numberofelements2d,numberofvertices2d,elements2d,numlayers,elements_width,iomodel->meshtype,num_procs);
84
85 /*Free elements2d: */
86 xDelete<int>(elements2d);
87
88 /*Deal with rifts, they have to be included into one partition only, not several: */
89 if(numrifts){
90 iomodel->FetchData(&riftinfo,&numrifts,NULL,RiftsRiftstructEnum);
91 for(i=0;i<numrifts;i++){
92 el1=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+2))-1; //matlab indexing to c indexing
93 el2=reCast<int>(*(riftinfo+RIFTINFOSIZE*i+3))-1; //matlab indexing to c indexing
94 epart[el2]=epart[el1]; //ensures that this pair of elements will be in the same partition, as well as the corresponding vertices;
95 }
96 iomodel->DeleteData(riftinfo,RiftsRiftstructEnum);
97 }
98
99 /*Used later on: */
100 my_vertices=xNewZeroInit<int>(iomodel->numberofvertices);
101 my_elements=xNewZeroInit<bool>(iomodel->numberofelements);
102
103 /*Start figuring out, out of the partition, which elements belong to this cpu: */
104 for (i=0;i<iomodel->numberofelements;i++){
105
106 /*!All elements have been partitioned above, only deal with elements for this cpu: */
107 if(my_rank==epart[i]){
108
109 my_elements[i]=true;
110
111 /*Now that we are here, we can also start building the list of vertices belonging to this cpu partition: we use
112 *the element index to do this. For each element n, we know index[n][0:2] holds the indices (matlab indexing)
113 into the vertices coordinates. If we start plugging 1 into my_vertices for each index[n][i] (i=0:2), then my_vertices
114 will hold which vertices belong to this partition*/
115 my_vertices[iomodel->elements[elements_width*i+0]-1]=1;
116 my_vertices[iomodel->elements[elements_width*i+1]-1]=1;
117 my_vertices[iomodel->elements[elements_width*i+2]-1]=1;
118
119 if(elements_width==6){
120 my_vertices[iomodel->elements[elements_width*i+3]-1]=1;
121 my_vertices[iomodel->elements[elements_width*i+4]-1]=1;
122 my_vertices[iomodel->elements[elements_width*i+5]-1]=1;
123 }
124 }
125 }
126
127 /*We might have vertex_pairing in which case, some vertices have to be cloned:
128 * penpair has 2 nodes that are poointing toward 2 vertices.
129 * The 2 vertices must be in the same cpu as the penpair*/
130 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,StressbalanceVertexPairingEnum);
131 for(i=0;i<numvertex_pairing;i++){
132 if(my_vertices[vertex_pairing[2*i+0]-1] && !my_vertices[vertex_pairing[2*i+1]-1]){
133 my_vertices[vertex_pairing[2*i+1]-1]=2; //to know that these elements are not on the partition
134 }
135 }
136 xDelete<int>(vertex_pairing);
137 iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
138 for(i=0;i<numvertex_pairing;i++){
139 if(my_vertices[vertex_pairing[2*i+0]-1] && !my_vertices[vertex_pairing[2*i+1]-1]){
140 my_vertices[vertex_pairing[2*i+1]-1]=2; //to know that these elements are not on the partition
141 }
142 }
143 xDelete<int>(vertex_pairing);
144
145 /*Free ressources:*/
146 xDelete<int>(npart);
147 xDelete<int>(epart);
148
149 /*Assign output pointers:*/
150 *pmy_elements=my_elements;
151 *pmy_vertices=my_vertices;
152}
Note: See TracBrowser for help on using the repository browser.