196dfb0cf81cfad0c91361a6e711d7fc0d5bd367
[charm.git] / src / libs / ck-libs / ParFUM-Tops / ParFUM_TOPS.C
1 /**
2
3    @file
4    @brief Implementation of ParFUM-TOPS layer, except for Iterators
5
6    @author Isaac Dooley
7
8    @todo add code to generate ghost layers
9    @todo Support multiple models
10
11    @note FEM_DATA+0 holds the elemAttr or nodeAtt data
12    @note FEM_DATA+1 holds the id
13    @note FEM_CONN holds the element connectivity
14    @note FEM_COORD holds the nodal coordinates
15
16 */
17
18 #include "ParFUM_TOPS.h"
19 #include "ParFUM.decl.h"
20 #include "ParFUM_internals.h"
21
22 TopModel* topModel_Create_Init(int elem_attr_sz, int node_attr_sz){
23   CkAssert(elem_attr_sz > 0);
24   CkAssert(node_attr_sz > 0);
25   TopModel *model = new TopModel;
26   model->elem_attr_size = elem_attr_sz;
27   model->node_attr_size = node_attr_sz;
28
29   // This only uses a single mesh, so better not create multiple ones of these
30   int which_mesh=FEM_Mesh_default_write();
31   model->mesh = FEM_Mesh_lookup(which_mesh,"TopModel::TopModel()");
32
33   char* temp_array = new char[16]; // just some junk array to use below
34
35   // Allocate node coords
36 #ifdef FP_TYPE_FLOAT
37   FEM_Mesh_data(which_mesh,FEM_NODE,FEM_COORD,temp_array, 0, 0, FEM_FLOAT, 3);
38 #else
39   FEM_Mesh_data(which_mesh,FEM_NODE,FEM_COORD,temp_array, 0, 0, FEM_DOUBLE, 3);
40 #endif
41
42   // Allocate element connectivity
43   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_CONN,temp_array, 0, 0, FEM_INDEX_0, 4);
44   // Allocate element attributes
45   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_DATA+0,temp_array, 0, 0, FEM_BYTE, model->elem_attr_size);
46   // Allocate element Id array
47   FEM_Mesh_data(which_mesh,FEM_ELEM+0,FEM_DATA+1,temp_array, 0, 0, FEM_INT, 1);
48
49
50
51   // Allocate node attributes
52   FEM_Mesh_data(which_mesh,FEM_NODE+0,FEM_DATA+0,temp_array, 0, 0, FEM_BYTE, model->node_attr_size);
53   // Allocate node Id array
54   FEM_Mesh_data(which_mesh,FEM_NODE+0,FEM_DATA+1,temp_array, 0, 0, FEM_INT, 1);
55
56   delete[] temp_array;
57
58   // Don't Allocate the Global Number attribute for the elements and nodes
59   // It will be automatically created upon calls to void FEM_Entity::setGlobalno(int r,int g) {
60
61   FEM_Mesh_allocate_valid_attr(which_mesh, FEM_NODE);
62   FEM_Mesh_allocate_valid_attr(which_mesh, FEM_ELEM+0);
63
64   // Setup ghost layers, assume Tet4
65   // const int tet4vertices[4] = {0,1,2,3};
66   //  FEM_Add_ghost_layer(1,1);
67   // FEM_Add_ghost_elem(0,4,tet4vertices);
68
69   return model;
70 }
71
72 TopModel* topModel_Create_Driver(int elem_attr_sz, int node_attr_sz, int model_attr_sz, void *mAtt){
73     // This only uses a single mesh, so don't create multiple TopModels of these
74     CkAssert(elem_attr_sz > 0);
75     CkAssert(node_attr_sz > 0);
76     int which_mesh=FEM_Mesh_default_read();
77     TopModel *model = new TopModel;
78     model->elem_attr_size = elem_attr_sz;
79     model->node_attr_size = node_attr_sz;
80     model->model_attr_size = model_attr_sz;
81
82     model->mesh = FEM_Mesh_lookup(which_mesh,"TopModel::TopModel()");
83
84     model->mAtt = mAtt;
85
86 #if CUDA
87     /** Copy element Attribute array to device global memory */
88     {
89         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+0,"topModel_Create_Driver");
90         AllocTable2d<unsigned char> &dataTable  = at->getChar();
91         unsigned char *ElemData = dataTable.getData();
92         int size = at->size();
93         cudaMalloc(size, (void**)&(model->ElemDataDevice));
94         cudaMemcpy(ElemDataDevice,ElemData,size,cudaMemcpyHostToDevice);
95         model->num_local_elem = size / elem_attr_sz;
96     }
97
98     /** Copy node Attribute array to device global memory */
99     {
100         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+0,"topModel_Create_Driver");
101         AllocTable2d<unsigned char> &dataTable  = at->getChar();
102         unsigned char *NodeData = dataTable.getData();
103         int size = at->size();
104         cudaMalloc(size, (void**)&(model->NodeDataDevice));
105         cudaMemcpy(NodeDataDevice,NodeData,size,cudaMemcpyHostToDevice);
106         model->num_local_node = size / node_attr_sz;
107     }
108
109     /** Copy elem connectivity array to device global memory */
110     {
111         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_CONN,"topModel_Create_Driver");
112         AllocTable2d<int> &dataTable  = at->getInt();
113         int *data = dataTable.getData();
114         int size = at->size()*sizeof(int);
115         cudaMalloc(size, (void**)&(model->ElemConnDevice));
116         cudaMemcpy(ElemConnDevice,data,size,cudaMemcpyHostToDevice);
117     }
118
119     /** Copy model Attribute to device global memory */
120     {
121         cudaMalloc(model->model_attr_size, (void**)&(model->mAttDevice));
122         cudaMemcpy(mAttDevice,mAtt,model->model_attr_size,cudaMemcpyHostToDevice);
123     }
124
125     /** Copy model to device global memory */
126     {
127         cudaMalloc(sizeof(TopModel), (void**)&(model->modelD);
128         cudaMemcpy(model->modelD,model,sizeof(TopModel),cudaMemcpyHostToDevice);
129     }
130
131 #endif
132
133     return model;
134 }
135
136 /** Copy node attribute array from CUDA device back to the ParFUM attribute */
137 void top_retreive_node_data(TopModel* m){
138 #if CUDA
139     FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+0,"top_retreive_node_data");
140     AllocTable2d<unsigned char> &dataTable  = at->getChar();
141     unsigned char *NodeData = dataTable.getData();
142     int size = at->size();
143
144     cudaMemcpy(NodeData,m->NodeDataDevice,size,cudaMemcpyDeviceToHost);
145 #endif
146 }
147
148
149 /** Copy element attribute array from CUDA device back to the ParFUM attribute */
150 void top_retreive_elem_data(TopModel* m){
151 #if CUDA
152     FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+0,"top_retreive_elem_data");
153     AllocTable2d<unsigned char> &dataTable  = at->getChar();
154     unsigned char *ElemData = dataTable.getData();
155     int size = at->size();
156
157     cudaMemcpy(ElemData,m->ElemDataDevice,size,cudaMemcpyDeviceToHost);
158 #endif
159 }
160
161
162
163 /** Cleanup a model */
164 void topModel_Destroy(TopModel* m){
165 #if CUDA
166     cudaFree(m->mAttDevice);
167     cudaFree(m->NodeDataDevice);
168     cudaFree(m->ElemDataDevice);
169 #endif
170     delete m;
171 }
172
173
174 TopNode topModel_InsertNode(TopModel* m, FP_TYPE x, FP_TYPE y, FP_TYPE z){
175   int newNode = FEM_add_node_local(m->mesh,false,false,false);
176   m->mesh->node.set_coord(newNode,x,y,z);
177   return newNode;
178 }
179
180
181 /** Set id of a node
182 @todo Make this work with ghosts
183 */
184 void topNode_SetId(TopModel* m, TopNode n, TopID id){
185   CkAssert(n>=0);
186   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+1,"topNode_SetId");
187   at->getInt()(n,0)=id;
188 }
189
190 /** Insert an element */
191 TopElement topModel_InsertElem(TopModel*m, TopElemType type, TopNode* nodes){
192   CkAssert(type ==  TOP_ELEMENT_TET4);
193   int conn[4];
194   conn[0] = nodes[0];
195   conn[1] = nodes[1];
196   conn[2] = nodes[2];
197   conn[3] = nodes[3];
198   int newEl = FEM_add_element_local(m->mesh, conn, 4, 0, 0, 0);
199   return newEl;
200 }
201
202 /** Set id of an element
203 @todo Make this work with ghosts
204 */
205 void topElement_SetId(TopModel* m, TopElement e, TopID id){
206   CkAssert(e>=0);
207   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+1,"topElement_SetId");
208   at->getInt()(e,0)=id;
209 }
210
211
212
213 /**
214         @brief Set attribute of a node
215
216         The attribute passed in must be a contiguous data structure with size equal to the value node_attr_sz passed into topModel_Create_Driver() and topModel_Create_Init()
217
218         The supplied attribute will be copied into the ParFUM attribute array "FEM_DATA+0". Then ParFUM will own this data. The function topNode_GetAttrib() will return a pointer to the copy owned by ParFUM. If a single material parameter attribute is used for multiple nodes, each node will get a separate copy of the array. Any subsequent modifications to the data will only be reflected at a single node.
219
220         The user is responsible for deallocating parameter d passed into this function.
221
222 */
223 void topNode_SetAttrib(TopModel* m, TopNode n, void* d){
224   CkAssert(n>=0);
225   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+0,"topNode_SetAttrib");
226   AllocTable2d<unsigned char> &dataTable  = at->getChar();
227   unsigned char *data = dataTable.getData();
228   memcpy(data + n*m->node_attr_size, d, m->node_attr_size);
229 }
230
231 /** @brief Set attribute of an element
232 See topNode_SetAttrib() for description
233 */
234 void topElement_SetAttrib(TopModel* m, TopElement e, void* d){
235   CkAssert(e>=0);
236   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+0,"topElem_SetAttrib");
237   AllocTable2d<unsigned char> &dataTable  = at->getChar();
238   unsigned char *data = dataTable.getData();
239   memcpy(data + e*m->elem_attr_size, d, m->elem_attr_size);
240 }
241
242
243 /** @brief Get elem attribute
244 See topNode_SetAttrib() for description
245 */
246 void* topElement_GetAttrib(TopModel* m, TopElement e){
247   if(! m->mesh->elem[0].is_valid_any_idx(e))
248         return NULL;
249
250   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+0,"topElem_GetAttrib");
251   AllocTable2d<unsigned char> &dataTable  = at->getChar();
252   unsigned char *data = dataTable.getData();
253   return (data + e*m->elem_attr_size);
254 }
255
256 /** @brief Get nodal attribute
257 See topNode_SetAttrib() for description
258 */
259 void* topNode_GetAttrib(TopModel* m, TopNode n){
260   if(! m->mesh->node.is_valid_any_idx(n))
261         return NULL;
262
263   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+0,"topNode_GetAttrib");
264   AllocTable2d<unsigned char> &dataTable  = at->getChar();
265   unsigned char *data = dataTable.getData();
266   return (data + n*m->node_attr_size);
267 }
268
269
270
271 /**
272         Get node via id
273         @todo Re-implement this function with some hash to make it fast.
274         @note In the explicit FEA example, this is just used during initialization, so speed is not too important.
275         @todo Does not work with ghosts yet.
276 */
277 TopNode topModel_GetNodeAtId(TopModel* m, TopID id){
278   // lookup node via global ID
279   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+1,"topModel_GetNodeAtId");
280   for(int i=0;i<at->getInt().size();++i){
281         if(at->getInt()(i,0)==id){
282           return i;
283         }
284   }
285   return -1;
286 }
287
288 /**
289         Get elem via id
290         @todo Implement this function
291         @note Is this function even supposed to exist?
292  */
293 TopElement topModel_GetElemAtId(TopModel*m,TopID id){
294   CkPrintf("Error: topModel_GetElemAtId() not yet implemented");
295   CkExit();
296 }
297
298
299
300 TopNode topElement_GetNode(TopModel* m,TopElement e,int idx){
301   CkAssert(e>=0);
302   const AllocTable2d<int> &conn = m->mesh->elem[0].getConn();
303   CkAssert(idx>=0 && idx<conn.width());
304   CkAssert(idx<conn.size());
305
306   int node = conn(e,idx);
307
308   return conn(e,idx);
309 }
310
311 int topNode_GetId(TopModel* m, TopNode n){
312   CkAssert(n>=0);
313   FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+1,"topNode_SetId");
314   return at->getInt()(n,0);
315 }
316
317
318 /** @todo handle ghost nodes as appropriate */
319 int topModel_GetNNodes(TopModel *model){
320   return model->mesh->node.count_valid();
321 }
322
323 /** @todo Fix to return the width of the conn array */
324 int topElement_GetNNodes(TopModel* model, TopElement elem){
325   return 4;
326 }
327
328 /** @todo make sure we are in a getting mesh */
329 void topNode_GetPosition(TopModel*model, TopNode node,FP_TYPE*x,FP_TYPE*y,FP_TYPE*z){
330   CkAssert(node>=0);
331
332   // lookup node via global ID
333   FEM_DataAttribute * at = (FEM_DataAttribute*) model->mesh->node.lookup(FEM_COORD,"topModel_GetNodeAtId");
334
335 #ifdef FP_TYPE_FLOAT
336   CkAssert(at->getFloat().width()==3);
337   CkAssert(node<at->getFloat().size());
338   *x = at->getFloat()(node,0);
339   *y = at->getFloat()(node,1);
340   *z = at->getFloat()(node,2);
341 #else
342   CkAssert(at->getDouble().width()==3);
343   CkAssert(node<at->getDouble().size());
344   *x = at->getDouble()(node,0);
345   *y = at->getDouble()(node,1);
346   *z = at->getDouble()(node,2);
347 #endif
348
349 }
350
351 void topModel_Sync(TopModel*m){
352   MPI_Barrier(MPI_COMM_WORLD);
353
354
355   //  CkPrintf("%d: %d local, %d ghost elements\n", FEM_My_partition(), m->mesh->elem[0].size(),m->mesh->elem[0].ghost->size() );
356   //  CkPrintf("%d: %d local, %d ghost valid elements\n", FEM_My_partition(), m->mesh->elem[0].count_valid(),m->mesh->elem[0].ghost->count_valid() );
357
358 }
359
360 void topModel_TestIterators(TopModel*m){
361   CkAssert(m->mesh->elem[0].ghost!=NULL);
362   CkAssert(m->mesh->node.ghost!=NULL);
363
364   int expected_elem_count = m->mesh->elem[0].count_valid() + m->mesh->elem[0].ghost->count_valid();
365   int iterated_elem_count = 0;
366
367   int expected_node_count = m->mesh->node.count_valid() + m->mesh->node.ghost->count_valid();
368   int iterated_node_count = 0;
369
370   int myId = FEM_My_partition();
371
372
373   TopNodeItr* itr = topModel_CreateNodeItr(m);
374   for(topNodeItr_Begin(itr);topNodeItr_IsValid(itr);topNodeItr_Next(itr)){
375         iterated_node_count++;
376         TopNode node = topNodeItr_GetCurr(itr);
377         void* na = topNode_GetAttrib(m,node);
378         CkAssert(na != NULL);
379   }
380
381   TopElemItr* e_itr = topModel_CreateElemItr(m);
382   for(topElemItr_Begin(e_itr);topElemItr_IsValid(e_itr);topElemItr_Next(e_itr)){
383         iterated_elem_count++;
384         TopElement elem = topElemItr_GetCurr(e_itr);
385         void* ea = topElement_GetAttrib(m,elem);
386         CkAssert(ea != NULL);
387   }
388
389   CkAssert(iterated_node_count == expected_node_count);
390   CkAssert(iterated_elem_count==expected_elem_count);
391
392 }
393
394
395 #include "ParFUM_TOPS.def.h"