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