/[escript]/branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/domexper/dudley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.82  
changed lines
  Added in v.3090

  ViewVC Help
Powered by ViewVC 1.1.26