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