/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1339 by ksteube, Wed Nov 7 01:53:12 2007 UTC revision 1455 by phornby, Thu Feb 28 17:19:44 2008 UTC
# Line 16  Line 16 
16  #include "MeshAdapter.h"  #include "MeshAdapter.h"
17  #include "escript/Data.h"  #include "escript/Data.h"
18  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
19    #ifdef USE_NETCDF
20    #include <netcdfcpp.h>
21    #endif
22  extern "C" {  extern "C" {
23  #include "escript/blocktimer.h"  #include "escript/blocktimer.h"
24  }  }
25  #include <vector>  #include <vector>
26    
27    #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256
28    
29  using namespace std;  using namespace std;
30  using namespace escript;  using namespace escript;
31    
# Line 100  void MeshAdapter::Print_Mesh_Info(const Line 105  void MeshAdapter::Print_Mesh_Info(const
105    
106  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const std::string& fileName) const
107  {  {
108    char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;  #ifdef USE_NETCDF
109    strcpy(fName,fileName.c_str());     const NcDim* ncdims[12];
110    Finley_Mesh_dump(m_finleyMesh.get(),fName);     NcVar *ids, *data;
111    checkFinleyError();     int *int_ptr;
112    TMPMEMFREE(fName);     Finley_Mesh *mesh = m_finleyMesh.get();
113       Finley_TagMap* tag_map;
114       int num_Tags = 0;
115       int mpi_size             = mesh->MPIInfo->size;
116       int mpi_rank             = mesh->MPIInfo->rank;
117       int numDim               = mesh->Nodes->numDim;
118       int numNodes             = mesh->Nodes->numNodes;
119       int num_Elements         = mesh->Elements->numElements;
120       int num_FaceElements         = mesh->FaceElements->numElements;
121       int num_ContactElements      = mesh->ContactElements->numElements;
122       int num_Points           = mesh->Points->numElements;
123       int num_Elements_numNodes        = mesh->Elements->numNodes;
124       int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
125       int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
126       char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);
127    
128       /* Figure out how much storage is required for tags */
129       tag_map = mesh->TagMap;
130       if (tag_map) {
131         while (tag_map) {
132        num_Tags++;
133            tag_map=tag_map->next;
134         }
135       }
136    
137       // NetCDF error handler
138       NcError err(NcError::verbose_nonfatal);
139       // Create the file.
140       NcFile dataFile(newFileName, NcFile::Replace);
141       // check if writing was successful
142       if (!dataFile.is_valid())
143            throw DataException("Error - MeshAdapter::dump: opening of NetCDF file for output failed: " + *newFileName);
144    
145       // Define dimensions (num_Elements and dim_Elements are identical, dim_Elements only appears if > 0)
146       if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
147            throw DataException("Error - MeshAdapter::dump: appending dimension numNodes to netCDF file failed: " + *newFileName);
148       if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
149            throw DataException("Error - MeshAdapter::dump: appending dimension numDim to netCDF file failed: " + *newFileName);
150       if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
151            throw DataException("Error - MeshAdapter::dump: appending dimension mpi_size to netCDF file failed: " + *newFileName);
152       if (num_Elements>0)
153          if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
154            throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements to netCDF file failed: " + *newFileName);
155       if (num_FaceElements>0)
156          if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
157            throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements to netCDF file failed: " + *newFileName);
158       if (num_ContactElements>0)
159          if (! (ncdims[5] = dataFile.add_dim("dim_ContactElements", num_ContactElements)) )
160            throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements to netCDF file failed: " + *newFileName);
161       if (num_Points>0)
162          if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
163            throw DataException("Error - MeshAdapter::dump: appending dimension dim_Points to netCDF file failed: " + *newFileName);
164       if (num_Elements>0)
165          if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
166            throw DataException("Error - MeshAdapter::dump: appending dimension dim_Elements_Nodes to netCDF file failed: " + *newFileName);
167       if (num_FaceElements>0)
168          if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
169            throw DataException("Error - MeshAdapter::dump: appending dimension dim_FaceElements_numNodes to netCDF file failed: " + *newFileName);
170       if (num_ContactElements>0)
171          if (! (ncdims[9] = dataFile.add_dim("dim_ContactElements_numNodes", num_ContactElements_numNodes)) )
172            throw DataException("Error - MeshAdapter::dump: appending dimension dim_ContactElements_numNodes to netCDF file failed: " + *newFileName);
173       if (num_Tags>0)
174          if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
175            throw DataException("Error - MeshAdapter::dump: appending dimension dim_Tags to netCDF file failed: " + *newFileName);
176    
177       // Attributes: MPI size, MPI rank, Name, order, reduced_order
178       if (!dataFile.add_att("mpi_size", mpi_size) )
179            throw DataException("Error - MeshAdapter::dump: appending mpi_size to NetCDF file failed: " + *newFileName);
180       if (!dataFile.add_att("mpi_rank", mpi_rank) )
181            throw DataException("Error - MeshAdapter::dump: appending mpi_rank to NetCDF file failed: " + *newFileName);
182       if (!dataFile.add_att("Name",mesh->Name) )
183            throw DataException("Error - MeshAdapter::dump: appending Name to NetCDF file failed: " + *newFileName);
184       if (!dataFile.add_att("numDim",numDim) )
185            throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
186       if (!dataFile.add_att("order",mesh->order) )
187            throw DataException("Error - MeshAdapter::dump: appending order to NetCDF file failed: " + *newFileName);
188       if (!dataFile.add_att("reduced_order",mesh->reduced_order) )
189            throw DataException("Error - MeshAdapter::dump: appending reduced_order to NetCDF file failed: " + *newFileName);
190       if (!dataFile.add_att("numNodes",numNodes) )
191            throw DataException("Error - MeshAdapter::dump: appending numNodes to NetCDF file failed: " + *newFileName);
192       if (!dataFile.add_att("num_Elements",num_Elements) )
193            throw DataException("Error - MeshAdapter::dump: appending num_Elements to NetCDF file failed: " + *newFileName);
194       if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
195            throw DataException("Error - MeshAdapter::dump: appending num_FaceElements to NetCDF file failed: " + *newFileName);
196       if (!dataFile.add_att("num_ContactElements",num_ContactElements) )
197            throw DataException("Error - MeshAdapter::dump: appending num_ContactElements to NetCDF file failed: " + *newFileName);
198       if (!dataFile.add_att("num_Points",num_Points) )
199            throw DataException("Error - MeshAdapter::dump: appending num_Points to NetCDF file failed: " + *newFileName);
200       if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
201            throw DataException("Error - MeshAdapter::dump: appending num_Elements_numNodes to NetCDF file failed: " + *newFileName);
202       if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
203            throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);
204       if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
205            throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);
206       if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )
207          throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
208       if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )
209          throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
210       if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )
211          throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
212       if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )
213          throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);
214       if (!dataFile.add_att("num_Tags", num_Tags) )
215          throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);
216    
217       // // // // // Nodes // // // // //
218    
219       // Only write nodes if non-empty because NetCDF doesn't like empty arrays (it treats them as NC_UNLIMITED)
220       if (numNodes>0) {
221    
222       // Nodes Id
223       if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
224            throw DataException("Error - MeshAdapter::dump: appending Nodes_Id to netCDF file failed: " + *newFileName);
225       int_ptr = &mesh->Nodes->Id[0];
226       if (! (ids->put(int_ptr, numNodes)) )
227            throw DataException("Error - MeshAdapter::dump: copy Nodes_Id to netCDF buffer failed: " + *newFileName);
228    
229       // Nodes Tag
230       if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
231            throw DataException("Error - MeshAdapter::dump: appending Nodes_Tag to netCDF file failed: " + *newFileName);
232       int_ptr = &mesh->Nodes->Tag[0];
233       if (! (ids->put(int_ptr, numNodes)) )
234            throw DataException("Error - MeshAdapter::dump: copy Nodes_Tag to netCDF buffer failed: " + *newFileName);
235    
236       // Nodes gDOF
237       if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
238            throw DataException("Error - MeshAdapter::dump: appending Nodes_gDOF to netCDF file failed: " + *newFileName);
239       int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
240       if (! (ids->put(int_ptr, numNodes)) )
241            throw DataException("Error - MeshAdapter::dump: copy Nodes_gDOF to netCDF buffer failed: " + *newFileName);
242    
243       // Nodes global node index
244       if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
245            throw DataException("Error - MeshAdapter::dump: appending Nodes_gNI to netCDF file failed: " + *newFileName);
246       int_ptr = &mesh->Nodes->globalNodesIndex[0];
247       if (! (ids->put(int_ptr, numNodes)) )
248            throw DataException("Error - MeshAdapter::dump: copy Nodes_gNI to netCDF buffer failed: " + *newFileName);
249    
250       // Nodes grDof
251       if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
252            throw DataException("Error - MeshAdapter::dump: appending Nodes_grDfI to netCDF file failed: " + *newFileName);
253       int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
254       if (! (ids->put(int_ptr, numNodes)) )
255            throw DataException("Error - MeshAdapter::dump: copy Nodes_grDfI to netCDF buffer failed: " + *newFileName);
256    
257       // Nodes grNI
258       if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
259            throw DataException("Error - MeshAdapter::dump: appending Nodes_grNI to netCDF file failed: " + *newFileName);
260       int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
261       if (! (ids->put(int_ptr, numNodes)) )
262            throw DataException("Error - MeshAdapter::dump: copy Nodes_grNI to netCDF buffer failed: " + *newFileName);
263    
264       // Nodes Coordinates
265       if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
266            throw DataException("Error - MeshAdapter::dump: appending Nodes_Coordinates to netCDF file failed: " + *newFileName);
267       if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
268            throw DataException("Error - MeshAdapter::dump: copy Nodes_Coordinates to netCDF buffer failed: " + *newFileName);
269    
270       // Nodes degreesOfFreedomDistribution
271       if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
272            throw DataException("Error - MeshAdapter::dump: appending Nodes_DofDistribution to netCDF file failed: " + *newFileName);
273       int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
274       if (! (ids->put(int_ptr, mpi_size+1)) )
275            throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
276    
277       }
278    
279       // // // // // Elements // // // // //
280    
281       if (num_Elements>0) {
282    
283         // Temp storage to gather node IDs
284         int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);
285    
286         // Elements_Id
287         if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
288            throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);
289         int_ptr = &mesh->Elements->Id[0];
290         if (! (ids->put(int_ptr, num_Elements)) )
291            throw DataException("Error - MeshAdapter::dump: copy Elements_Id to netCDF buffer failed: " + *newFileName);
292    
293         // Elements_Tag
294         if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
295            throw DataException("Error - MeshAdapter::dump: appending Elements_Tag to netCDF file failed: " + *newFileName);
296         int_ptr = &mesh->Elements->Tag[0];
297         if (! (ids->put(int_ptr, num_Elements)) )
298            throw DataException("Error - MeshAdapter::dump: copy Elements_Tag to netCDF buffer failed: " + *newFileName);
299    
300         // Elements_Owner
301         if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
302            throw DataException("Error - MeshAdapter::dump: appending Elements_Owner to netCDF file failed: " + *newFileName);
303         int_ptr = &mesh->Elements->Owner[0];
304         if (! (ids->put(int_ptr, num_Elements)) )
305            throw DataException("Error - MeshAdapter::dump: copy Elements_Owner to netCDF buffer failed: " + *newFileName);
306    
307         // Elements_Color
308         if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
309            throw DataException("Error - MeshAdapter::dump: appending Elements_Color to netCDF file failed: " + *newFileName);
310         int_ptr = &mesh->Elements->Color[0];
311         if (! (ids->put(int_ptr, num_Elements)) )
312            throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
313    
314         // Elements_Nodes
315         for (int i=0; i<num_Elements; i++)
316           for (int j=0; j<num_Elements_numNodes; j++)
317             Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];
318         if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
319            throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
320         if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )
321            throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
322    
323         TMPMEMFREE(Elements_Nodes);
324    
325       }
326    
327       // // // // // Face_Elements // // // // //
328    
329       if (num_FaceElements>0) {
330    
331         // Temp storage to gather node IDs
332         int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);
333    
334         // FaceElements_Id
335         if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
336            throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);
337         int_ptr = &mesh->FaceElements->Id[0];
338         if (! (ids->put(int_ptr, num_FaceElements)) )
339            throw DataException("Error - MeshAdapter::dump: copy FaceElements_Id to netCDF buffer failed: " + *newFileName);
340    
341         // FaceElements_Tag
342         if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
343            throw DataException("Error - MeshAdapter::dump: appending FaceElements_Tag to netCDF file failed: " + *newFileName);
344         int_ptr = &mesh->FaceElements->Tag[0];
345         if (! (ids->put(int_ptr, num_FaceElements)) )
346            throw DataException("Error - MeshAdapter::dump: copy FaceElements_Tag to netCDF buffer failed: " + *newFileName);
347    
348         // FaceElements_Owner
349         if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
350            throw DataException("Error - MeshAdapter::dump: appending FaceElements_Owner to netCDF file failed: " + *newFileName);
351         int_ptr = &mesh->FaceElements->Owner[0];
352         if (! (ids->put(int_ptr, num_FaceElements)) )
353            throw DataException("Error - MeshAdapter::dump: copy FaceElements_Owner to netCDF buffer failed: " + *newFileName);
354    
355         // FaceElements_Color
356         if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
357            throw DataException("Error - MeshAdapter::dump: appending FaceElements_Color to netCDF file failed: " + *newFileName);
358         int_ptr = &mesh->FaceElements->Color[0];
359         if (! (ids->put(int_ptr, num_FaceElements)) )
360            throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
361    
362         // FaceElements_Nodes
363         for (int i=0; i<num_FaceElements; i++)
364           for (int j=0; j<num_FaceElements_numNodes; j++)
365             FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];
366         if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
367            throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
368         if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
369            throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
370    
371         TMPMEMFREE(FaceElements_Nodes);
372    
373       }
374    
375       // // // // // Contact_Elements // // // // //
376    
377       if (num_ContactElements>0) {
378    
379         // Temp storage to gather node IDs
380         int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);
381    
382         // ContactElements_Id
383         if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
384            throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);
385         int_ptr = &mesh->ContactElements->Id[0];
386         if (! (ids->put(int_ptr, num_ContactElements)) )
387            throw DataException("Error - MeshAdapter::dump: copy ContactElements_Id to netCDF buffer failed: " + *newFileName);
388    
389         // ContactElements_Tag
390         if (! ( ids = dataFile.add_var("ContactElements_Tag", ncInt, ncdims[5])) )
391            throw DataException("Error - MeshAdapter::dump: appending ContactElements_Tag to netCDF file failed: " + *newFileName);
392         int_ptr = &mesh->ContactElements->Tag[0];
393         if (! (ids->put(int_ptr, num_ContactElements)) )
394            throw DataException("Error - MeshAdapter::dump: copy ContactElements_Tag to netCDF buffer failed: " + *newFileName);
395    
396         // ContactElements_Owner
397         if (! ( ids = dataFile.add_var("ContactElements_Owner", ncInt, ncdims[5])) )
398            throw DataException("Error - MeshAdapter::dump: appending ContactElements_Owner to netCDF file failed: " + *newFileName);
399         int_ptr = &mesh->ContactElements->Owner[0];
400         if (! (ids->put(int_ptr, num_ContactElements)) )
401            throw DataException("Error - MeshAdapter::dump: copy ContactElements_Owner to netCDF buffer failed: " + *newFileName);
402    
403         // ContactElements_Color
404         if (! ( ids = dataFile.add_var("ContactElements_Color", ncInt, ncdims[5])) )
405            throw DataException("Error - MeshAdapter::dump: appending ContactElements_Color to netCDF file failed: " + *newFileName);
406         int_ptr = &mesh->ContactElements->Color[0];
407         if (! (ids->put(int_ptr, num_ContactElements)) )
408            throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
409    
410         // ContactElements_Nodes
411         for (int i=0; i<num_ContactElements; i++)
412           for (int j=0; j<num_ContactElements_numNodes; j++)
413             ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];
414         if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
415            throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
416         if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
417            throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
418    
419         TMPMEMFREE(ContactElements_Nodes);
420    
421       }
422    
423       // // // // // Points // // // // //
424    
425       if (num_Points>0) {
426    
427         fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
428    
429         // Temp storage to gather node IDs
430         int *Points_Nodes = TMPMEMALLOC(num_Points,int);
431    
432         // Points_Id
433         if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
434            throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
435         int_ptr = &mesh->Points->Id[0];
436         if (! (ids->put(int_ptr, num_Points)) )
437            throw DataException("Error - MeshAdapter::dump: copy Points_Id to netCDF buffer failed: " + *newFileName);
438    
439         // Points_Tag
440         if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
441            throw DataException("Error - MeshAdapter::dump: appending Points_Tag to netCDF file failed: " + *newFileName);
442         int_ptr = &mesh->Points->Tag[0];
443         if (! (ids->put(int_ptr, num_Points)) )
444            throw DataException("Error - MeshAdapter::dump: copy Points_Tag to netCDF buffer failed: " + *newFileName);
445    
446         // Points_Owner
447         if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
448            throw DataException("Error - MeshAdapter::dump: appending Points_Owner to netCDF file failed: " + *newFileName);
449         int_ptr = &mesh->Points->Owner[0];
450         if (! (ids->put(int_ptr, num_Points)) )
451            throw DataException("Error - MeshAdapter::dump: copy Points_Owner to netCDF buffer failed: " + *newFileName);
452    
453         // Points_Color
454         if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
455            throw DataException("Error - MeshAdapter::dump: appending Points_Color to netCDF file failed: " + *newFileName);
456         int_ptr = &mesh->Points->Color[0];
457         if (! (ids->put(int_ptr, num_Points)) )
458            throw DataException("Error - MeshAdapter::dump: copy Points_Color to netCDF buffer failed: " + *newFileName);
459    
460         // Points_Nodes
461         // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
462         for (int i=0; i<num_Points; i++)
463           Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];
464         if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
465            throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
466         if (! (ids->put(&(Points_Nodes[0]), num_Points)) )
467            throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
468    
469         TMPMEMFREE(Points_Nodes);
470    
471       }
472    
473       // // // // // TagMap // // // // //
474    
475       if (num_Tags>0) {
476    
477         // Temp storage to gather node IDs
478         int *Tags_keys = TMPMEMALLOC(num_Tags, int);
479         char name_temp[4096];
480    
481         /* Copy tag data into temp arrays */
482         tag_map = mesh->TagMap;
483         if (tag_map) {
484           int i = 0;
485           while (tag_map) {
486        Tags_keys[i++] = tag_map->tag_key;
487            tag_map=tag_map->next;
488           }
489         }
490    
491         // Tags_keys
492         if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
493            throw DataException("Error - MeshAdapter::dump: appending Tags_keys to netCDF file failed: " + *newFileName);
494         int_ptr = &Tags_keys[0];
495         if (! (ids->put(int_ptr, num_Tags)) )
496            throw DataException("Error - MeshAdapter::dump: copy Tags_keys to netCDF buffer failed: " + *newFileName);
497    
498         // Tags_names_*
499         // This is an array of strings, it should be stored as an array but instead I have hacked in one attribute per string
500         // because the NetCDF manual doesn't tell how to do an array of strings
501         tag_map = mesh->TagMap;
502         if (tag_map) {
503           int i = 0;
504           while (tag_map) {
505             sprintf(name_temp, "Tags_name_%d", i);
506             if (!dataFile.add_att(name_temp, tag_map->name) )
507               throw DataException("Error - MeshAdapter::dump: appending Tags_names_ to NetCDF file failed: " + *newFileName);
508             tag_map=tag_map->next;
509         i++;
510           }
511         }
512    
513         TMPMEMFREE(Tags_keys);
514    
515       }
516    
517    
518       // NetCDF file is closed by destructor of NcFile object
519    #else
520       Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
521    #endif  /* USE_NETCDF */
522       checkFinleyError();
523  }  }
524    
525  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
# Line 396  void MeshAdapter::addPDEToRHS( escript:: Line 811  void MeshAdapter::addPDEToRHS( escript::
811     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
812     checkFinleyError();     checkFinleyError();
813  }  }
814    //
815    // adds PDE of second order into a transport problem
816    //
817    void MeshAdapter::addPDEToTransportProblem(
818                         TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
819                         const escript::Data& A, const escript::Data& B, const escript::Data& C,
820                         const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
821                         const escript::Data& d, const escript::Data& y,
822                         const escript::Data& d_contact,const escript::Data& y_contact) const
823    {
824       DataArrayView::ShapeType shape;
825       source.expand();
826       escriptDataC _source=source.getDataC();
827       escriptDataC _M=M.getDataC();
828       escriptDataC _A=A.getDataC();
829       escriptDataC _B=B.getDataC();
830       escriptDataC _C=C.getDataC();
831       escriptDataC _D=D.getDataC();
832       escriptDataC _X=X.getDataC();
833       escriptDataC _Y=Y.getDataC();
834       escriptDataC _d=d.getDataC();
835       escriptDataC _y=y.getDataC();
836       escriptDataC _d_contact=d_contact.getDataC();
837       escriptDataC _y_contact=y_contact.getDataC();
838    
839       Finley_Mesh* mesh=m_finleyMesh.get();
840       Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();
841      
842       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
843       checkFinleyError();
844    
845       Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
846       checkFinleyError();
847    
848       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
849       checkFinleyError();
850    
851       Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
852       checkFinleyError();
853    }
854    
855  //  //
856  // interpolates data between different function spaces:  // interpolates data between different function spaces:
# Line 504  void MeshAdapter::interpolateOnDomain(es Line 959  void MeshAdapter::interpolateOnDomain(es
959             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
960          } else {          } else {
961             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
962         }          }
963         break;          break;
964       case(ReducedFaceElements):       case(ReducedFaceElements):
965          if (target.getFunctionSpace().getTypeCode()==FaceElements) {          if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
966             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
967          } else {          } else {
968             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");             throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
# Line 528  void MeshAdapter::interpolateOnDomain(es Line 983  void MeshAdapter::interpolateOnDomain(es
983             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
984          } else {          } else {
985             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
            break;  
986          }          }
987          break;          break;
988       case(ReducedContactElementsZero):       case(ReducedContactElementsZero):
# Line 537  void MeshAdapter::interpolateOnDomain(es Line 991  void MeshAdapter::interpolateOnDomain(es
991             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
992          } else {          } else {
993             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");             throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
            break;  
994          }          }
995          break;          break;
996       case(DegreesOfFreedom):             case(DegreesOfFreedom):      
# Line 992  void MeshAdapter::setNewX(const escript: Line 1445  void MeshAdapter::setNewX(const escript:
1445  // saves a data array in openDX format:  // saves a data array in openDX format:
1446  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1447  {  {
1448      int MAX_namelength=256;      unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;
1449      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1450    /* win32 refactor */    /* win32 refactor */
1451    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
# Line 1040  void MeshAdapter::saveDX(const std::stri Line 1493  void MeshAdapter::saveDX(const std::stri
1493  // saves a data array in openVTK format:  // saves a data array in openVTK format:
1494  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1495  {  {
1496      int MAX_namelength=256;      unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;
1497      const int num_data=boost::python::extract<int>(arg.attr("__len__")());      const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1498    /* win32 refactor */    /* win32 refactor */
1499    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;    char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
# Line 1136  SystemMatrixAdapter MeshAdapter::newSyst Line 1589  SystemMatrixAdapter MeshAdapter::newSyst
1589      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1590      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1591  }  }
1592    // creates a TransportProblemAdapter
1593    TransportProblemAdapter MeshAdapter::newTransportProblem(
1594                          const double theta,
1595                          const int blocksize,
1596                          const escript::FunctionSpace& functionspace,
1597                          const int type) const
1598    {
1599        int reduceOrder=0;
1600        // is the domain right?
1601        const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());
1602        if (domain!=*this)
1603              throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1604        // is the function space type right
1605        if (functionspace.getTypeCode()==DegreesOfFreedom) {
1606            reduceOrder=0;
1607        } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
1608            reduceOrder=1;
1609        } else {
1610            throw FinleyAdapterException("Error - illegal function space type for system matrix rows.");
1611        }
1612        // generate matrix:
1613        
1614        Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1615        checkFinleyError();
1616        Paso_FCTransportProblem* transportProblem;
1617        transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);
1618        checkPasoError();
1619        Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1620        return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);
1621    }
1622    
1623  //  //
1624  // vtkObject MeshAdapter::createVtkObject() const  // vtkObject MeshAdapter::createVtkObject() const
# Line 1370  escript::Data MeshAdapter::getSize() con Line 1853  escript::Data MeshAdapter::getSize() con
1853    
1854  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1855  {  {
1856    int *out=0,i;    int *out = NULL;
1857    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
1858    switch (functionSpaceType) {    switch (functionSpaceType) {
1859      case(Nodes):      case(Nodes):

Legend:
Removed from v.1339  
changed lines
  Added in v.1455

  ViewVC Help
Powered by ViewVC 1.1.26