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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

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

Legend:
Removed from v.532  
changed lines
  Added in v.3490

  ViewVC Help
Powered by ViewVC 1.1.26