A working version of the layer.
[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 element Attribute array to device global memory */
110     {
111         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->elem[0].lookup(FEM_DATA+0,"topModel_Create_Driver");
112         AllocTable2d<unsigned char> &dataTable  = at->getChar();
113         unsigned char *ElemData = dataTable.getData();
114         int size = at->size();
115         cudaMalloc(size, (void**)&(model->ElemDataDevice));
116         cudaMemcpy(ElemDataDevice,ElemData,size,cudaMemcpyHostToDevice);
117     }
118
119     /** Copy node Attribute array to device global memory */
120     {
121         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_DATA+0,"topModel_Create_Driver");
122         AllocTable2d<unsigned char> &dataTable  = at->getChar();
123         unsigned char *NodeData = dataTable.getData();
124         int size = at->size();
125         cudaMalloc(size, (void**)&(model->NodeDataDevice));
126         cudaMemcpy(NodeDataDevice,NodeData,size,cudaMemcpyHostToDevice);
127     }
128
129     /** Copy elem connectivity array to device global memory */
130     {
131         FEM_DataAttribute * at = (FEM_DataAttribute*) m->mesh->node.lookup(FEM_CONN,"topModel_Create_Driver");
132         AllocTable2d<int> &dataTable  = at->getInt();
133         int *data = dataTable.getData();
134         int size = at->size()*sizeof(int);
135         cudaMalloc(size, (void**)&(model->ElemConnDevice));
136         cudaMemcpy(ElemConnDevice,data,size,cudaMemcpyHostToDevice);
137     }
138
139     /** Copy model Attribute to device global memory */
140     {
141         cudaMalloc(model->model_attr_size, (void**)&(model->mAttDevice));
142         cudaMemcpy(mAttDevice,mAtt,model->model_attr_size,cudaMemcpyHostToDevice);
143     }
144
145     /** Copy model to device global memory */
146     {
147         cudaMalloc(sizeof(TopModel), (void**)&(model->modelD);
148         cudaMemcpy(model->modelD,model,sizeof(TopModel),cudaMemcpyHostToDevice);
149     }
150 #endif
151
152     return model;
153 }
154
155 /** Copy node attribute array from CUDA device back to the ParFUM attribute */
156 void top_retreive_node_data(TopModel* m){
157 #if CUDA
158     cudaMemcpy(m->NodeData,m->NodeDataDevice,size,cudaMemcpyDeviceToHost);
159 #endif
160 }
161
162
163 /** Copy element attribute array from CUDA device back to the ParFUM attribute */
164 void top_retreive_elem_data(TopModel* m){
165 #if CUDA
166     cudaMemcpy(m->ElemData,m->ElemDataDevice,size,cudaMemcpyDeviceToHost);
167 #endif
168 }
169
170
171
172 /** Cleanup a model */
173 void topModel_Destroy(TopModel* m){
174 #if CUDA
175     cudaFree(m->mAttDevice);
176     cudaFree(m->NodeDataDevice);
177     cudaFree(m->ElemDataDevice);
178 #endif
179     delete m;
180 }
181
182
183 TopNode topModel_InsertNode(TopModel* m, FP_TYPE x, FP_TYPE y, FP_TYPE z){
184   int newNode = FEM_add_node_local(m->mesh,false,false,false);
185   (*m->coord_T)(newNode,0)=x;
186   (*m->coord_T)(newNode,1)=y;
187   (*m->coord_T)(newNode,2)=z;
188   return newNode;
189 }
190
191
192 /** Set id of a node
193 @todo Make this work with ghosts
194 */
195 void topNode_SetId(TopModel* m, TopNode n, TopID id){
196   CkAssert(n>=0);
197   (*m->node_id_T)(n,0)=id;
198 }
199
200 /** Insert an element */
201 TopElement topModel_InsertElem(TopModel*m, TopElemType type, TopNode* nodes){
202   CkAssert(type ==  TOP_ELEMENT_TET4);
203   int conn[4];
204   conn[0] = nodes[0];
205   conn[1] = nodes[1];
206   conn[2] = nodes[2];
207   conn[3] = nodes[3];
208   int newEl = FEM_add_element_local(m->mesh, conn, 4, 0, 0, 0);
209   return newEl;
210 }
211
212 /** Set id of an element
213 @todo Make this work with ghosts
214 */
215 void topElement_SetId(TopModel* m, TopElement e, TopID id){
216   CkAssert(e>=0);
217   (*m->elem_id_T)(e,0)=id;
218 }
219
220
221
222 /**
223         @brief Set attribute of a node
224
225         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()
226
227         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.
228
229         The user is responsible for deallocating parameter d passed into this function.
230
231 */
232 void topNode_SetAttrib(TopModel* m, TopNode n, void* d){
233   CkAssert(n>=0);
234   unsigned char *data = m->NodeData_T->getData();
235   memcpy(data + n*m->node_attr_size, d, m->node_attr_size);
236 }
237
238 /** @brief Set attribute of an element
239 See topNode_SetAttrib() for description
240 */
241 void topElement_SetAttrib(TopModel* m, TopElement e, void* d){
242   CkAssert(e>=0);
243   unsigned char *data = m->ElemData_T->getData();
244   memcpy(data + e*m->elem_attr_size, d, m->elem_attr_size);
245 }
246
247
248 /** @brief Get elem attribute
249 See topNode_SetAttrib() for description
250 */
251 void* topElement_GetAttrib(TopModel* m, TopElement e){
252   if(! m->mesh->elem[0].is_valid_any_idx(e))
253         return NULL;
254   unsigned char *data = m->ElemData_T->getData();
255   return (data + e*m->elem_attr_size);
256 }
257
258 /** @brief Get nodal attribute
259 See topNode_SetAttrib() for description
260 */
261 void* topNode_GetAttrib(TopModel* m, TopNode n){
262   if(! m->mesh->node.is_valid_any_idx(n))
263         return NULL;
264   unsigned char *data = m->NodeData_T->getData();
265   return (data + n*m->node_attr_size);
266 }
267
268
269
270 /**
271         Get node via id
272         @todo Re-implement this function with some hash to make it fast.
273         @note In the explicit FEA example, this is just used during initialization, so speed is not too important.
274         @todo Does not work with ghosts yet.
275 */
276 TopNode topModel_GetNodeAtId(TopModel* m, TopID id){
277   // lookup node via global ID
278   for(int i=0;i<m->node_id_T->size();++i){
279         if((*m->node_id_T)(i,0)==id){
280           return i;
281         }
282   }
283   return -1;
284 }
285
286 /**
287         Get elem via id
288         @todo Implement this function
289         @note Is this function even supposed to exist?
290  */
291 TopElement topModel_GetElemAtId(TopModel*m,TopID id){
292   CkPrintf("Error: topModel_GetElemAtId() not yet implemented");
293   CkExit();
294 }
295
296
297
298 TopNode topElement_GetNode(TopModel* m,TopElement e,int idx){
299   CkAssert(e>=0);
300   const AllocTable2d<int> &conn = m->mesh->elem[0].getConn();
301   CkAssert(idx>=0 && idx<conn.width());
302   CkAssert(idx<conn.size());
303
304   int node = conn(e,idx);
305
306   return conn(e,idx);
307 }
308
309 int topNode_GetId(TopModel* m, TopNode n){
310   CkAssert(n>=0);
311   return (*m->node_id_T)(n,0);
312 }
313
314
315 /** @todo handle ghost nodes as appropriate */
316 int topModel_GetNNodes(TopModel *model){
317   return model->mesh->node.count_valid();
318 }
319
320 /** @todo Fix to return the width of the conn array */
321 int topElement_GetNNodes(TopModel* model, TopElement elem){
322   return 4;
323 }
324
325 /** @todo make sure we are in a getting mesh */
326 void topNode_GetPosition(TopModel*model, TopNode node,FP_TYPE*x,FP_TYPE*y,FP_TYPE*z){
327   CkAssert(node>=0);
328   *x = (*model->coord_T)(node,0);
329   *y = (*model->coord_T)(node,1);
330   *z = (*model->coord_T)(node,2);
331 }
332
333 void topModel_Sync(TopModel*m){
334   MPI_Barrier(MPI_COMM_WORLD);
335
336
337   //  CkPrintf("%d: %d local, %d ghost elements\n", FEM_My_partition(), m->mesh->elem[0].size(),m->mesh->elem[0].ghost->size() );
338   //  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() );
339
340 }
341
342 void topModel_TestIterators(TopModel*m){
343   CkAssert(m->mesh->elem[0].ghost!=NULL);
344   CkAssert(m->mesh->node.ghost!=NULL);
345
346   int expected_elem_count = m->mesh->elem[0].count_valid() + m->mesh->elem[0].ghost->count_valid();
347   int iterated_elem_count = 0;
348
349   int expected_node_count = m->mesh->node.count_valid() + m->mesh->node.ghost->count_valid();
350   int iterated_node_count = 0;
351
352   int myId = FEM_My_partition();
353
354
355   TopNodeItr* itr = topModel_CreateNodeItr(m);
356   for(topNodeItr_Begin(itr);topNodeItr_IsValid(itr);topNodeItr_Next(itr)){
357         iterated_node_count++;
358         TopNode node = topNodeItr_GetCurr(itr);
359         void* na = topNode_GetAttrib(m,node);
360         CkAssert(na != NULL);
361   }
362
363   TopElemItr* e_itr = topModel_CreateElemItr(m);
364   for(topElemItr_Begin(e_itr);topElemItr_IsValid(e_itr);topElemItr_Next(e_itr)){
365         iterated_elem_count++;
366         TopElement elem = topElemItr_GetCurr(e_itr);
367         void* ea = topElement_GetAttrib(m,elem);
368         CkAssert(ea != NULL);
369   }
370
371   CkAssert(iterated_node_count == expected_node_count);
372   CkAssert(iterated_elem_count==expected_elem_count);
373
374 }
375
376
377 #include "ParFUM_TOPS.def.h"