← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/src/CPPAdapter/MeshAdapter.cpp

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/src/CPPAdapter/MeshAdapter.cpp
 
1
 
 
2
/*******************************************************
 
3
*
 
4
* Copyright (c) 2003-2010 by University of Queensland
 
5
* Earth Systems Science Computational Center (ESSCC)
 
6
* http://www.uq.edu.au/esscc
 
7
*
 
8
* Primary Business: Queensland, Australia
 
9
* Licensed under the Open Software License version 3.0
 
10
* http://www.opensource.org/licenses/osl-3.0.php
 
11
*
 
12
*******************************************************/
 
13
 
 
14
 
 
15
#include "MeshAdapter.h"
 
16
#include "escript/Data.h"
 
17
#include "escript/DataFactory.h"
 
18
#ifdef USE_NETCDF
 
19
#include <netcdfcpp.h>
 
20
#endif
 
21
#ifdef ESYS_MPI
 
22
#include "esysUtils/Esys_MPI.h"
 
23
#endif
 
24
extern "C" {
 
25
#include "esysUtils/blocktimer.h"
 
26
}
 
27
 
 
28
using namespace std;
 
29
using namespace escript;
 
30
 
 
31
namespace dudley {
 
32
 
 
33
//
 
34
// define the static constants
 
35
MeshAdapter::FunctionSpaceNamesMapType MeshAdapter::m_functionSpaceTypeNames;
 
36
const int MeshAdapter::DegreesOfFreedom=DUDLEY_DEGREES_OF_FREEDOM;
 
37
const int MeshAdapter::ReducedDegreesOfFreedom=DUDLEY_REDUCED_DEGREES_OF_FREEDOM;
 
38
const int MeshAdapter::Nodes=DUDLEY_NODES;
 
39
const int MeshAdapter::ReducedNodes=DUDLEY_REDUCED_NODES;
 
40
const int MeshAdapter::Elements=DUDLEY_ELEMENTS;
 
41
const int MeshAdapter::ReducedElements=DUDLEY_REDUCED_ELEMENTS;
 
42
const int MeshAdapter::FaceElements=DUDLEY_FACE_ELEMENTS;
 
43
const int MeshAdapter::ReducedFaceElements=DUDLEY_REDUCED_FACE_ELEMENTS;
 
44
const int MeshAdapter::Points=DUDLEY_POINTS;
 
45
 
 
46
MeshAdapter::MeshAdapter(Dudley_Mesh* dudleyMesh)
 
47
{
 
48
   setFunctionSpaceTypeNames();
 
49
   //
 
50
   // need to use a null_deleter as Dudley_Mesh_free deletes the pointer
 
51
   // for us.
 
52
   m_dudleyMesh.reset(dudleyMesh,null_deleter());
 
53
}
 
54
 
 
55
//
 
56
// The copy constructor should just increment the use count
 
57
MeshAdapter::MeshAdapter(const MeshAdapter& in):
 
58
m_dudleyMesh(in.m_dudleyMesh)
 
59
{
 
60
   setFunctionSpaceTypeNames();
 
61
}
 
62
 
 
63
MeshAdapter::~MeshAdapter()
 
64
{
 
65
   //
 
66
   // I hope the case for the pointer being zero has been taken care of.
 
67
   //  cout << "In MeshAdapter destructor." << endl;
 
68
   if (m_dudleyMesh.unique()) {
 
69
      Dudley_Mesh_free(m_dudleyMesh.get());
 
70
   }
 
71
}
 
72
 
 
73
int MeshAdapter::getMPISize() const
 
74
{
 
75
   return m_dudleyMesh.get()->MPIInfo->size;
 
76
}
 
77
int MeshAdapter::getMPIRank() const
 
78
{
 
79
   return m_dudleyMesh.get()->MPIInfo->rank;
 
80
}
 
81
void MeshAdapter::MPIBarrier() const
 
82
{
 
83
#ifdef ESYS_MPI
 
84
   MPI_Barrier(m_dudleyMesh.get()->MPIInfo->comm);
 
85
#endif
 
86
   return;
 
87
}
 
88
bool MeshAdapter::onMasterProcessor() const
 
89
{
 
90
   return m_dudleyMesh.get()->MPIInfo->rank == 0;
 
91
}
 
92
 
 
93
 
 
94
#ifdef ESYS_MPI
 
95
  MPI_Comm
 
96
#else
 
97
  unsigned int
 
98
#endif
 
99
MeshAdapter::getMPIComm() const
 
100
{
 
101
#ifdef ESYS_MPI
 
102
        return m_dudleyMesh->MPIInfo->comm;
 
103
#else
 
104
        return 0;
 
105
#endif
 
106
}
 
107
 
 
108
 
 
109
Dudley_Mesh* MeshAdapter::getDudley_Mesh() const {
 
110
   return m_dudleyMesh.get();
 
111
}
 
112
 
 
113
void MeshAdapter::write(const string& fileName) const
 
114
{
 
115
   char *fName = (fileName.size()+1>0) ? TMPMEMALLOC(fileName.size()+1,char) : (char*)NULL;
 
116
   strcpy(fName,fileName.c_str());
 
117
   Dudley_Mesh_write(m_dudleyMesh.get(),fName);
 
118
   checkDudleyError();
 
119
   TMPMEMFREE(fName);
 
120
}
 
121
 
 
122
void MeshAdapter::Print_Mesh_Info(const bool full) const
 
123
{
 
124
   Dudley_PrintMesh_Info(m_dudleyMesh.get(), full);
 
125
}
 
126
 
 
127
void MeshAdapter::dump(const string& fileName) const
 
128
{
 
129
#ifdef USE_NETCDF
 
130
   const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
 
131
   NcVar *ids;
 
132
   int *int_ptr;
 
133
   Dudley_Mesh *mesh = m_dudleyMesh.get();
 
134
   Dudley_TagMap* tag_map;
 
135
   int num_Tags = 0;
 
136
   int mpi_size                         = mesh->MPIInfo->size;
 
137
   int mpi_rank                         = mesh->MPIInfo->rank;
 
138
   int numDim                           = mesh->Nodes->numDim;
 
139
   int numNodes                         = mesh->Nodes->numNodes;
 
140
   int num_Elements                     = mesh->Elements->numElements;
 
141
   int num_FaceElements                 = mesh->FaceElements->numElements;
 
142
   int num_Points                       = mesh->Points->numElements;
 
143
   int num_Elements_numNodes            = mesh->Elements->numNodes;
 
144
   int num_FaceElements_numNodes        = mesh->FaceElements->numNodes;
 
145
#ifdef ESYS_MPI
 
146
   MPI_Status status;
 
147
#endif
 
148
 
 
149
/* Incoming token indicates it's my turn to write */
 
150
#ifdef ESYS_MPI
 
151
   if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
 
152
#endif
 
153
 
 
154
   char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
 
155
                                                     mpi_size, mpi_rank);
 
156
 
 
157
   /* Figure out how much storage is required for tags */
 
158
   tag_map = mesh->TagMap;
 
159
   num_Tags = 0;
 
160
   while (tag_map) {
 
161
      num_Tags++;
 
162
      tag_map=tag_map->next;
 
163
   }
 
164
 
 
165
   // NetCDF error handler
 
166
   NcError err(NcError::verbose_nonfatal);
 
167
   // Create the file.
 
168
   NcFile dataFile(newFileName, NcFile::Replace);
 
169
   string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
 
170
   // check if writing was successful
 
171
   if (!dataFile.is_valid())
 
172
      throw DataException(msgPrefix+"Open file for output");
 
173
 
 
174
   // Define dimensions (num_Elements and dim_Elements are identical,
 
175
   // dim_Elements only appears if > 0)
 
176
   if (! (ncdims[0] = dataFile.add_dim("numNodes", numNodes)) )
 
177
      throw DataException(msgPrefix+"add_dim(numNodes)");
 
178
   if (! (ncdims[1] = dataFile.add_dim("numDim", numDim)) )
 
179
      throw DataException(msgPrefix+"add_dim(numDim)");
 
180
   if (! (ncdims[2] = dataFile.add_dim("mpi_size_plus_1", mpi_size+1)) )
 
181
      throw DataException(msgPrefix+"add_dim(mpi_size)");
 
182
   if (num_Elements>0)
 
183
      if (! (ncdims[3] = dataFile.add_dim("dim_Elements", num_Elements)) )
 
184
         throw DataException(msgPrefix+"add_dim(dim_Elements)");
 
185
   if (num_FaceElements>0)
 
186
      if (! (ncdims[4] = dataFile.add_dim("dim_FaceElements", num_FaceElements)) )
 
187
         throw DataException(msgPrefix+"add_dim(dim_FaceElements)");
 
188
   if (num_Points>0)
 
189
      if (! (ncdims[6] = dataFile.add_dim("dim_Points", num_Points)) )
 
190
         throw DataException(msgPrefix+"add_dim(dim_Points)");
 
191
   if (num_Elements>0)
 
192
      if (! (ncdims[7] = dataFile.add_dim("dim_Elements_Nodes", num_Elements_numNodes)) )
 
193
         throw DataException(msgPrefix+"add_dim(dim_Elements_Nodes)");
 
194
   if (num_FaceElements>0)
 
195
      if (! (ncdims[8] = dataFile.add_dim("dim_FaceElements_numNodes", num_FaceElements_numNodes)) )
 
196
         throw DataException(msgPrefix+"add_dim(dim_FaceElements_numNodes)");
 
197
   if (num_Tags>0)
 
198
      if (! (ncdims[10] = dataFile.add_dim("dim_Tags", num_Tags)) )
 
199
         throw DataException(msgPrefix+"add_dim(dim_Tags)");
 
200
 
 
201
   // Attributes: MPI size, MPI rank, Name, order, reduced_order
 
202
   if (!dataFile.add_att("mpi_size", mpi_size) )
 
203
      throw DataException(msgPrefix+"add_att(mpi_size)");
 
204
   if (!dataFile.add_att("mpi_rank", mpi_rank) )
 
205
      throw DataException(msgPrefix+"add_att(mpi_rank)");
 
206
   if (!dataFile.add_att("Name",mesh->Name) )
 
207
      throw DataException(msgPrefix+"add_att(Name)");
 
208
   if (!dataFile.add_att("numDim",numDim) )
 
209
      throw DataException(msgPrefix+"add_att(order)");
 
210
   if (!dataFile.add_att("order",mesh->integrationOrder) )
 
211
      throw DataException(msgPrefix+"add_att(order)");
 
212
   if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
 
213
      throw DataException(msgPrefix+"add_att(reduced_order)");
 
214
   if (!dataFile.add_att("numNodes",numNodes) )
 
215
      throw DataException(msgPrefix+"add_att(numNodes)");
 
216
   if (!dataFile.add_att("num_Elements",num_Elements) )
 
217
      throw DataException(msgPrefix+"add_att(num_Elements)");
 
218
   if (!dataFile.add_att("num_FaceElements",num_FaceElements) )
 
219
      throw DataException(msgPrefix+"add_att(num_FaceElements)");
 
220
   if (!dataFile.add_att("num_Points",num_Points) )
 
221
      throw DataException(msgPrefix+"add_att(num_Points)");
 
222
   if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )
 
223
      throw DataException(msgPrefix+"add_att(num_Elements_numNodes)");
 
224
   if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
 
225
      throw DataException(msgPrefix+"add_att(num_FaceElements_numNodes)");
 
226
   if (!dataFile.add_att("Elements_TypeId", mesh->Elements->etype) )
 
227
      throw DataException(msgPrefix+"add_att(Elements_TypeId)");
 
228
   if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->etype) )
 
229
      throw DataException(msgPrefix+"add_att(FaceElements_TypeId)");
 
230
   if (!dataFile.add_att("Points_TypeId", mesh->Points->etype) )
 
231
      throw DataException(msgPrefix+"add_att(Points_TypeId)");
 
232
   if (!dataFile.add_att("num_Tags", num_Tags) )
 
233
      throw DataException(msgPrefix+"add_att(num_Tags)");
 
234
 
 
235
   // // // // // Nodes // // // // //
 
236
 
 
237
   // Nodes nodeDistribution
 
238
   if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
 
239
      throw DataException(msgPrefix+"add_var(Nodes_NodeDistribution)");
 
240
   int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
 
241
   if (! (ids->put(int_ptr, mpi_size+1)) )
 
242
      throw DataException(msgPrefix+"put(Nodes_NodeDistribution)");
 
243
 
 
244
   // Nodes degreesOfFreedomDistribution
 
245
   if (! ( ids = dataFile.add_var("Nodes_DofDistribution", ncInt, ncdims[2])) )
 
246
      throw DataException(msgPrefix+"add_var(Nodes_DofDistribution)");
 
247
   int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
 
248
   if (! (ids->put(int_ptr, mpi_size+1)) )
 
249
      throw DataException(msgPrefix+"put(Nodes_DofDistribution)");
 
250
 
 
251
   // Only write nodes if non-empty because NetCDF doesn't like empty arrays
 
252
   // (it treats them as NC_UNLIMITED)
 
253
   if (numNodes>0) {
 
254
 
 
255
      // Nodes Id
 
256
      if (! ( ids = dataFile.add_var("Nodes_Id", ncInt, ncdims[0])) )
 
257
         throw DataException(msgPrefix+"add_var(Nodes_Id)");
 
258
      int_ptr = &mesh->Nodes->Id[0];
 
259
      if (! (ids->put(int_ptr, numNodes)) )
 
260
         throw DataException(msgPrefix+"put(Nodes_Id)");
 
261
 
 
262
      // Nodes Tag
 
263
      if (! ( ids = dataFile.add_var("Nodes_Tag", ncInt, ncdims[0])) )
 
264
         throw DataException(msgPrefix+"add_var(Nodes_Tag)");
 
265
      int_ptr = &mesh->Nodes->Tag[0];
 
266
      if (! (ids->put(int_ptr, numNodes)) )
 
267
         throw DataException(msgPrefix+"put(Nodes_Tag)");
 
268
 
 
269
      // Nodes gDOF
 
270
      if (! ( ids = dataFile.add_var("Nodes_gDOF", ncInt, ncdims[0])) )
 
271
         throw DataException(msgPrefix+"add_var(Nodes_gDOF)");
 
272
      int_ptr = &mesh->Nodes->globalDegreesOfFreedom[0];
 
273
      if (! (ids->put(int_ptr, numNodes)) )
 
274
         throw DataException(msgPrefix+"put(Nodes_gDOF)");
 
275
 
 
276
      // Nodes global node index
 
277
      if (! ( ids = dataFile.add_var("Nodes_gNI", ncInt, ncdims[0])) )
 
278
         throw DataException(msgPrefix+"add_var(Nodes_gNI)");
 
279
      int_ptr = &mesh->Nodes->globalNodesIndex[0];
 
280
      if (! (ids->put(int_ptr, numNodes)) )
 
281
         throw DataException(msgPrefix+"put(Nodes_gNI)");
 
282
 
 
283
      // Nodes grDof
 
284
      if (! ( ids = dataFile.add_var("Nodes_grDfI", ncInt, ncdims[0])) )
 
285
         throw DataException(msgPrefix+"add_var(Nodes_grDfI)");
 
286
      int_ptr = &mesh->Nodes->globalReducedDOFIndex[0];
 
287
      if (! (ids->put(int_ptr, numNodes)) )
 
288
         throw DataException(msgPrefix+"put(Nodes_grDfI)");
 
289
 
 
290
      // Nodes grNI
 
291
      if (! ( ids = dataFile.add_var("Nodes_grNI", ncInt, ncdims[0])) )
 
292
         throw DataException(msgPrefix+"add_var(Nodes_grNI)");
 
293
      int_ptr = &mesh->Nodes->globalReducedNodesIndex[0];
 
294
      if (! (ids->put(int_ptr, numNodes)) )
 
295
         throw DataException(msgPrefix+"put(Nodes_grNI)");
 
296
 
 
297
      // Nodes Coordinates
 
298
      if (! ( ids = dataFile.add_var("Nodes_Coordinates", ncDouble, ncdims[0], ncdims[1]) ) )
 
299
         throw DataException(msgPrefix+"add_var(Nodes_Coordinates)");
 
300
      if (! (ids->put(&(mesh->Nodes->Coordinates[INDEX2(0,0,numDim)]), numNodes, numDim)) )
 
301
         throw DataException(msgPrefix+"put(Nodes_Coordinates)");
 
302
 
 
303
   }
 
304
 
 
305
   // // // // // Elements // // // // //
 
306
 
 
307
   if (num_Elements>0) {
 
308
 
 
309
      // Elements_Id
 
310
      if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
 
311
         throw DataException(msgPrefix+"add_var(Elements_Id)");
 
312
      int_ptr = &mesh->Elements->Id[0];
 
313
      if (! (ids->put(int_ptr, num_Elements)) )
 
314
         throw DataException(msgPrefix+"put(Elements_Id)");
 
315
 
 
316
      // Elements_Tag
 
317
      if (! ( ids = dataFile.add_var("Elements_Tag", ncInt, ncdims[3])) )
 
318
         throw DataException(msgPrefix+"add_var(Elements_Tag)");
 
319
      int_ptr = &mesh->Elements->Tag[0];
 
320
      if (! (ids->put(int_ptr, num_Elements)) )
 
321
         throw DataException(msgPrefix+"put(Elements_Tag)");
 
322
 
 
323
      // Elements_Owner
 
324
      if (! ( ids = dataFile.add_var("Elements_Owner", ncInt, ncdims[3])) )
 
325
         throw DataException(msgPrefix+"add_var(Elements_Owner)");
 
326
      int_ptr = &mesh->Elements->Owner[0];
 
327
      if (! (ids->put(int_ptr, num_Elements)) )
 
328
         throw DataException(msgPrefix+"put(Elements_Owner)");
 
329
 
 
330
      // Elements_Color
 
331
      if (! ( ids = dataFile.add_var("Elements_Color", ncInt, ncdims[3])) )
 
332
         throw DataException(msgPrefix+"add_var(Elements_Color)");
 
333
      int_ptr = &mesh->Elements->Color[0];
 
334
      if (! (ids->put(int_ptr, num_Elements)) )
 
335
         throw DataException(msgPrefix+"put(Elements_Color)");
 
336
 
 
337
      // Elements_Nodes
 
338
      if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
 
339
         throw DataException(msgPrefix+"add_var(Elements_Nodes)");
 
340
      if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
 
341
         throw DataException(msgPrefix+"put(Elements_Nodes)");
 
342
 
 
343
   }
 
344
 
 
345
   // // // // // Face_Elements // // // // //
 
346
 
 
347
   if (num_FaceElements>0) {
 
348
 
 
349
      // FaceElements_Id
 
350
      if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
 
351
         throw DataException(msgPrefix+"add_var(FaceElements_Id)");
 
352
      int_ptr = &mesh->FaceElements->Id[0];
 
353
      if (! (ids->put(int_ptr, num_FaceElements)) )
 
354
         throw DataException(msgPrefix+"put(FaceElements_Id)");
 
355
 
 
356
      // FaceElements_Tag
 
357
      if (! ( ids = dataFile.add_var("FaceElements_Tag", ncInt, ncdims[4])) )
 
358
         throw DataException(msgPrefix+"add_var(FaceElements_Tag)");
 
359
      int_ptr = &mesh->FaceElements->Tag[0];
 
360
      if (! (ids->put(int_ptr, num_FaceElements)) )
 
361
         throw DataException(msgPrefix+"put(FaceElements_Tag)");
 
362
 
 
363
      // FaceElements_Owner
 
364
      if (! ( ids = dataFile.add_var("FaceElements_Owner", ncInt, ncdims[4])) )
 
365
         throw DataException(msgPrefix+"add_var(FaceElements_Owner)");
 
366
      int_ptr = &mesh->FaceElements->Owner[0];
 
367
      if (! (ids->put(int_ptr, num_FaceElements)) )
 
368
         throw DataException(msgPrefix+"put(FaceElements_Owner)");
 
369
 
 
370
      // FaceElements_Color
 
371
      if (! ( ids = dataFile.add_var("FaceElements_Color", ncInt, ncdims[4])) )
 
372
         throw DataException(msgPrefix+"add_var(FaceElements_Color)");
 
373
      int_ptr = &mesh->FaceElements->Color[0];
 
374
      if (! (ids->put(int_ptr, num_FaceElements)) )
 
375
         throw DataException(msgPrefix+"put(FaceElements_Color)");
 
376
 
 
377
      // FaceElements_Nodes
 
378
      if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
 
379
         throw DataException(msgPrefix+"add_var(FaceElements_Nodes)");
 
380
      if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
 
381
         throw DataException(msgPrefix+"put(FaceElements_Nodes)");
 
382
 
 
383
   }
 
384
 
 
385
   // // // // // Points // // // // //
 
386
 
 
387
   if (num_Points>0) {
 
388
 
 
389
      fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
 
390
 
 
391
      // Points_Id
 
392
      if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
 
393
         throw DataException(msgPrefix+"add_var(Points_Id)");
 
394
      int_ptr = &mesh->Points->Id[0];
 
395
      if (! (ids->put(int_ptr, num_Points)) )
 
396
         throw DataException(msgPrefix+"put(Points_Id)");
 
397
 
 
398
      // Points_Tag
 
399
      if (! ( ids = dataFile.add_var("Points_Tag", ncInt, ncdims[6])) )
 
400
         throw DataException(msgPrefix+"add_var(Points_Tag)");
 
401
      int_ptr = &mesh->Points->Tag[0];
 
402
      if (! (ids->put(int_ptr, num_Points)) )
 
403
         throw DataException(msgPrefix+"put(Points_Tag)");
 
404
 
 
405
      // Points_Owner
 
406
      if (! ( ids = dataFile.add_var("Points_Owner", ncInt, ncdims[6])) )
 
407
         throw DataException(msgPrefix+"add_var(Points_Owner)");
 
408
      int_ptr = &mesh->Points->Owner[0];
 
409
      if (! (ids->put(int_ptr, num_Points)) )
 
410
         throw DataException(msgPrefix+"put(Points_Owner)");
 
411
 
 
412
      // Points_Color
 
413
      if (! ( ids = dataFile.add_var("Points_Color", ncInt, ncdims[6])) )
 
414
         throw DataException(msgPrefix+"add_var(Points_Color)");
 
415
      int_ptr = &mesh->Points->Color[0];
 
416
      if (! (ids->put(int_ptr, num_Points)) )
 
417
         throw DataException(msgPrefix+"put(Points_Color)");
 
418
 
 
419
      // Points_Nodes
 
420
      // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
 
421
      if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
 
422
         throw DataException(msgPrefix+"add_var(Points_Nodes)");
 
423
      if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
 
424
         throw DataException(msgPrefix+"put(Points_Nodes)");
 
425
 
 
426
   }
 
427
 
 
428
   // // // // // TagMap // // // // //
 
429
 
 
430
   if (num_Tags>0) {
 
431
 
 
432
      // Temp storage to gather node IDs
 
433
      int *Tags_keys = TMPMEMALLOC(num_Tags, int);
 
434
      char name_temp[4096];
 
435
 
 
436
      /* Copy tag data into temp arrays */
 
437
      tag_map = mesh->TagMap;
 
438
      if (tag_map) {
 
439
         int i = 0;
 
440
         while (tag_map) {
 
441
            Tags_keys[i++] = tag_map->tag_key;
 
442
            tag_map=tag_map->next;
 
443
         }
 
444
      }
 
445
 
 
446
      // Tags_keys
 
447
      if (! ( ids = dataFile.add_var("Tags_keys", ncInt, ncdims[10])) )
 
448
         throw DataException(msgPrefix+"add_var(Tags_keys)");
 
449
      int_ptr = &Tags_keys[0];
 
450
      if (! (ids->put(int_ptr, num_Tags)) )
 
451
         throw DataException(msgPrefix+"put(Tags_keys)");
 
452
 
 
453
      // Tags_names_*
 
454
      // This is an array of strings, it should be stored as an array but
 
455
      // instead I have hacked in one attribute per string because the NetCDF
 
456
      // manual doesn't tell how to do an array of strings
 
457
      tag_map = mesh->TagMap;
 
458
      if (tag_map) {
 
459
         int i = 0;
 
460
         while (tag_map) {
 
461
            sprintf(name_temp, "Tags_name_%d", i);
 
462
            if (!dataFile.add_att(name_temp, tag_map->name) )
 
463
               throw DataException(msgPrefix+"add_att(Tags_names_XX)");
 
464
            tag_map=tag_map->next;
 
465
            i++;
 
466
         }
 
467
      }
 
468
 
 
469
      TMPMEMFREE(Tags_keys);
 
470
   }
 
471
 
 
472
/* Send token to next MPI process so he can take his turn */
 
473
#ifdef ESYS_MPI
 
474
   if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
 
475
#endif
 
476
 
 
477
   // NetCDF file is closed by destructor of NcFile object
 
478
 
 
479
#else
 
480
   Dudley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
 
481
#endif  /* USE_NETCDF */
 
482
   checkDudleyError();
 
483
}
 
484
 
 
485
string MeshAdapter::getDescription() const
 
486
{
 
487
   return "DudleyMesh";
 
488
}
 
489
 
 
490
string MeshAdapter::functionSpaceTypeAsString(int functionSpaceType) const
 
491
{
 
492
   FunctionSpaceNamesMapType::iterator loc;
 
493
   loc=m_functionSpaceTypeNames.find(functionSpaceType);
 
494
   if (loc==m_functionSpaceTypeNames.end()) {
 
495
      return "Invalid function space type code.";
 
496
   } else {
 
497
      return loc->second;
 
498
   }
 
499
}
 
500
 
 
501
bool MeshAdapter::isValidFunctionSpaceType(int functionSpaceType) const
 
502
{
 
503
   FunctionSpaceNamesMapType::iterator loc;
 
504
   loc=m_functionSpaceTypeNames.find(functionSpaceType);
 
505
   return (loc!=m_functionSpaceTypeNames.end());
 
506
}
 
507
 
 
508
void MeshAdapter::setFunctionSpaceTypeNames()
 
509
{
 
510
   m_functionSpaceTypeNames.insert
 
511
   (FunctionSpaceNamesMapType::value_type(DegreesOfFreedom,"Dudley_DegreesOfFreedom"));
 
512
   m_functionSpaceTypeNames.insert
 
513
   (FunctionSpaceNamesMapType::value_type(ReducedDegreesOfFreedom,"Dudley_ReducedDegreesOfFreedom"));
 
514
   m_functionSpaceTypeNames.insert
 
515
   (FunctionSpaceNamesMapType::value_type(Nodes,"Dudley_Nodes"));
 
516
   m_functionSpaceTypeNames.insert
 
517
   (FunctionSpaceNamesMapType::value_type(ReducedNodes,"Dudley_Reduced_Nodes"));
 
518
   m_functionSpaceTypeNames.insert
 
519
   (FunctionSpaceNamesMapType::value_type(Elements,"Dudley_Elements"));
 
520
   m_functionSpaceTypeNames.insert
 
521
   (FunctionSpaceNamesMapType::value_type(ReducedElements,"Dudley_Reduced_Elements"));
 
522
   m_functionSpaceTypeNames.insert
 
523
   (FunctionSpaceNamesMapType::value_type(FaceElements,"Dudley_Face_Elements"));
 
524
   m_functionSpaceTypeNames.insert
 
525
   (FunctionSpaceNamesMapType::value_type(ReducedFaceElements,"Dudley_Reduced_Face_Elements"));
 
526
   m_functionSpaceTypeNames.insert
 
527
   (FunctionSpaceNamesMapType::value_type(Points,"Dudley_Points"));
 
528
}
 
529
 
 
530
int MeshAdapter::getContinuousFunctionCode() const
 
531
{
 
532
   return Nodes;
 
533
}
 
534
int MeshAdapter::getReducedContinuousFunctionCode() const
 
535
{
 
536
   return ReducedNodes;
 
537
}
 
538
 
 
539
int MeshAdapter::getFunctionCode() const
 
540
{
 
541
   return Elements;
 
542
}
 
543
int MeshAdapter::getReducedFunctionCode() const
 
544
{
 
545
   return ReducedElements;
 
546
}
 
547
 
 
548
int MeshAdapter::getFunctionOnBoundaryCode() const
 
549
{
 
550
   return FaceElements;
 
551
}
 
552
int MeshAdapter::getReducedFunctionOnBoundaryCode() const
 
553
{
 
554
   return ReducedFaceElements;
 
555
}
 
556
 
 
557
int MeshAdapter::getFunctionOnContactZeroCode() const
 
558
{
 
559
   throw DudleyAdapterException("Dudley does not support contact elements.");
 
560
}
 
561
 
 
562
int MeshAdapter::getReducedFunctionOnContactZeroCode() const
 
563
{
 
564
   throw DudleyAdapterException("Dudley does not support contact elements.");
 
565
}
 
566
 
 
567
int MeshAdapter::getFunctionOnContactOneCode() const
 
568
{
 
569
   throw DudleyAdapterException("Dudley does not support contact elements.");
 
570
}
 
571
 
 
572
int MeshAdapter::getReducedFunctionOnContactOneCode() const
 
573
{
 
574
   throw DudleyAdapterException("Dudley does not support contact elements.");
 
575
}
 
576
 
 
577
int MeshAdapter::getSolutionCode() const
 
578
{
 
579
   return DegreesOfFreedom;
 
580
}
 
581
 
 
582
int MeshAdapter::getReducedSolutionCode() const
 
583
{
 
584
   return ReducedDegreesOfFreedom;
 
585
}
 
586
 
 
587
int MeshAdapter::getDiracDeltaFunctionCode() const
 
588
{
 
589
   return Points;
 
590
}
 
591
 
 
592
//
 
593
// return the spatial dimension of the Mesh:
 
594
//
 
595
int MeshAdapter::getDim() const
 
596
{
 
597
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
598
   int numDim=Dudley_Mesh_getDim(mesh);
 
599
   checkDudleyError();
 
600
   return numDim;
 
601
}
 
602
 
 
603
//
 
604
// Return the number of data points summed across all MPI processes
 
605
//
 
606
int MeshAdapter::getNumDataPointsGlobal() const
 
607
{
 
608
   return Dudley_NodeFile_getGlobalNumNodes(m_dudleyMesh.get()->Nodes);
 
609
}
 
610
 
 
611
//
 
612
// return the number of data points per sample and the number of samples
 
613
// needed to represent data on a parts of the mesh.
 
614
//
 
615
pair<int,int> MeshAdapter::getDataShape(int functionSpaceCode) const
 
616
{
 
617
   int numDataPointsPerSample=0;
 
618
   int numSamples=0;
 
619
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
620
   switch (functionSpaceCode) {
 
621
   case(Nodes):
 
622
   numDataPointsPerSample=1;
 
623
   numSamples=Dudley_NodeFile_getNumNodes(mesh->Nodes);
 
624
   break;
 
625
   case(ReducedNodes):
 
626
   numDataPointsPerSample=1;
 
627
   numSamples=Dudley_NodeFile_getNumReducedNodes(mesh->Nodes);
 
628
   break;
 
629
   case(Elements):
 
630
   if (mesh->Elements!=NULL) {
 
631
      numSamples=mesh->Elements->numElements;
 
632
      numDataPointsPerSample=mesh->Elements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
 
633
   }
 
634
   break;
 
635
   case(ReducedElements):
 
636
   if (mesh->Elements!=NULL) {
 
637
      numSamples=mesh->Elements->numElements;
 
638
      numDataPointsPerSample=(mesh->Elements->numLocalDim==0)?0:1;
 
639
   }
 
640
   break;
 
641
   case(FaceElements):
 
642
   if (mesh->FaceElements!=NULL) {
 
643
      numDataPointsPerSample=mesh->FaceElements->numLocalDim+1/*referenceElementSet->referenceElement->BasisFunctions->numQuadNodes*/;
 
644
      numSamples=mesh->FaceElements->numElements;
 
645
   }
 
646
   break;
 
647
   case(ReducedFaceElements):
 
648
   if (mesh->FaceElements!=NULL) {
 
649
      numDataPointsPerSample=(mesh->FaceElements->numLocalDim==0)?0:1/*referenceElementSet->referenceElementReducedQuadrature->BasisFunctions->numQuadNodes*/;
 
650
      numSamples=mesh->FaceElements->numElements;
 
651
   }
 
652
   break;
 
653
   case(Points):
 
654
   if (mesh->Points!=NULL) {
 
655
      numDataPointsPerSample=1;
 
656
      numSamples=mesh->Points->numElements;
 
657
   }
 
658
   break;
 
659
   case(DegreesOfFreedom):
 
660
   if (mesh->Nodes!=NULL) {
 
661
      numDataPointsPerSample=1;
 
662
      numSamples=Dudley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
 
663
   }
 
664
   break;
 
665
   case(ReducedDegreesOfFreedom):
 
666
   if (mesh->Nodes!=NULL) {
 
667
      numDataPointsPerSample=1;
 
668
      numSamples=Dudley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
 
669
   }
 
670
   break;
 
671
   default:
 
672
      stringstream temp;
 
673
      temp << "Error - Invalid function space type: " << functionSpaceCode << " for domain: " << getDescription();
 
674
      throw DudleyAdapterException(temp.str());
 
675
      break;
 
676
   }
 
677
   return pair<int,int>(numDataPointsPerSample,numSamples);
 
678
}
 
679
 
 
680
//
 
681
// adds linear PDE of second order into a given stiffness matrix and right hand side:
 
682
//
 
683
void MeshAdapter::addPDEToSystem(
 
684
                                 SystemMatrixAdapter& mat, escript::Data& rhs,
 
685
                                 const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
 
686
                                 const escript::Data& d, const escript::Data& y) const
 
687
{
 
688
   escriptDataC _rhs=rhs.getDataC();
 
689
   escriptDataC _A  =A.getDataC();
 
690
   escriptDataC _B=B.getDataC();
 
691
   escriptDataC _C=C.getDataC();
 
692
   escriptDataC _D=D.getDataC();
 
693
   escriptDataC _X=X.getDataC();
 
694
   escriptDataC _Y=Y.getDataC();
 
695
   escriptDataC _d=d.getDataC();
 
696
   escriptDataC _y=y.getDataC();
 
697
 
 
698
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
699
 
 
700
   Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
 
701
   checkDudleyError();
 
702
 
 
703
 
 
704
   Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
 
705
   checkDudleyError();
 
706
 
 
707
 
 
708
   checkDudleyError();
 
709
}
 
710
 
 
711
void  MeshAdapter::addPDEToLumpedSystem(
 
712
                                        escript::Data& mat,
 
713
                                        const escript::Data& D,
 
714
                                        const escript::Data& d) const
 
715
{
 
716
   escriptDataC _mat=mat.getDataC();
 
717
   escriptDataC _D=D.getDataC();
 
718
   escriptDataC _d=d.getDataC();
 
719
 
 
720
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
721
 
 
722
   Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
 
723
   checkDudleyError();
 
724
   
 
725
   Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
 
726
   checkDudleyError();
 
727
 
 
728
}
 
729
 
 
730
 
 
731
//
 
732
// adds linear PDE of second order into the right hand side only
 
733
//
 
734
void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y) const
 
735
{
 
736
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
737
 
 
738
   escriptDataC _rhs=rhs.getDataC();
 
739
   escriptDataC _X=X.getDataC();
 
740
   escriptDataC _Y=Y.getDataC();
 
741
   escriptDataC _y=y.getDataC();
 
742
 
 
743
   Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
 
744
   checkDudleyError();
 
745
 
 
746
   Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y );
 
747
   checkDudleyError();
 
748
 
 
749
   checkDudleyError();
 
750
}
 
751
//
 
752
// adds PDE of second order into a transport problem
 
753
//
 
754
void MeshAdapter::addPDEToTransportProblem(
 
755
                                           TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,
 
756
                                           const escript::Data& A, const escript::Data& B, const escript::Data& C,
 
757
                                           const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
 
758
                                           const escript::Data& d, const escript::Data& y, 
 
759
                                           const escript::Data& d_contact,const escript::Data& y_contact) const
 
760
{
 
761
   DataTypes::ShapeType shape;
 
762
   source.expand();
 
763
   escriptDataC _source=source.getDataC();
 
764
   escriptDataC _M=M.getDataC();
 
765
   escriptDataC _A=A.getDataC();
 
766
   escriptDataC _B=B.getDataC();
 
767
   escriptDataC _C=C.getDataC();
 
768
   escriptDataC _D=D.getDataC();
 
769
   escriptDataC _X=X.getDataC();
 
770
   escriptDataC _Y=Y.getDataC();
 
771
   escriptDataC _d=d.getDataC();
 
772
   escriptDataC _y=y.getDataC();
 
773
   escriptDataC _d_contact=d_contact.getDataC();
 
774
   escriptDataC _y_contact=y_contact.getDataC();
 
775
 
 
776
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
777
   Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();
 
778
 
 
779
   Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
 
780
   checkDudleyError();
 
781
 
 
782
   Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &_source, &_A, &_B, &_C, &_D, &_X, &_Y );
 
783
   checkDudleyError();
 
784
 
 
785
   Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &_source, 0, 0, 0, &_d, 0, &_y );
 
786
   checkDudleyError();
 
787
 
 
788
   checkDudleyError();
 
789
}
 
790
 
 
791
//
 
792
// interpolates data between different function spaces:
 
793
//
 
794
void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
 
795
{
 
796
   const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
 
797
   const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
 
798
   if (inDomain!=*this)  
 
799
      throw DudleyAdapterException("Error - Illegal domain of interpolant.");
 
800
   if (targetDomain!=*this) 
 
801
      throw DudleyAdapterException("Error - Illegal domain of interpolation target.");
 
802
 
 
803
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
804
   escriptDataC _target=target.getDataC();
 
805
   escriptDataC _in=in.getDataC();
 
806
   switch(in.getFunctionSpace().getTypeCode()) {
 
807
   case(Nodes):
 
808
      switch(target.getFunctionSpace().getTypeCode()) {
 
809
      case(Nodes):
 
810
      case(ReducedNodes):
 
811
      case(DegreesOfFreedom):
 
812
      case(ReducedDegreesOfFreedom):
 
813
      Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
814
      break;
 
815
      case(Elements):
 
816
      case(ReducedElements):
 
817
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
 
818
      break;
 
819
      case(FaceElements):
 
820
      case(ReducedFaceElements):
 
821
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
 
822
      break;
 
823
      case(Points):
 
824
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
 
825
      break;
 
826
      default:
 
827
         stringstream temp;
 
828
         temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
 
829
         throw DudleyAdapterException(temp.str());
 
830
         break;
 
831
      }
 
832
      break;
 
833
   case(ReducedNodes):
 
834
      switch(target.getFunctionSpace().getTypeCode()) {
 
835
      case(Nodes):
 
836
      case(ReducedNodes):
 
837
      case(DegreesOfFreedom):
 
838
      case(ReducedDegreesOfFreedom):
 
839
      Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
840
      break;
 
841
      case(Elements):
 
842
      case(ReducedElements):
 
843
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
 
844
      break;
 
845
      case(FaceElements):
 
846
      case(ReducedFaceElements):
 
847
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
 
848
      break;
 
849
      case(Points):
 
850
      Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
 
851
      break;
 
852
      default:
 
853
         stringstream temp;
 
854
         temp << "Error - Interpolation on Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
 
855
         throw DudleyAdapterException(temp.str());
 
856
         break;
 
857
      }
 
858
      break;
 
859
   case(Elements):
 
860
      if (target.getFunctionSpace().getTypeCode()==Elements) {
 
861
         Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
 
862
      } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
 
863
         Dudley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
 
864
      } else {
 
865
         throw DudleyAdapterException("Error - No interpolation with data on elements possible.");
 
866
      }
 
867
      break;
 
868
   case(ReducedElements):
 
869
      if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
 
870
         Dudley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
 
871
      } else {
 
872
         throw DudleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
 
873
      }
 
874
      break;
 
875
   case(FaceElements):
 
876
      if (target.getFunctionSpace().getTypeCode()==FaceElements) {
 
877
         Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
 
878
      } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
 
879
         Dudley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
 
880
      } else {
 
881
         throw DudleyAdapterException("Error - No interpolation with data on face elements possible.");
 
882
      }
 
883
      break;
 
884
   case(ReducedFaceElements):
 
885
      if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
 
886
         Dudley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
 
887
      } else {
 
888
         throw DudleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
 
889
      }
 
890
      break;
 
891
   case(Points):
 
892
      if (target.getFunctionSpace().getTypeCode()==Points) {
 
893
         Dudley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
 
894
      } else {
 
895
         throw DudleyAdapterException("Error - No interpolation with data on points possible.");
 
896
      }
 
897
      break;
 
898
   case(DegreesOfFreedom):      
 
899
      switch(target.getFunctionSpace().getTypeCode()) {
 
900
      case(ReducedDegreesOfFreedom):
 
901
      case(DegreesOfFreedom):
 
902
      Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
903
      break;
 
904
   
 
905
      case(Nodes):
 
906
      case(ReducedNodes):
 
907
      if (getMPISize()>1) {
 
908
         escript::Data temp=escript::Data(in);
 
909
         temp.expand();
 
910
         escriptDataC _in2 = temp.getDataC();
 
911
         Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
 
912
      } else {
 
913
         Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
914
      }
 
915
      break;
 
916
      case(Elements):
 
917
      case(ReducedElements):
 
918
      if (getMPISize()>1) {
 
919
         escript::Data temp=escript::Data( in,  continuousFunction(*this) );
 
920
         escriptDataC _in2 = temp.getDataC();
 
921
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
 
922
      } else {
 
923
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
 
924
      }
 
925
      break;
 
926
      case(FaceElements):
 
927
      case(ReducedFaceElements):
 
928
      if (getMPISize()>1) {
 
929
         escript::Data temp=escript::Data( in,  continuousFunction(*this) );
 
930
         escriptDataC _in2 = temp.getDataC();
 
931
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
 
932
   
 
933
      } else {
 
934
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
 
935
      }
 
936
      break;
 
937
      case(Points):
 
938
      if (getMPISize()>1) {
 
939
         escript::Data temp=escript::Data( in,  continuousFunction(*this) );
 
940
         escriptDataC _in2 = temp.getDataC();
 
941
      } else {
 
942
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
 
943
      }
 
944
      break;
 
945
      default:
 
946
         stringstream temp;
 
947
         temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
 
948
         throw DudleyAdapterException(temp.str());
 
949
         break;
 
950
      }
 
951
      break;
 
952
   case(ReducedDegreesOfFreedom):
 
953
      switch(target.getFunctionSpace().getTypeCode()) {
 
954
      case(Nodes):
 
955
      throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to mesh nodes.");
 
956
      break;
 
957
      case(ReducedNodes):
 
958
      if (getMPISize()>1) {
 
959
         escript::Data temp=escript::Data(in);
 
960
         temp.expand();
 
961
         escriptDataC _in2 = temp.getDataC();
 
962
         Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
 
963
      } else {
 
964
         Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
965
      }
 
966
      break;
 
967
      case(DegreesOfFreedom):
 
968
      throw DudleyAdapterException("Error - Dudley does not support interpolation from reduced degrees of freedom to degrees of freedom");
 
969
      break;
 
970
      case(ReducedDegreesOfFreedom):
 
971
      Dudley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
 
972
      break;
 
973
      case(Elements):
 
974
      case(ReducedElements):
 
975
      if (getMPISize()>1) {
 
976
         escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
 
977
         escriptDataC _in2 = temp.getDataC();
 
978
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
 
979
      } else {
 
980
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
 
981
      }
 
982
      break;
 
983
      case(FaceElements):
 
984
      case(ReducedFaceElements):
 
985
      if (getMPISize()>1) {
 
986
         escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
 
987
         escriptDataC _in2 = temp.getDataC();
 
988
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
 
989
      } else {
 
990
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
 
991
      }
 
992
      break;
 
993
      case(Points):
 
994
      if (getMPISize()>1) {
 
995
         escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
 
996
         escriptDataC _in2 = temp.getDataC();
 
997
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
 
998
      } else {
 
999
         Dudley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
 
1000
      }
 
1001
      break;
 
1002
      default:
 
1003
         stringstream temp;
 
1004
         temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
 
1005
         throw DudleyAdapterException(temp.str());
 
1006
         break;
 
1007
      }
 
1008
      break;
 
1009
   default:
 
1010
      stringstream temp;
 
1011
      temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
 
1012
      throw DudleyAdapterException(temp.str());
 
1013
      break;
 
1014
   }
 
1015
   checkDudleyError();
 
1016
}
 
1017
 
 
1018
//
 
1019
// copies the locations of sample points into x:
 
1020
//
 
1021
void MeshAdapter::setToX(escript::Data& arg) const
 
1022
{
 
1023
   const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
 
1024
   if (argDomain!=*this) 
 
1025
      throw DudleyAdapterException("Error - Illegal domain of data point locations");
 
1026
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1027
   // in case of values node coordinates we can do the job directly:
 
1028
   if (arg.getFunctionSpace().getTypeCode()==Nodes) {
 
1029
      escriptDataC _arg=arg.getDataC();
 
1030
      Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
 
1031
   } else {
 
1032
      escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
 
1033
      escriptDataC _tmp_data=tmp_data.getDataC();
 
1034
      Dudley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
 
1035
      // this is then interpolated onto arg:
 
1036
      interpolateOnDomain(arg,tmp_data);
 
1037
   }
 
1038
   checkDudleyError();
 
1039
}
 
1040
 
 
1041
//
 
1042
// return the normal vectors at the location of data points as a Data object:
 
1043
//
 
1044
void MeshAdapter::setToNormal(escript::Data& normal) const
 
1045
{
 
1046
/*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
 
1047
   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
 
1048
   if (normalDomain!=*this) 
 
1049
      throw DudleyAdapterException("Error - Illegal domain of normal locations");
 
1050
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1051
   escriptDataC _normal=normal.getDataC();
 
1052
   switch(normal.getFunctionSpace().getTypeCode()) {
 
1053
   case(Nodes):
 
1054
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for nodes");
 
1055
   break;
 
1056
   case(ReducedNodes):
 
1057
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced nodes");
 
1058
   break;
 
1059
   case(Elements):
 
1060
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements");
 
1061
   break;
 
1062
   case(ReducedElements):
 
1063
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for elements with reduced integration order");
 
1064
   break;
 
1065
   case (FaceElements):
 
1066
   Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
 
1067
   break;
 
1068
   case (ReducedFaceElements):
 
1069
   Dudley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);
 
1070
   break;
 
1071
   case(Points):
 
1072
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for point elements");
 
1073
   break;
 
1074
   case(DegreesOfFreedom):
 
1075
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for degrees of freedom.");
 
1076
   break;
 
1077
   case(ReducedDegreesOfFreedom):
 
1078
   throw DudleyAdapterException("Error - Dudley does not support surface normal vectors for reduced degrees of freedom.");
 
1079
   break;
 
1080
   default:
 
1081
      stringstream temp;
 
1082
      temp << "Error - Normal Vectors: Dudley does not know anything about function space type " << normal.getFunctionSpace().getTypeCode();
 
1083
      throw DudleyAdapterException(temp.str());
 
1084
      break;
 
1085
   }
 
1086
   checkDudleyError();
 
1087
}
 
1088
 
 
1089
//
 
1090
// interpolates data to other domain:
 
1091
//
 
1092
void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
 
1093
{
 
1094
   const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
 
1095
   const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
 
1096
   if (targetDomain!=this) 
 
1097
      throw DudleyAdapterException("Error - Illegal domain of interpolation target");
 
1098
 
 
1099
   throw DudleyAdapterException("Error - Dudley does not allow interpolation across domains yet.");
 
1100
}
 
1101
 
 
1102
//
 
1103
// calculates the integral of a function defined of arg:
 
1104
//
 
1105
void MeshAdapter::setToIntegrals(vector<double>& integrals,const escript::Data& arg) const
 
1106
{
 
1107
   const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
 
1108
   if (argDomain!=*this) 
 
1109
      throw DudleyAdapterException("Error - Illegal domain of integration kernel");
 
1110
 
 
1111
   double blocktimer_start = blocktimer_time();
 
1112
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1113
   escriptDataC _temp;
 
1114
   escript::Data temp;
 
1115
   escriptDataC _arg=arg.getDataC();
 
1116
   switch(arg.getFunctionSpace().getTypeCode()) {
 
1117
   case(Nodes):
 
1118
   temp=escript::Data( arg, escript::function(*this) );
 
1119
   _temp=temp.getDataC();
 
1120
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
 
1121
   break;
 
1122
   case(ReducedNodes):
 
1123
   temp=escript::Data( arg, escript::function(*this) );
 
1124
   _temp=temp.getDataC();
 
1125
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
 
1126
   break;
 
1127
   case(Elements):
 
1128
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
 
1129
   break;
 
1130
   case(ReducedElements):
 
1131
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
 
1132
   break;
 
1133
   case(FaceElements):
 
1134
   Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
 
1135
   break;
 
1136
   case(ReducedFaceElements):
 
1137
   Dudley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);
 
1138
   break;
 
1139
   case(Points):
 
1140
   throw DudleyAdapterException("Error - Integral of data on points is not supported.");
 
1141
   break;
 
1142
   case(DegreesOfFreedom):
 
1143
   temp=escript::Data( arg, escript::function(*this) );
 
1144
   _temp=temp.getDataC();
 
1145
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
 
1146
   break;
 
1147
   case(ReducedDegreesOfFreedom):
 
1148
   temp=escript::Data( arg, escript::function(*this) );
 
1149
   _temp=temp.getDataC();
 
1150
   Dudley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
 
1151
   break;
 
1152
   default:
 
1153
      stringstream temp;
 
1154
      temp << "Error - Integrals: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
 
1155
      throw DudleyAdapterException(temp.str());
 
1156
      break;
 
1157
   }
 
1158
   checkDudleyError();
 
1159
   blocktimer_increment("integrate()", blocktimer_start);
 
1160
}
 
1161
 
 
1162
//
 
1163
// calculates the gradient of arg:
 
1164
//
 
1165
void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
 
1166
{
 
1167
   const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
 
1168
   if (argDomain!=*this)
 
1169
      throw DudleyAdapterException("Error - Illegal domain of gradient argument");
 
1170
   const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
 
1171
   if (gradDomain!=*this)
 
1172
      throw DudleyAdapterException("Error - Illegal domain of gradient");
 
1173
 
 
1174
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1175
   escriptDataC _grad=grad.getDataC();
 
1176
   escriptDataC nodeDataC;
 
1177
   escript::Data temp;
 
1178
   if (getMPISize()>1) {
 
1179
      if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
 
1180
         temp=escript::Data( arg,  continuousFunction(*this) );
 
1181
         nodeDataC = temp.getDataC();
 
1182
      } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
 
1183
         temp=escript::Data( arg,  reducedContinuousFunction(*this) );
 
1184
         nodeDataC = temp.getDataC();
 
1185
      } else {
 
1186
         nodeDataC = arg.getDataC();
 
1187
      }
 
1188
   } else {
 
1189
      nodeDataC = arg.getDataC();
 
1190
   }
 
1191
   switch(grad.getFunctionSpace().getTypeCode()) {
 
1192
   case(Nodes):
 
1193
   throw DudleyAdapterException("Error - Gradient at nodes is not supported.");
 
1194
   break;
 
1195
   case(ReducedNodes):
 
1196
   throw DudleyAdapterException("Error - Gradient at reduced nodes is not supported.");
 
1197
   break;
 
1198
   case(Elements):
 
1199
   Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
 
1200
   break;
 
1201
   case(ReducedElements):
 
1202
   Dudley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);
 
1203
   break;
 
1204
   case(FaceElements):
 
1205
   Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
 
1206
   break;
 
1207
   case(ReducedFaceElements):
 
1208
   Dudley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);
 
1209
   break;
 
1210
   case(Points):
 
1211
   throw DudleyAdapterException("Error - Gradient at points is not supported.");
 
1212
   break;
 
1213
   case(DegreesOfFreedom):
 
1214
   throw DudleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
 
1215
   break;
 
1216
   case(ReducedDegreesOfFreedom):
 
1217
   throw DudleyAdapterException("Error - Gradient at reduced degrees of freedom is not supported.");
 
1218
   break;
 
1219
   default:
 
1220
      stringstream temp;
 
1221
      temp << "Error - Gradient: Dudley does not know anything about function space type " << arg.getFunctionSpace().getTypeCode();
 
1222
      throw DudleyAdapterException(temp.str());
 
1223
      break;
 
1224
   }
 
1225
   checkDudleyError();
 
1226
}
 
1227
 
 
1228
//
 
1229
// returns the size of elements:
 
1230
//
 
1231
void MeshAdapter::setToSize(escript::Data& size) const
 
1232
{
 
1233
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1234
   escriptDataC tmp=size.getDataC();
 
1235
   switch(size.getFunctionSpace().getTypeCode()) {
 
1236
   case(Nodes):
 
1237
   throw DudleyAdapterException("Error - Size of nodes is not supported.");
 
1238
   break;
 
1239
   case(ReducedNodes):
 
1240
   throw DudleyAdapterException("Error - Size of reduced nodes is not supported.");
 
1241
   break;
 
1242
   case(Elements):
 
1243
   Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
 
1244
   break;
 
1245
   case(ReducedElements):
 
1246
   Dudley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);
 
1247
   break;
 
1248
   case(FaceElements):
 
1249
   Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
 
1250
   break;
 
1251
   case(ReducedFaceElements):
 
1252
   Dudley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);
 
1253
   break;
 
1254
   case(Points):
 
1255
   throw DudleyAdapterException("Error - Size of point elements is not supported.");
 
1256
   break;
 
1257
   case(DegreesOfFreedom):
 
1258
   throw DudleyAdapterException("Error - Size of degrees of freedom is not supported.");
 
1259
   break;
 
1260
   case(ReducedDegreesOfFreedom):
 
1261
   throw DudleyAdapterException("Error - Size of reduced degrees of freedom is not supported.");
 
1262
   break;
 
1263
   default:
 
1264
      stringstream temp;
 
1265
      temp << "Error - Element size: Dudley does not know anything about function space type " << size.getFunctionSpace().getTypeCode();
 
1266
      throw DudleyAdapterException(temp.str());
 
1267
      break;
 
1268
   }
 
1269
   checkDudleyError();
 
1270
}
 
1271
 
 
1272
//
 
1273
// sets the location of nodes
 
1274
//
 
1275
void MeshAdapter::setNewX(const escript::Data& new_x)
 
1276
{
 
1277
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1278
   escriptDataC tmp;
 
1279
   const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
 
1280
   if (newDomain!=*this) 
 
1281
      throw DudleyAdapterException("Error - Illegal domain of new point locations");
 
1282
   if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
 
1283
       tmp = new_x.getDataC();
 
1284
       Dudley_Mesh_setCoordinates(mesh,&tmp);
 
1285
   } else {
 
1286
       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
 
1287
       tmp = new_x_inter.getDataC();
 
1288
       Dudley_Mesh_setCoordinates(mesh,&tmp);
 
1289
   }
 
1290
   checkDudleyError();
 
1291
}
 
1292
 
 
1293
//
 
1294
// Helper for the save* methods. Extracts optional data variable names and the
 
1295
// corresponding pointers from python dictionary. Caller must free arrays.
 
1296
//
 
1297
void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
 
1298
{
 
1299
   numData = boost::python::extract<int>(arg.attr("__len__")());
 
1300
   /* win32 refactor */
 
1301
   names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
 
1302
   data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
 
1303
   dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
 
1304
 
 
1305
   boost::python::list keys=arg.keys();
 
1306
   for (int i=0; i<numData; ++i) {
 
1307
      string n=boost::python::extract<string>(keys[i]);
 
1308
      escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
 
1309
      if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this) 
 
1310
         throw DudleyAdapterException("Error: Data must be defined on same Domain");
 
1311
      data[i] = d.getDataC();
 
1312
      dataPtr[i] = &(data[i]);
 
1313
      names[i] = TMPMEMALLOC(n.length()+1, char);
 
1314
      strcpy(names[i], n.c_str());
 
1315
   }
 
1316
}
 
1317
 
 
1318
//
 
1319
// saves mesh and optionally data arrays in openDX format
 
1320
//
 
1321
void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const
 
1322
{
 
1323
   int num_data;
 
1324
   char **names;
 
1325
   escriptDataC *data;
 
1326
   escriptDataC **ptr_data;
 
1327
 
 
1328
   extractArgsFromDict(arg, num_data, names, data, ptr_data);
 
1329
   Dudley_Mesh_saveDX(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data);
 
1330
   checkDudleyError();
 
1331
 
 
1332
   /* win32 refactor */
 
1333
   TMPMEMFREE(data);
 
1334
   TMPMEMFREE(ptr_data);
 
1335
   for(int i=0; i<num_data; i++)
 
1336
   {
 
1337
      TMPMEMFREE(names[i]);
 
1338
   }
 
1339
   TMPMEMFREE(names);
 
1340
 
 
1341
   return;
 
1342
}
 
1343
 
 
1344
//
 
1345
// saves mesh and optionally data arrays in VTK format
 
1346
//
 
1347
void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
 
1348
{
 
1349
   int num_data;
 
1350
   char **names;
 
1351
   escriptDataC *data;
 
1352
   escriptDataC **ptr_data;
 
1353
 
 
1354
   extractArgsFromDict(arg, num_data, names, data, ptr_data);
 
1355
   Dudley_Mesh_saveVTK(filename.c_str(), m_dudleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
 
1356
   checkDudleyError();
 
1357
 
 
1358
   /* win32 refactor */
 
1359
   TMPMEMFREE(data);
 
1360
   TMPMEMFREE(ptr_data);
 
1361
   for(int i=0; i<num_data; i++)
 
1362
   {
 
1363
      TMPMEMFREE(names[i]);
 
1364
   }
 
1365
   TMPMEMFREE(names);
 
1366
}
 
1367
 
 
1368
bool MeshAdapter::ownSample(int fs_code, index_t id) const
 
1369
{
 
1370
#ifdef ESYS_MPI
 
1371
    index_t myFirstNode=0, myLastNode=0, k=0;
 
1372
    index_t* globalNodeIndex=0;
 
1373
    Dudley_Mesh* mesh_p=m_dudleyMesh.get();
 
1374
    if (fs_code == DUDLEY_REDUCED_NODES) 
 
1375
    {
 
1376
        myFirstNode = Dudley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
 
1377
        myLastNode = Dudley_NodeFile_getLastReducedNode(mesh_p->Nodes);
 
1378
        globalNodeIndex = Dudley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
 
1379
    }
 
1380
    else
 
1381
    {
 
1382
        myFirstNode = Dudley_NodeFile_getFirstNode(mesh_p->Nodes);
 
1383
        myLastNode = Dudley_NodeFile_getLastNode(mesh_p->Nodes);
 
1384
        globalNodeIndex = Dudley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
 
1385
    }
 
1386
    k=globalNodeIndex[id];
 
1387
    return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
 
1388
#endif
 
1389
    return true;
 
1390
}
 
1391
 
 
1392
 
 
1393
 
 
1394
//
 
1395
// creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
 
1396
//
 
1397
ASM_ptr MeshAdapter::newSystemMatrix(
 
1398
                                                 const int row_blocksize,
 
1399
                                                 const escript::FunctionSpace& row_functionspace,
 
1400
                                                 const int column_blocksize,
 
1401
                                                 const escript::FunctionSpace& column_functionspace,
 
1402
                                                 const int type) const
 
1403
{
 
1404
   int reduceRowOrder=0;
 
1405
   int reduceColOrder=0;
 
1406
   // is the domain right?
 
1407
   const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
 
1408
   if (row_domain!=*this) 
 
1409
      throw DudleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
 
1410
   const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
 
1411
   if (col_domain!=*this) 
 
1412
      throw DudleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
 
1413
   // is the function space type right 
 
1414
   if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
 
1415
      reduceRowOrder=0;
 
1416
   } else if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
 
1417
      reduceRowOrder=1;
 
1418
   } else {
 
1419
      throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
 
1420
   }
 
1421
   if (column_functionspace.getTypeCode()==DegreesOfFreedom) {
 
1422
      reduceColOrder=0;
 
1423
   } else if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
 
1424
      reduceColOrder=1;
 
1425
   } else {
 
1426
      throw DudleyAdapterException("Error - illegal function space type for system matrix columns.");
 
1427
   }
 
1428
   // generate matrix:
 
1429
 
 
1430
   Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder);
 
1431
   checkDudleyError();
 
1432
   Paso_SystemMatrix* fsystemMatrix;
 
1433
   int trilinos = 0;
 
1434
   if (trilinos) {
 
1435
#ifdef TRILINOS
 
1436
      /* Allocation Epetra_VrbMatrix here */
 
1437
#endif
 
1438
   }
 
1439
   else {
 
1440
      fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
 
1441
   }
 
1442
   checkPasoError();
 
1443
   Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
 
1444
   SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace, column_blocksize,column_functionspace);
 
1445
   return ASM_ptr(sma);
 
1446
}
 
1447
 
 
1448
//
 
1449
// creates a TransportProblemAdapter
 
1450
//
 
1451
ATP_ptr MeshAdapter::newTransportProblem(
 
1452
                                                         const bool useBackwardEuler,
 
1453
                                                         const int blocksize,
 
1454
                                                         const escript::FunctionSpace& functionspace,
 
1455
                                                         const int type) const
 
1456
{
 
1457
   int reduceOrder=0;
 
1458
   // is the domain right?
 
1459
   const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
 
1460
   if (domain!=*this) 
 
1461
      throw DudleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
 
1462
   // is the function space type right 
 
1463
   if (functionspace.getTypeCode()==DegreesOfFreedom) {
 
1464
      reduceOrder=0;
 
1465
   } else if (functionspace.getTypeCode()==ReducedDegreesOfFreedom) {
 
1466
      reduceOrder=1;
 
1467
   } else {
 
1468
      throw DudleyAdapterException("Error - illegal function space type for system matrix rows.");
 
1469
   }
 
1470
   // generate matrix:
 
1471
 
 
1472
   Paso_SystemMatrixPattern* fsystemMatrixPattern=Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder);
 
1473
   checkDudleyError();
 
1474
   Paso_TransportProblem* transportProblem;
 
1475
   transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
 
1476
   checkPasoError();
 
1477
   Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
 
1478
   AbstractTransportProblem* atp=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
 
1479
   return ATP_ptr(atp);
 
1480
}
 
1481
 
 
1482
//
 
1483
// vtkObject MeshAdapter::createVtkObject() const
 
1484
// TODO:
 
1485
//
 
1486
 
 
1487
//
 
1488
// returns true if data at the atom_type is considered as being cell centered:
 
1489
bool MeshAdapter::isCellOriented(int functionSpaceCode) const
 
1490
{
 
1491
   switch(functionSpaceCode) {
 
1492
   case(Nodes):
 
1493
   case(DegreesOfFreedom):
 
1494
   case(ReducedDegreesOfFreedom):
 
1495
   return false;
 
1496
   break;
 
1497
   case(Elements):
 
1498
   case(FaceElements):
 
1499
   case(Points):
 
1500
   case(ReducedElements):
 
1501
   case(ReducedFaceElements):
 
1502
   return true;
 
1503
   break;
 
1504
   default:
 
1505
      stringstream temp;
 
1506
      temp << "Error - Cell: Dudley does not know anything about function space type " << functionSpaceCode;
 
1507
      throw DudleyAdapterException(temp.str());
 
1508
      break;
 
1509
   }
 
1510
   checkDudleyError();
 
1511
   return false;
 
1512
}
 
1513
 
 
1514
bool
 
1515
MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
 
1516
{
 
1517
   /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
 
1518
        class 1: DOF <-> Nodes
 
1519
        class 2: ReducedDOF <-> ReducedNodes
 
1520
        class 3: Points
 
1521
        class 4: Elements
 
1522
        class 5: ReducedElements
 
1523
        class 6: FaceElements
 
1524
        class 7: ReducedFaceElements
 
1525
        class 8: ContactElementZero <-> ContactElementOne
 
1526
        class 9: ReducedContactElementZero <-> ReducedContactElementOne
 
1527
 
 
1528
   There is also a set of lines. Interpolation is possible down a line but not between lines.
 
1529
   class 1 and 2 belong to all lines so aren't considered.
 
1530
        line 0: class 3
 
1531
        line 1: class 4,5
 
1532
        line 2: class 6,7
 
1533
        line 3: class 8,9
 
1534
 
 
1535
   For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
 
1536
   eg hasnodes is true if we have at least one instance of Nodes.
 
1537
   */
 
1538
    if (fs.empty())
 
1539
    {
 
1540
        return false;
 
1541
    }
 
1542
    vector<int> hasclass(10);
 
1543
    vector<int> hasline(4);     
 
1544
    bool hasnodes=false;
 
1545
    bool hasrednodes=false;
 
1546
    for (int i=0;i<fs.size();++i)
 
1547
    {
 
1548
        switch(fs[i])
 
1549
        {
 
1550
        case(Nodes):   hasnodes=true;   // no break is deliberate
 
1551
        case(DegreesOfFreedom):
 
1552
                hasclass[1]=1;
 
1553
                break;
 
1554
        case(ReducedNodes):    hasrednodes=true;        // no break is deliberate
 
1555
        case(ReducedDegreesOfFreedom):
 
1556
                hasclass[2]=1;
 
1557
                break;
 
1558
        case(Points):
 
1559
                hasline[0]=1;
 
1560
                hasclass[3]=1;
 
1561
                break;
 
1562
        case(Elements):
 
1563
                hasclass[4]=1;
 
1564
                hasline[1]=1;
 
1565
                break;
 
1566
        case(ReducedElements):
 
1567
                hasclass[5]=1;
 
1568
                hasline[1]=1;
 
1569
                break;
 
1570
        case(FaceElements):
 
1571
                hasclass[6]=1;
 
1572
                hasline[2]=1;
 
1573
                break;
 
1574
        case(ReducedFaceElements):
 
1575
                hasclass[7]=1;
 
1576
                hasline[2]=1;
 
1577
                break;
 
1578
        default:
 
1579
                return false;
 
1580
        }
 
1581
    }
 
1582
    int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
 
1583
    // fail if we have more than one leaf group
 
1584
 
 
1585
    if (totlines>1)
 
1586
    {
 
1587
        return false;   // there are at least two branches we can't interpolate between
 
1588
    }
 
1589
    else if (totlines==1)
 
1590
    {
 
1591
        if (hasline[0]==1)              // we have points
 
1592
        {
 
1593
            resultcode=Points;
 
1594
        }
 
1595
        else if (hasline[1]==1)
 
1596
        {
 
1597
            if (hasclass[5]==1)
 
1598
            {
 
1599
                resultcode=ReducedElements;
 
1600
            }
 
1601
            else
 
1602
            {
 
1603
                resultcode=Elements;
 
1604
            }
 
1605
        }
 
1606
        else if (hasline[2]==1)
 
1607
        {
 
1608
            if (hasclass[7]==1)
 
1609
            {
 
1610
                resultcode=ReducedFaceElements;
 
1611
            }
 
1612
            else
 
1613
            {
 
1614
                resultcode=FaceElements;
 
1615
            }
 
1616
        }
 
1617
        else    // so we must be in line3
 
1618
        {
 
1619
 
 
1620
            throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
 
1621
 
 
1622
        }
 
1623
    }
 
1624
    else        // totlines==0
 
1625
    {
 
1626
        if (hasclass[2]==1)
 
1627
        {
 
1628
                // something from class 2
 
1629
                resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
 
1630
        }
 
1631
        else
 
1632
        {       // something from class 1
 
1633
                resultcode=(hasnodes?Nodes:DegreesOfFreedom);
 
1634
        }
 
1635
    }
 
1636
    return true;
 
1637
}
 
1638
 
 
1639
bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
 
1640
{
 
1641
   switch(functionSpaceType_source) {
 
1642
   case(Nodes):
 
1643
        switch(functionSpaceType_target) {
 
1644
        case(Nodes):
 
1645
        case(ReducedNodes):
 
1646
        case(ReducedDegreesOfFreedom):
 
1647
        case(DegreesOfFreedom):
 
1648
        case(Elements):
 
1649
        case(ReducedElements):
 
1650
        case(FaceElements):
 
1651
        case(ReducedFaceElements):
 
1652
        case(Points):
 
1653
        return true;
 
1654
        default:
 
1655
              stringstream temp;
 
1656
              temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
 
1657
              throw DudleyAdapterException(temp.str());
 
1658
   }
 
1659
   break;
 
1660
   case(ReducedNodes):
 
1661
        switch(functionSpaceType_target) {
 
1662
        case(ReducedNodes):
 
1663
        case(ReducedDegreesOfFreedom):
 
1664
        case(Elements):
 
1665
        case(ReducedElements):
 
1666
        case(FaceElements):
 
1667
        case(ReducedFaceElements):
 
1668
        case(Points):
 
1669
        return true;
 
1670
        case(Nodes):
 
1671
        case(DegreesOfFreedom):
 
1672
        return false;
 
1673
        default:
 
1674
                stringstream temp;
 
1675
                temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
 
1676
                throw DudleyAdapterException(temp.str());
 
1677
   }
 
1678
   break;
 
1679
   case(Elements):
 
1680
        if (functionSpaceType_target==Elements) {
 
1681
          return true;
 
1682
        } else if (functionSpaceType_target==ReducedElements) {
 
1683
          return true;
 
1684
        } else {
 
1685
          return false;
 
1686
        }
 
1687
   case(ReducedElements):
 
1688
        if (functionSpaceType_target==ReducedElements) {
 
1689
          return true;
 
1690
        } else {
 
1691
          return false;
 
1692
        }
 
1693
   case(FaceElements):
 
1694
        if (functionSpaceType_target==FaceElements) {
 
1695
                return true;
 
1696
        } else if (functionSpaceType_target==ReducedFaceElements) {
 
1697
                return true;
 
1698
        } else {
 
1699
                return false;
 
1700
        }
 
1701
   case(ReducedFaceElements):
 
1702
        if (functionSpaceType_target==ReducedFaceElements) {
 
1703
                return true;
 
1704
        } else {
 
1705
                return false;
 
1706
        }
 
1707
   case(Points):
 
1708
        if (functionSpaceType_target==Points) {
 
1709
                return true;
 
1710
        } else {
 
1711
                return false;
 
1712
        }
 
1713
   case(DegreesOfFreedom):
 
1714
        switch(functionSpaceType_target) {
 
1715
        case(ReducedDegreesOfFreedom):
 
1716
        case(DegreesOfFreedom):
 
1717
        case(Nodes):
 
1718
        case(ReducedNodes):
 
1719
        case(Elements):
 
1720
        case(ReducedElements):
 
1721
        case(Points):
 
1722
        case(FaceElements):
 
1723
        case(ReducedFaceElements):
 
1724
        return true;
 
1725
        default:
 
1726
                stringstream temp;
 
1727
                temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
 
1728
                throw DudleyAdapterException(temp.str());
 
1729
        }
 
1730
        break;
 
1731
   case(ReducedDegreesOfFreedom):
 
1732
   switch(functionSpaceType_target) {
 
1733
        case(ReducedDegreesOfFreedom):
 
1734
        case(ReducedNodes):
 
1735
        case(Elements):
 
1736
        case(ReducedElements):
 
1737
        case(FaceElements):
 
1738
        case(ReducedFaceElements):
 
1739
        case(Points):
 
1740
        return true;
 
1741
        case(Nodes):
 
1742
        case(DegreesOfFreedom):
 
1743
        return false;
 
1744
        default:
 
1745
                stringstream temp;
 
1746
                temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
 
1747
                throw DudleyAdapterException(temp.str());
 
1748
        }
 
1749
        break;
 
1750
   default:
 
1751
      stringstream temp;
 
1752
      temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
 
1753
      throw DudleyAdapterException(temp.str());
 
1754
      break;
 
1755
   }
 
1756
   checkDudleyError();
 
1757
   return false;
 
1758
}
 
1759
 
 
1760
bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
 
1761
{
 
1762
   return false;
 
1763
}
 
1764
 
 
1765
bool MeshAdapter::operator==(const AbstractDomain& other) const
 
1766
{
 
1767
   const MeshAdapter* temp=dynamic_cast<const MeshAdapter*>(&other);
 
1768
   if (temp!=0) {
 
1769
      return (m_dudleyMesh==temp->m_dudleyMesh);
 
1770
   } else {
 
1771
      return false;
 
1772
   }
 
1773
}
 
1774
 
 
1775
bool MeshAdapter::operator!=(const AbstractDomain& other) const
 
1776
{
 
1777
   return !(operator==(other));
 
1778
}
 
1779
 
 
1780
int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
 
1781
{
 
1782
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1783
   int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
 
1784
   checkPasoError();
 
1785
   return out;
 
1786
}
 
1787
int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
 
1788
{
 
1789
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1790
   int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
 
1791
   checkPasoError();
 
1792
   return out;
 
1793
}
 
1794
 
 
1795
escript::Data MeshAdapter::getX() const
 
1796
{
 
1797
   return continuousFunction(*this).getX();
 
1798
}
 
1799
 
 
1800
escript::Data MeshAdapter::getNormal() const
 
1801
{
 
1802
   return functionOnBoundary(*this).getNormal();
 
1803
}
 
1804
 
 
1805
escript::Data MeshAdapter::getSize() const
 
1806
{
 
1807
   return escript::function(*this).getSize();
 
1808
}
 
1809
 
 
1810
const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
 
1811
{
 
1812
   int *out = NULL;
 
1813
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1814
   switch (functionSpaceType) {
 
1815
   case(Nodes):
 
1816
   out=mesh->Nodes->Id;
 
1817
   break;
 
1818
   case(ReducedNodes):
 
1819
   out=mesh->Nodes->reducedNodesId;
 
1820
   break;
 
1821
   case(Elements):
 
1822
   out=mesh->Elements->Id;
 
1823
   break;
 
1824
   case(ReducedElements):
 
1825
   out=mesh->Elements->Id;
 
1826
   break;
 
1827
   case(FaceElements):
 
1828
   out=mesh->FaceElements->Id;
 
1829
   break;
 
1830
   case(ReducedFaceElements):
 
1831
   out=mesh->FaceElements->Id;
 
1832
   break;
 
1833
   case(Points):
 
1834
   out=mesh->Points->Id;
 
1835
   break;
 
1836
   case(DegreesOfFreedom):
 
1837
   out=mesh->Nodes->degreesOfFreedomId;
 
1838
   break;
 
1839
   case(ReducedDegreesOfFreedom):
 
1840
   out=mesh->Nodes->reducedDegreesOfFreedomId;
 
1841
   break;
 
1842
   default:
 
1843
      stringstream temp;
 
1844
      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
 
1845
      throw DudleyAdapterException(temp.str());
 
1846
      break;
 
1847
   }
 
1848
   return out;
 
1849
}
 
1850
int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
 
1851
{
 
1852
   int out=0;
 
1853
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1854
   switch (functionSpaceType) {
 
1855
   case(Nodes):
 
1856
   out=mesh->Nodes->Tag[sampleNo];
 
1857
   break;
 
1858
   case(ReducedNodes):
 
1859
   throw DudleyAdapterException(" Error - ReducedNodes does not support tags.");
 
1860
   break;
 
1861
   case(Elements):
 
1862
   out=mesh->Elements->Tag[sampleNo];
 
1863
   break;
 
1864
   case(ReducedElements):
 
1865
   out=mesh->Elements->Tag[sampleNo];
 
1866
   break;
 
1867
   case(FaceElements):
 
1868
   out=mesh->FaceElements->Tag[sampleNo];
 
1869
   break;
 
1870
   case(ReducedFaceElements):
 
1871
   out=mesh->FaceElements->Tag[sampleNo];
 
1872
   break;
 
1873
   case(Points):
 
1874
   out=mesh->Points->Tag[sampleNo];
 
1875
   break;
 
1876
   case(DegreesOfFreedom):
 
1877
   throw DudleyAdapterException(" Error - DegreesOfFreedom does not support tags.");
 
1878
   break;
 
1879
   case(ReducedDegreesOfFreedom):
 
1880
   throw DudleyAdapterException(" Error - ReducedDegreesOfFreedom does not support tags.");
 
1881
   break;
 
1882
   default:
 
1883
      stringstream temp;
 
1884
      temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
 
1885
      throw DudleyAdapterException(temp.str());
 
1886
      break;
 
1887
   }
 
1888
   return out;
 
1889
}
 
1890
 
 
1891
 
 
1892
void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
 
1893
{
 
1894
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1895
   escriptDataC tmp=mask.getDataC();
 
1896
   switch(functionSpaceType) {
 
1897
   case(Nodes):
 
1898
   Dudley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);
 
1899
   break;
 
1900
   case(ReducedNodes):
 
1901
   throw DudleyAdapterException("Error - ReducedNodes does not support tags");
 
1902
   break;
 
1903
   case(DegreesOfFreedom):
 
1904
   throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
 
1905
   break;
 
1906
   case(ReducedDegreesOfFreedom):
 
1907
   throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
 
1908
   break;
 
1909
   case(Elements):
 
1910
   Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
 
1911
   break;
 
1912
   case(ReducedElements):
 
1913
   Dudley_ElementFile_setTags(mesh->Elements,newTag,&tmp);
 
1914
   break;
 
1915
   case(FaceElements):
 
1916
   Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
 
1917
   break;
 
1918
   case(ReducedFaceElements):
 
1919
   Dudley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);
 
1920
   break;
 
1921
   case(Points):
 
1922
   Dudley_ElementFile_setTags(mesh->Points,newTag,&tmp);
 
1923
   break;
 
1924
   default:
 
1925
      stringstream temp;
 
1926
      temp << "Error - Dudley does not know anything about function space type " << functionSpaceType;
 
1927
      throw DudleyAdapterException(temp.str());
 
1928
   }
 
1929
   checkDudleyError();
 
1930
   return;
 
1931
}
 
1932
 
 
1933
void MeshAdapter::setTagMap(const string& name,  int tag)
 
1934
{
 
1935
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1936
   Dudley_Mesh_addTagMap(mesh, name.c_str(),tag);
 
1937
   checkPasoError();
 
1938
   // throwStandardException("MeshAdapter::set TagMap is not implemented.");
 
1939
}
 
1940
 
 
1941
int MeshAdapter::getTag(const string& name) const
 
1942
{
 
1943
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1944
   int tag=0;
 
1945
   tag=Dudley_Mesh_getTag(mesh, name.c_str());
 
1946
   checkPasoError();
 
1947
   // throwStandardException("MeshAdapter::getTag is not implemented.");
 
1948
   return tag;
 
1949
}
 
1950
 
 
1951
bool MeshAdapter::isValidTagName(const string& name) const
 
1952
{
 
1953
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1954
   return Dudley_Mesh_isValidTagName(mesh,name.c_str());
 
1955
}
 
1956
 
 
1957
string MeshAdapter::showTagNames() const
 
1958
{
 
1959
   stringstream temp;
 
1960
   Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1961
   Dudley_TagMap* tag_map=mesh->TagMap;
 
1962
   while (tag_map) {
 
1963
      temp << tag_map->name;
 
1964
      tag_map=tag_map->next;
 
1965
      if (tag_map) temp << ", ";
 
1966
   }
 
1967
   return temp.str();
 
1968
}
 
1969
 
 
1970
int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
 
1971
{
 
1972
  Dudley_Mesh* mesh=m_dudleyMesh.get();
 
1973
  dim_t numTags=0;
 
1974
  switch(functionSpaceCode) {
 
1975
   case(Nodes):
 
1976
          numTags=mesh->Nodes->numTagsInUse;
 
1977
          break;
 
1978
   case(ReducedNodes):
 
1979
          throw DudleyAdapterException("Error - ReducedNodes does not support tags");
 
1980
          break;
 
1981
   case(DegreesOfFreedom):
 
1982
          throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
 
1983
          break;
 
1984
   case(ReducedDegreesOfFreedom):
 
1985
          throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
 
1986
          break;
 
1987
   case(Elements):
 
1988
   case(ReducedElements):
 
1989
          numTags=mesh->Elements->numTagsInUse;
 
1990
          break;
 
1991
   case(FaceElements):
 
1992
   case(ReducedFaceElements):
 
1993
          numTags=mesh->FaceElements->numTagsInUse;
 
1994
          break;
 
1995
   case(Points):
 
1996
          numTags=mesh->Points->numTagsInUse;
 
1997
          break;
 
1998
   default:
 
1999
      stringstream temp;
 
2000
      temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
 
2001
      throw DudleyAdapterException(temp.str());
 
2002
  }
 
2003
  return numTags;
 
2004
}
 
2005
 
 
2006
const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
 
2007
{
 
2008
  Dudley_Mesh* mesh=m_dudleyMesh.get();
 
2009
  index_t* tags=NULL;
 
2010
  switch(functionSpaceCode) {
 
2011
   case(Nodes):
 
2012
          tags=mesh->Nodes->tagsInUse;
 
2013
          break;
 
2014
   case(ReducedNodes):
 
2015
          throw DudleyAdapterException("Error - ReducedNodes does not support tags");
 
2016
          break;
 
2017
   case(DegreesOfFreedom):
 
2018
          throw DudleyAdapterException("Error - DegreesOfFreedom does not support tags");
 
2019
          break;
 
2020
   case(ReducedDegreesOfFreedom):
 
2021
          throw DudleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
 
2022
          break;
 
2023
   case(Elements):
 
2024
   case(ReducedElements):
 
2025
          tags=mesh->Elements->tagsInUse;
 
2026
          break;
 
2027
   case(FaceElements):
 
2028
   case(ReducedFaceElements):
 
2029
          tags=mesh->FaceElements->tagsInUse;
 
2030
          break;
 
2031
   case(Points):
 
2032
          tags=mesh->Points->tagsInUse;
 
2033
          break;
 
2034
   default:
 
2035
      stringstream temp;
 
2036
      temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
 
2037
      throw DudleyAdapterException(temp.str());
 
2038
  }
 
2039
  return tags;
 
2040
}
 
2041
 
 
2042
 
 
2043
bool MeshAdapter::canTag(int functionSpaceCode) const
 
2044
{
 
2045
  switch(functionSpaceCode) {
 
2046
   case(Nodes):
 
2047
   case(Elements):
 
2048
   case(ReducedElements):
 
2049
   case(FaceElements):
 
2050
   case(ReducedFaceElements):
 
2051
   case(Points):
 
2052
          return true;
 
2053
   case(ReducedNodes):
 
2054
   case(DegreesOfFreedom):
 
2055
   case(ReducedDegreesOfFreedom):
 
2056
          return false;
 
2057
   default:
 
2058
        return false;
 
2059
  }
 
2060
}
 
2061
 
 
2062
AbstractDomain::StatusType MeshAdapter::getStatus() const
 
2063
{
 
2064
  Dudley_Mesh* mesh=m_dudleyMesh.get();
 
2065
  return Dudley_Mesh_getStatus(mesh);
 
2066
}
 
2067
 
 
2068
int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
 
2069
{
 
2070
   
 
2071
  Dudley_Mesh* mesh=m_dudleyMesh.get();
 
2072
  int order =-1;
 
2073
  switch(functionSpaceCode) {
 
2074
   case(Nodes):
 
2075
   case(DegreesOfFreedom):
 
2076
          order=mesh->approximationOrder;
 
2077
          break;
 
2078
   case(ReducedNodes):
 
2079
   case(ReducedDegreesOfFreedom):
 
2080
          order=mesh->reducedApproximationOrder;
 
2081
          break;
 
2082
   case(Elements):
 
2083
   case(FaceElements):
 
2084
   case(Points):
 
2085
          order=mesh->integrationOrder;
 
2086
          break;
 
2087
   case(ReducedElements):
 
2088
   case(ReducedFaceElements):
 
2089
          order=mesh->reducedIntegrationOrder;
 
2090
          break;
 
2091
   default:
 
2092
      stringstream temp;
 
2093
      temp << "Error - Dudley does not know anything about function space type " << functionSpaceCode;
 
2094
      throw DudleyAdapterException(temp.str());
 
2095
  }
 
2096
  return order;
 
2097
}
 
2098
 
 
2099
 
 
2100
bool MeshAdapter::supportsContactElements() const
 
2101
{
 
2102
    return false;
 
2103
}
 
2104
 
 
2105
}  // end of namespace