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