*** empty log message ***
[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, double x, double y, double 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 TopNode topModel_InsertNode(TopModel* m, float x, float y, float z){
192   int newNode = FEM_add_node_local(m->mesh,false,false,false);
193   (*m->coord_T)(newNode,0)=x;
194   (*m->coord_T)(newNode,1)=y;
195   (*m->coord_T)(newNode,2)=z;
196   return newNode;
197 }
198
199
200 /** Set id of a node
201 @todo Make this work with ghosts
202 */
203 void topNode_SetId(TopModel* m, TopNode n, TopID id){
204   CkAssert(n>=0);
205   (*m->node_id_T)(n,0)=id;
206 }
207
208 /** Insert an element */
209 TopElement topModel_InsertElem(TopModel*m, TopElemType type, TopNode* nodes){
210   CkAssert(type ==  TOP_ELEMENT_TET4);
211   int conn[4];
212   conn[0] = nodes[0];
213   conn[1] = nodes[1];
214   conn[2] = nodes[2];
215   conn[3] = nodes[3];
216   int newEl = FEM_add_element_local(m->mesh, conn, 4, 0, 0, 0);
217   return newEl;
218 }
219
220 /** Set id of an element
221 @todo Make this work with ghosts
222 */
223 void topElement_SetId(TopModel* m, TopElement e, TopID id){
224   CkAssert(e>=0);
225   (*m->elem_id_T)(e,0)=id;
226 }
227
228
229
230 /**
231         @brief Set attribute of a node
232
233         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()
234
235         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.
236
237         The user is responsible for deallocating parameter d passed into this function.
238
239 */
240 void topNode_SetAttrib(TopModel* m, TopNode n, void* d){
241   CkAssert(n>=0);
242   unsigned char *data = m->NodeData_T->getData();
243   memcpy(data + n*m->node_attr_size, d, m->node_attr_size);
244 }
245
246 /** @brief Set attribute of an element
247 See topNode_SetAttrib() for description
248 */
249 void topElement_SetAttrib(TopModel* m, TopElement e, void* d){
250   CkAssert(e>=0);
251   unsigned char *data = m->ElemData_T->getData();
252   memcpy(data + e*m->elem_attr_size, d, m->elem_attr_size);
253 }
254
255
256 /** @brief Get elem attribute
257 See topNode_SetAttrib() for description
258 */
259 void* topElement_GetAttrib(TopModel* m, TopElement e){
260   if(! m->mesh->elem[0].is_valid_any_idx(e))
261         return NULL;
262   unsigned char *data = m->ElemData_T->getData();
263   return (data + e*m->elem_attr_size);
264 }
265
266 /** @brief Get nodal attribute
267 See topNode_SetAttrib() for description
268 */
269 void* topNode_GetAttrib(TopModel* m, TopNode n){
270   if(! m->mesh->node.is_valid_any_idx(n))
271         return NULL;
272   unsigned char *data = m->NodeData_T->getData();
273   return (data + n*m->node_attr_size);
274 }
275
276
277
278 /**
279         Get node via id
280         @todo Re-implement this function with some hash to make it fast.
281         @note In the explicit FEA example, this is just used during initialization, so speed is not too important.
282         @todo Does not work with ghosts yet.
283 */
284 TopNode topModel_GetNodeAtId(TopModel* m, TopID id){
285   // lookup node via global ID
286   for(int i=0;i<m->node_id_T->size();++i){
287         if((*m->node_id_T)(i,0)==id){
288           return i;
289         }
290   }
291   return -1;
292 }
293
294 /**
295         Get elem via id
296         @todo Implement this function
297         @note Is this function even supposed to exist?
298  */
299 TopElement topModel_GetElemAtId(TopModel*m,TopID id){
300   CkPrintf("Error: topModel_GetElemAtId() not yet implemented");
301   CkExit();
302 }
303
304
305
306 TopNode topElement_GetNode(TopModel* m,TopElement e,int idx){
307   CkAssert(e>=0);
308   const AllocTable2d<int> &conn = m->mesh->elem[0].getConn();
309   CkAssert(idx>=0 && idx<conn.width());
310   CkAssert(idx<conn.size());
311
312   int node = conn(e,idx);
313
314   return conn(e,idx);
315 }
316
317 int topNode_GetId(TopModel* m, TopNode n){
318   CkAssert(n>=0);
319   return (*m->node_id_T)(n,0);
320 }
321
322
323 /** @todo handle ghost nodes as appropriate */
324 int topModel_GetNNodes(TopModel *model){
325   return model->mesh->node.count_valid();
326 }
327
328 /** @todo Fix to return the width of the conn array */
329 int topElement_GetNNodes(TopModel* model, TopElement elem){
330   return 4;
331 }
332
333 /** @todo make sure we are in a getting mesh */
334 void topNode_GetPosition(TopModel*model, TopNode node,double*x,double*y,double*z){
335   CkAssert(node>=0);
336   *x = (*model->coord_T)(node,0);
337   *y = (*model->coord_T)(node,1);
338   *z = (*model->coord_T)(node,2);
339 }
340
341 /** @todo make sure we are in a getting mesh */
342 void topNode_GetPosition(TopModel*model, TopNode node,float*x,float*y,float*z){
343   CkAssert(node>=0);
344   *x = (*model->coord_T)(node,0);
345   *y = (*model->coord_T)(node,1);
346   *z = (*model->coord_T)(node,2);
347 }
348
349 void topModel_Sync(TopModel*m){
350   MPI_Barrier(MPI_COMM_WORLD);
351
352
353   //  CkPrintf("%d: %d local, %d ghost elements\n", FEM_My_partition(), m->mesh->elem[0].size(),m->mesh->elem[0].ghost->size() );
354   //  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() );
355
356 }
357
358 void topModel_TestIterators(TopModel*m){
359   CkAssert(m->mesh->elem[0].ghost!=NULL);
360   CkAssert(m->mesh->node.ghost!=NULL);
361
362   int expected_elem_count = m->mesh->elem[0].count_valid() + m->mesh->elem[0].ghost->count_valid();
363   int iterated_elem_count = 0;
364
365   int expected_node_count = m->mesh->node.count_valid() + m->mesh->node.ghost->count_valid();
366   int iterated_node_count = 0;
367
368   int myId = FEM_My_partition();
369
370
371   TopNodeItr* itr = topModel_CreateNodeItr(m);
372   for(topNodeItr_Begin(itr);topNodeItr_IsValid(itr);topNodeItr_Next(itr)){
373         iterated_node_count++;
374         TopNode node = topNodeItr_GetCurr(itr);
375         void* na = topNode_GetAttrib(m,node);
376         CkAssert(na != NULL);
377   }
378
379   TopElemItr* e_itr = topModel_CreateElemItr(m);
380   for(topElemItr_Begin(e_itr);topElemItr_IsValid(e_itr);topElemItr_Next(e_itr)){
381         iterated_elem_count++;
382         TopElement elem = topElemItr_GetCurr(e_itr);
383         void* ea = topElement_GetAttrib(m,elem);
384         CkAssert(ea != NULL);
385   }
386
387   CkAssert(iterated_node_count == expected_node_count);
388   CkAssert(iterated_elem_count==expected_elem_count);
389
390 }
391
392
393 #include "ParFUM_TOPS.def.h"