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

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

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

revision 4408 by caltinay, Tue May 14 06:58:43 2013 UTC revision 4800 by caltinay, Wed Mar 26 01:50:04 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
# Line 19  Line 20 
20  #include "escript/Data.h"  #include "escript/Data.h"
21  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23    #include "esysUtils/EsysRandom.h"
24    
25  #include <boost/python/import.hpp>  #include <boost/python/import.hpp>
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 46  const int MeshAdapter::ReducedContactEle Line 48  const int MeshAdapter::ReducedContactEle
48  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
49  const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;  const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Mesh* finleyMesh)
52  {  {
53      setFunctionSpaceTypeNames();      setFunctionSpaceTypeNames();
54      // need to use a null_deleter as Finley_Mesh_free deletes the pointer      // need to use a null_deleter as Finley_Mesh_free deletes the pointer
# Line 66  MeshAdapter::~MeshAdapter() Line 68  MeshAdapter::~MeshAdapter()
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      if (m_finleyMesh.unique()) {      if (m_finleyMesh.unique()) {
71          Finley_Mesh_free(m_finleyMesh.get());          delete m_finleyMesh.get();
72      }      }
73  }  }
74    
# Line 106  MeshAdapter::getMPIComm() const Line 108  MeshAdapter::getMPIComm() const
108  #endif  #endif
109  }  }
110    
111  Finley_Mesh* MeshAdapter::getFinley_Mesh() const  Mesh* MeshAdapter::getFinley_Mesh() const
112  {  {
113      return m_finleyMesh.get();      return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118      Finley_Mesh_write(m_finleyMesh.get(), fileName.c_str());      m_finleyMesh.get()->write(fileName);
119      checkFinleyError();      checkFinleyError();
120  }  }
121    
122  void MeshAdapter::Print_Mesh_Info(const bool full) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
123  {  {
124      Finley_PrintMesh_Info(m_finleyMesh.get(), full);      m_finleyMesh.get()->printInfo(full);
125  }  }
126    
127  void MeshAdapter::dump(const string& fileName) const  void MeshAdapter::dump(const string& fileName) const
# Line 128  void MeshAdapter::dump(const string& fil Line 130  void MeshAdapter::dump(const string& fil
130      const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};      const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
131      NcVar *ids;      NcVar *ids;
132      int *int_ptr;      int *int_ptr;
133      Finley_Mesh *mesh = m_finleyMesh.get();      Mesh *mesh = m_finleyMesh.get();
     Finley_TagMap* tag_map;  
134      int num_Tags = 0;      int num_Tags = 0;
135      int mpi_size                         = mesh->MPIInfo->size;      int mpi_size                         = mesh->MPIInfo->size;
136      int mpi_rank                         = mesh->MPIInfo->rank;      int mpi_rank                         = mesh->MPIInfo->rank;
# Line 153  void MeshAdapter::dump(const string& fil Line 154  void MeshAdapter::dump(const string& fil
154      }      }
155  #endif  #endif
156    
157      char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),      const std::string newFileName(esysUtils::appendRankToFileName(fileName,
158                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank));
159    
160      // Figure out how much storage is required for tags      // Figure out how much storage is required for tags
161      tag_map = mesh->TagMap;      num_Tags = mesh->tagMap.size();
     num_Tags = 0;  
     while (tag_map) {  
         num_Tags++;  
         tag_map=tag_map->next;  
     }  
162    
163      // NetCDF error handler      // NetCDF error handler
164      NcError err(NcError::verbose_nonfatal);      NcError err(NcError::verbose_nonfatal);
165      // Create the file.      // Create the file.
166      NcFile dataFile(newFileName, NcFile::Replace);      NcFile dataFile(newFileName.c_str(), NcFile::Replace);
167      string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");      string msgPrefix("Error in MeshAdapter::dump: NetCDF operation failed - ");
168      // check if writing was successful      // check if writing was successful
169      if (!dataFile.is_valid())      if (!dataFile.is_valid())
# Line 207  void MeshAdapter::dump(const string& fil Line 203  void MeshAdapter::dump(const string& fil
203              throw FinleyAdapterException(msgPrefix+"add_dim(dim_Tags)");              throw FinleyAdapterException(msgPrefix+"add_dim(dim_Tags)");
204    
205      // Attributes: MPI size, MPI rank, Name, order, reduced_order      // Attributes: MPI size, MPI rank, Name, order, reduced_order
206      if (!dataFile.add_att("mpi_size", mpi_size) )      if (!dataFile.add_att("mpi_size", mpi_size))
207          throw FinleyAdapterException(msgPrefix+"add_att(mpi_size)");          throw FinleyAdapterException(msgPrefix+"add_att(mpi_size)");
208      if (!dataFile.add_att("mpi_rank", mpi_rank) )      if (!dataFile.add_att("mpi_rank", mpi_rank))
209          throw FinleyAdapterException(msgPrefix+"add_att(mpi_rank)");          throw FinleyAdapterException(msgPrefix+"add_att(mpi_rank)");
210      if (!dataFile.add_att("Name",mesh->Name) )      if (!dataFile.add_att("Name",mesh->m_name.c_str()))
211          throw FinleyAdapterException(msgPrefix+"add_att(Name)");          throw FinleyAdapterException(msgPrefix+"add_att(Name)");
212      if (!dataFile.add_att("numDim",numDim) )      if (!dataFile.add_att("numDim",numDim))
213          throw FinleyAdapterException(msgPrefix+"add_att(order)");          throw FinleyAdapterException(msgPrefix+"add_att(order)");
214      if (!dataFile.add_att("order",mesh->integrationOrder) )      if (!dataFile.add_att("order",mesh->integrationOrder))
215          throw FinleyAdapterException(msgPrefix+"add_att(order)");          throw FinleyAdapterException(msgPrefix+"add_att(order)");
216      if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )      if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder))
217          throw FinleyAdapterException(msgPrefix+"add_att(reduced_order)");          throw FinleyAdapterException(msgPrefix+"add_att(reduced_order)");
218      if (!dataFile.add_att("numNodes",numNodes) )      if (!dataFile.add_att("numNodes",numNodes))
219          throw FinleyAdapterException(msgPrefix+"add_att(numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(numNodes)");
220      if (!dataFile.add_att("num_Elements",num_Elements) )      if (!dataFile.add_att("num_Elements",num_Elements))
221          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements)");
222      if (!dataFile.add_att("num_FaceElements",num_FaceElements) )      if (!dataFile.add_att("num_FaceElements",num_FaceElements))
223          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements)");
224      if (!dataFile.add_att("num_ContactElements",num_ContactElements) )      if (!dataFile.add_att("num_ContactElements",num_ContactElements))
225          throw FinleyAdapterException(msgPrefix+"add_att(num_ContactElements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_ContactElements)");
226      if (!dataFile.add_att("num_Points",num_Points) )      if (!dataFile.add_att("num_Points",num_Points))
227          throw FinleyAdapterException(msgPrefix+"add_att(num_Points)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Points)");
228      if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )      if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes))
229          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements_numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements_numNodes)");
230      if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )      if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
231          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements_numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements_numNodes)");
# Line 470  void MeshAdapter::dump(const string& fil Line 466  void MeshAdapter::dump(const string& fil
466          vector<int> Tags_keys;          vector<int> Tags_keys;
467    
468          // Copy tag data into temp arrays          // Copy tag data into temp arrays
469          tag_map = mesh->TagMap;          TagMap::const_iterator it;
470          while (tag_map) {          for (it=mesh->tagMap.begin(); it!=mesh->tagMap.end(); it++) {
471              Tags_keys.push_back(tag_map->tag_key);              Tags_keys.push_back(it->second);
             tag_map=tag_map->next;  
472          }          }
473    
474          // Tags_keys          // Tags_keys
# Line 487  void MeshAdapter::dump(const string& fil Line 482  void MeshAdapter::dump(const string& fil
482          // This is an array of strings, it should be stored as an array but          // This is an array of strings, it should be stored as an array but
483          // instead I have hacked in one attribute per string because the NetCDF          // instead I have hacked in one attribute per string because the NetCDF
484          // manual doesn't tell how to do an array of strings          // manual doesn't tell how to do an array of strings
         tag_map = mesh->TagMap;  
485          int i = 0;          int i = 0;
486          while (tag_map) {          for (it=mesh->tagMap.begin(); it!=mesh->tagMap.end(); it++, i++) {
487              stringstream tagnamestream;              stringstream tagnamestream;
488              tagnamestream << "Tags_name_" << i;              tagnamestream << "Tags_name_" << i;
489              const string tagname = tagnamestream.str();              const string tagname = tagnamestream.str();
490              if (!dataFile.add_att(tagname.c_str(), tag_map->name) )              if (!dataFile.add_att(tagname.c_str(), it->first.c_str()))
491                  throw FinleyAdapterException(msgPrefix+"add_att(Tags_names_XX)");                  throw FinleyAdapterException(msgPrefix+"add_att(Tags_names_X)");
             tag_map=tag_map->next;  
             i++;  
492          }          }
493      }      }
494    
# Line 510  void MeshAdapter::dump(const string& fil Line 502  void MeshAdapter::dump(const string& fil
502     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
503    
504  #else // USE_NETCDF  #else // USE_NETCDF
505      Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with netCDF. Please contact your installation manager.");      setError(IO_ERROR, "MeshAdapter::dump: not configured with netCDF. Please contact your installation manager.");
506  #endif // USE_NETCDF  #endif // USE_NETCDF
507      checkFinleyError();      checkFinleyError();
508  }  }
# Line 638  int MeshAdapter::getDiracDeltaFunctionsC Line 630  int MeshAdapter::getDiracDeltaFunctionsC
630  //  //
631  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
632  {  {
633      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
634      int numDim=Finley_Mesh_getDim(mesh);      return mesh->getDim();
     checkFinleyError();  
     return numDim;  
635  }  }
636    
637  //  //
# Line 649  int MeshAdapter::getDim() const Line 639  int MeshAdapter::getDim() const
639  //  //
640  int MeshAdapter::getNumDataPointsGlobal() const  int MeshAdapter::getNumDataPointsGlobal() const
641  {  {
642      return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);      return m_finleyMesh.get()->Nodes->getGlobalNumNodes();
643  }  }
644    
645  //  //
# Line 660  pair<int,int> MeshAdapter::getDataShape( Line 650  pair<int,int> MeshAdapter::getDataShape(
650  {  {
651      int numDataPointsPerSample=0;      int numDataPointsPerSample=0;
652      int numSamples=0;      int numSamples=0;
653      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
654      switch (functionSpaceCode) {      switch (functionSpaceCode) {
655          case Nodes:          case Nodes:
656              numDataPointsPerSample=1;              numDataPointsPerSample=1;
657              numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);              numSamples=mesh->Nodes->getNumNodes();
658              break;              break;
659          case ReducedNodes:          case ReducedNodes:
660              numDataPointsPerSample=1;              numDataPointsPerSample=1;
661              numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);              numSamples=mesh->Nodes->getNumReducedNodes();
662              break;              break;
663          case Elements:          case Elements:
664              if (mesh->Elements!=NULL) {              if (mesh->Elements!=NULL) {
# Line 727  pair<int,int> MeshAdapter::getDataShape( Line 717  pair<int,int> MeshAdapter::getDataShape(
717          case DegreesOfFreedom:          case DegreesOfFreedom:
718              if (mesh->Nodes!=NULL) {              if (mesh->Nodes!=NULL) {
719                  numDataPointsPerSample=1;                  numDataPointsPerSample=1;
720                  numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);                  numSamples=mesh->Nodes->getNumDegreesOfFreedom();
721              }              }
722              break;              break;
723          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
724              if (mesh->Nodes!=NULL) {              if (mesh->Nodes!=NULL) {
725                  numDataPointsPerSample=1;                  numDataPointsPerSample=1;
726                  numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);                  numSamples=mesh->Nodes->getNumReducedDegreesOfFreedom();
727              }              }
728              break;              break;
729          default:          default:
# Line 761  void MeshAdapter::addPDEToSystem( Line 751  void MeshAdapter::addPDEToSystem(
751      if (!smat)      if (!smat)
752          throw FinleyAdapterException("finley only supports Paso system matrices.");          throw FinleyAdapterException("finley only supports Paso system matrices.");
753    
754      escriptDataC _rhs=rhs.getDataC();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC _A  =A.getDataC();  
     escriptDataC _B=B.getDataC();  
     escriptDataC _C=C.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _d_contact=d_contact.getDataC();  
     escriptDataC _y_contact=y_contact.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Finley_Mesh* mesh=m_finleyMesh.get();  
755      Paso_SystemMatrix* S = smat->getPaso_SystemMatrix();      Paso_SystemMatrix* S = smat->getPaso_SystemMatrix();
756        Assemble_PDE(mesh->Nodes, mesh->Elements, S, rhs, A, B, C, D, X, Y);
     Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, S, &_rhs, &_A, &_B, &_C,  
                         &_D, &_X, &_Y);  
757      checkFinleyError();      checkFinleyError();
758    
759      Finley_Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, &_rhs, 0, 0, 0,  
760                          &_d, 0, &_y);      Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs,
761                escript::Data(), escript::Data(), escript::Data(), d,
762                escript::Data(), y);
763      checkFinleyError();      checkFinleyError();
764    
765      Finley_Assemble_PDE(mesh->Nodes, mesh->ContactElements, S, &_rhs, 0, 0, 0,      Assemble_PDE(mesh->Nodes, mesh->ContactElements, S, rhs,
766              &_d_contact, 0, &_y_contact);              escript::Data(), escript::Data(), escript::Data(), d_contact,
767                escript::Data(), y_contact);
768      checkFinleyError();      checkFinleyError();
769    
770      Finley_Assemble_PDE(mesh->Nodes, mesh->Points, S, &_rhs, 0, 0, 0,      Assemble_PDE(mesh->Nodes, mesh->Points, S, rhs, escript::Data(),
771                          &_d_dirac, 0, &_y_dirac );              escript::Data(), escript::Data(), d_dirac, escript::Data(), y_dirac);
772      checkFinleyError();      checkFinleyError();
773  }  }
774    
# Line 799  void MeshAdapter::addPDEToLumpedSystem(e Line 776  void MeshAdapter::addPDEToLumpedSystem(e
776          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
777          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
778  {  {
779      escriptDataC _mat=mat.getDataC();      Mesh* mesh=m_finleyMesh.get();
780      escriptDataC _D=D.getDataC();      Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
     escriptDataC _d=d.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
   
     Finley_Mesh* mesh=m_finleyMesh.get();  
   
     Finley_Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, &_mat, &_D, useHRZ);  
781      checkFinleyError();      checkFinleyError();
782    
783      Finley_Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, &_mat, &_d, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->FaceElements, mat, d, useHRZ);
784      checkFinleyError();      checkFinleyError();
785    
786      Finley_Assemble_LumpedSystem(mesh->Nodes, mesh->Points, &_mat, &_d_dirac, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->Points, mat, d_dirac, useHRZ);
787      checkFinleyError();      checkFinleyError();
788  }  }
789    
# Line 823  void MeshAdapter::addPDEToRHS(escript::D Line 794  void MeshAdapter::addPDEToRHS(escript::D
794          const escript::Data& Y, const escript::Data& y,          const escript::Data& Y, const escript::Data& y,
795          const escript::Data& y_contact, const escript::Data& y_dirac) const          const escript::Data& y_contact, const escript::Data& y_dirac) const
796  {  {
797      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
798        Assemble_PDE(mesh->Nodes, mesh->Elements, 0, rhs,
799      escriptDataC _rhs=rhs.getDataC();              escript::Data(), escript::Data(), escript::Data(), escript::Data(),
800      escriptDataC _X=X.getDataC();              X, Y);
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _y_contact=y_contact.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
   
     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y);  
801      checkFinleyError();      checkFinleyError();
802    
803      Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y);      Assemble_PDE(mesh->Nodes, mesh->FaceElements, 0, rhs,
804                escript::Data(), escript::Data(), escript::Data(), escript::Data(),
805                escript::Data(), y);
806      checkFinleyError();      checkFinleyError();
807    
808      Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs, 0, 0, 0, 0, 0, &_y_contact);      Assemble_PDE(mesh->Nodes, mesh->ContactElements, 0, rhs,
809                escript::Data(), escript::Data(), escript::Data(),
810                escript::Data(), escript::Data(), y_contact);
811      checkFinleyError();      checkFinleyError();
812    
813      Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs, 0, 0, 0, 0, 0, &_y_dirac);      Assemble_PDE(mesh->Nodes, mesh->Points, 0, rhs,
814                escript::Data(), escript::Data(), escript::Data(), escript::Data(),
815                escript::Data(), y_dirac);
816      checkFinleyError();      checkFinleyError();
817  }  }
818    
# Line 861  void MeshAdapter::addPDEToTransportProbl Line 832  void MeshAdapter::addPDEToTransportProbl
832          throw FinleyAdapterException("finley only supports Paso transport problems.");          throw FinleyAdapterException("finley only supports Paso transport problems.");
833    
834      source.expand();      source.expand();
     escriptDataC _source=source.getDataC();  
     escriptDataC _M=M.getDataC();  
     escriptDataC _A=A.getDataC();  
     escriptDataC _B=B.getDataC();  
     escriptDataC _C=C.getDataC();  
     escriptDataC _D=D.getDataC();  
     escriptDataC _X=X.getDataC();  
     escriptDataC _Y=Y.getDataC();  
     escriptDataC _d=d.getDataC();  
     escriptDataC _y=y.getDataC();  
     escriptDataC _d_contact=d_contact.getDataC();  
     escriptDataC _y_contact=y_contact.getDataC();  
     escriptDataC _d_dirac=d_dirac.getDataC();  
     escriptDataC _y_dirac=y_dirac.getDataC();  
835    
836      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
837      Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();      Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
838    
839      Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix,      Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix, source,
840                          &_source, 0, 0, 0, &_M, 0, 0);                          escript::Data(), escript::Data(), escript::Data(),
841                            M, escript::Data(), escript::Data());
842      checkFinleyError();      checkFinleyError();
843    
844      Finley_Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->transport_matrix,      Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->transport_matrix,
845                          &_source, &_A, &_B, &_C, &_D, &_X, &_Y);                          source, A, B, C, D, X, Y);
846      checkFinleyError();      checkFinleyError();
847    
848      Finley_Assemble_PDE(mesh->Nodes, mesh->FaceElements, _tp->transport_matrix,      Assemble_PDE(mesh->Nodes, mesh->FaceElements, _tp->transport_matrix,
849                          &_source, 0, 0, 0, &_d, 0, &_y);                          source, escript::Data(), escript::Data(),
850                            escript::Data(), d, escript::Data(), y);
851      checkFinleyError();      checkFinleyError();
852    
853      Finley_Assemble_PDE(mesh->Nodes, mesh->ContactElements,      Assemble_PDE(mesh->Nodes, mesh->ContactElements,
854                          _tp->transport_matrix, &_source, 0, 0, 0, &_d_contact, 0, &_y_contact);                          _tp->transport_matrix, source, escript::Data(),
855                            escript::Data(), escript::Data(), d_contact,
856                            escript::Data(), y_contact);
857      checkFinleyError();      checkFinleyError();
858    
859      Finley_Assemble_PDE(mesh->Nodes, mesh->Points, _tp->transport_matrix,      Assemble_PDE(mesh->Nodes, mesh->Points, _tp->transport_matrix,
860                          &_source, 0, 0, 0, &_d_dirac, 0, &_y_dirac);                          source, escript::Data(), escript::Data(),
861                            escript::Data(), d_dirac, escript::Data(), y_dirac);
862      checkFinleyError();      checkFinleyError();
863  }  }
864    
# Line 912  void MeshAdapter::interpolateOnDomain(es Line 874  void MeshAdapter::interpolateOnDomain(es
874      if (targetDomain!=*this)      if (targetDomain!=*this)
875          throw FinleyAdapterException("Error - Illegal domain of interpolation target.");          throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
876    
877      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC _target=target.getDataC();  
     escriptDataC _in=in.getDataC();  
878      switch(in.getFunctionSpace().getTypeCode()) {      switch(in.getFunctionSpace().getTypeCode()) {
879          case Nodes:          case Nodes:
880              switch(target.getFunctionSpace().getTypeCode()) {              switch(target.getFunctionSpace().getTypeCode()) {
# Line 922  void MeshAdapter::interpolateOnDomain(es Line 882  void MeshAdapter::interpolateOnDomain(es
882                  case ReducedNodes:                  case ReducedNodes:
883                  case DegreesOfFreedom:                  case DegreesOfFreedom:
884                  case ReducedDegreesOfFreedom:                  case ReducedDegreesOfFreedom:
885                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      Assemble_CopyNodalData(mesh->Nodes, target, in);
886                      break;                      break;
887                  case Elements:                  case Elements:
888                  case ReducedElements:                  case ReducedElements:
889                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->Elements, in,target);
890                      break;                      break;
891                  case FaceElements:                  case FaceElements:
892                  case ReducedFaceElements:                  case ReducedFaceElements:
893                      Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
894                      break;                      break;
895                  case Points:                  case Points:
896                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
897                      break;                      break;
898                  case ContactElementsZero:                  case ContactElementsZero:
899                  case ReducedContactElementsZero:                  case ReducedContactElementsZero:
900                  case ContactElementsOne:                  case ContactElementsOne:
901                  case ReducedContactElementsOne:                  case ReducedContactElementsOne:
902                      Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
903                      break;                      break;
904                  default:                  default:
905                      stringstream temp;                      stringstream temp;
# Line 954  void MeshAdapter::interpolateOnDomain(es Line 914  void MeshAdapter::interpolateOnDomain(es
914                  case ReducedNodes:                  case ReducedNodes:
915                  case DegreesOfFreedom:                  case DegreesOfFreedom:
916                  case ReducedDegreesOfFreedom:                  case ReducedDegreesOfFreedom:
917                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      Assemble_CopyNodalData(mesh->Nodes, target, in);
918                      break;                      break;
919                  case Elements:                  case Elements:
920                  case ReducedElements:                  case ReducedElements:
921                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
922                      break;                      break;
923                  case FaceElements:                  case FaceElements:
924                  case ReducedFaceElements:                  case ReducedFaceElements:
925                      Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
926                      break;                      break;
927                  case Points:                  case Points:
928                      Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                      Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
929                      break;                      break;
930                  case ContactElementsZero:                  case ContactElementsZero:
931                  case ReducedContactElementsZero:                  case ReducedContactElementsZero:
932                      Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                  case ContactElementsOne:
933                    case ReducedContactElementsOne:
934                        Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
935                      break;                      break;
936                  default:                  default:
937                      stringstream temp;                      stringstream temp;
# Line 980  void MeshAdapter::interpolateOnDomain(es Line 942  void MeshAdapter::interpolateOnDomain(es
942              break;              break;
943          case Elements:          case Elements:
944              if (target.getFunctionSpace().getTypeCode()==Elements) {              if (target.getFunctionSpace().getTypeCode()==Elements) {
945                  Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);                  Assemble_CopyElementData(mesh->Elements, target, in);
946              } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {              } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
947                  Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);                  Assemble_AverageElementData(mesh->Elements, target, in);
948              } else {              } else {
949                  throw FinleyAdapterException("Error - No interpolation with data on elements possible.");                  throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
950              }              }
951              break;              break;
952          case ReducedElements:          case ReducedElements:
953              if (target.getFunctionSpace().getTypeCode()==ReducedElements) {              if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
954                  Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);                  Assemble_CopyElementData(mesh->Elements, target, in);
955              } else {              } else {
956                  throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");                  throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
957              }              }
958              break;              break;
959          case FaceElements:          case FaceElements:
960              if (target.getFunctionSpace().getTypeCode()==FaceElements) {              if (target.getFunctionSpace().getTypeCode()==FaceElements) {
961                  Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);                  Assemble_CopyElementData(mesh->FaceElements, target, in);
962              } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {              } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
963                  Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);                  Assemble_AverageElementData(mesh->FaceElements, target, in);
964              } else {              } else {
965                  throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");                  throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
966              }              }
967              break;              break;
968          case ReducedFaceElements:          case ReducedFaceElements:
969              if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {              if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
970                  Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);                  Assemble_CopyElementData(mesh->FaceElements, target, in);
971              } else {              } else {
972                  throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");                  throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
973              }              }
974              break;              break;
975          case Points:          case Points:
976              if (target.getFunctionSpace().getTypeCode()==Points) {              if (target.getFunctionSpace().getTypeCode()==Points) {
977                  Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);                  Assemble_CopyElementData(mesh->Points, target, in);
978              } else {              } else {
979                  throw FinleyAdapterException("Error - No interpolation with data on points possible.");                  throw FinleyAdapterException("Error - No interpolation with data on points possible.");
980              }              }
# Line 1020  void MeshAdapter::interpolateOnDomain(es Line 982  void MeshAdapter::interpolateOnDomain(es
982          case ContactElementsZero:          case ContactElementsZero:
983          case ContactElementsOne:          case ContactElementsOne:
984              if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {              if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
985                  Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);                  Assemble_CopyElementData(mesh->ContactElements, target, in);
986              } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {              } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
987                  Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);                  Assemble_AverageElementData(mesh->ContactElements, target, in);
988              } else {              } else {
989                  throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");                  throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
990              }              }
# Line 1030  void MeshAdapter::interpolateOnDomain(es Line 992  void MeshAdapter::interpolateOnDomain(es
992          case ReducedContactElementsZero:          case ReducedContactElementsZero:
993          case ReducedContactElementsOne:          case ReducedContactElementsOne:
994              if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {              if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
995                  Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);                  Assemble_CopyElementData(mesh->ContactElements, target, in);
996              } else {              } else {
997                  throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");                  throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
998              }              }
# Line 1039  void MeshAdapter::interpolateOnDomain(es Line 1001  void MeshAdapter::interpolateOnDomain(es
1001              switch(target.getFunctionSpace().getTypeCode()) {              switch(target.getFunctionSpace().getTypeCode()) {
1002                  case ReducedDegreesOfFreedom:                  case ReducedDegreesOfFreedom:
1003                  case DegreesOfFreedom:                  case DegreesOfFreedom:
1004                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      Assemble_CopyNodalData(mesh->Nodes, target, in);
1005                      break;                      break;
1006    
1007                  case Nodes:                  case Nodes:
1008                  case ReducedNodes:                  case ReducedNodes:
1009                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1010                          escript::Data temp=escript::Data(in);                          escript::Data in2=escript::Data(in);
1011                          temp.expand();                          in2.expand();
1012                          escriptDataC _in2 = temp.getDataC();                          Assemble_CopyNodalData(mesh->Nodes, target, in2);
                         Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
1013                      } else {                      } else {
1014                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                          Assemble_CopyNodalData(mesh->Nodes, target, in);
1015                      }                      }
1016                      break;                      break;
1017                  case Elements:                  case Elements:
1018                  case ReducedElements:                  case ReducedElements:
1019                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1020                          escript::Data temp=escript::Data(in, continuousFunction(*this));                          escript::Data in2=escript::Data(in, continuousFunction(*this));
1021                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
1022                      } else {                      } else {
1023                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1024                      }                      }
1025                      break;                      break;
1026                  case FaceElements:                  case FaceElements:
1027                  case ReducedFaceElements:                  case ReducedFaceElements:
1028                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1029                          escript::Data temp=escript::Data(in, continuousFunction(*this));                          escript::Data in2=escript::Data(in, continuousFunction(*this));
1030                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
1031                      } else {                      } else {
1032                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1033                      }                      }
1034                      break;                      break;
1035                  case Points:                  case Points:
1036                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1037                          //escript::Data temp=escript::Data(in, continuousFunction(*this) );                          //escript::Data in2=escript::Data(in, continuousFunction(*this) );
                         //escriptDataC _in2 = temp.getDataC();  
1038                      } else {                      } else {
1039                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1040                      }                      }
1041                      break;                      break;
1042                  case ContactElementsZero:                  case ContactElementsZero:
# Line 1086  void MeshAdapter::interpolateOnDomain(es Line 1044  void MeshAdapter::interpolateOnDomain(es
1044                  case ReducedContactElementsZero:                  case ReducedContactElementsZero:
1045                  case ReducedContactElementsOne:                  case ReducedContactElementsOne:
1046                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1047                          escript::Data temp=escript::Data(in, continuousFunction(*this));                          escript::Data in2=escript::Data(in, continuousFunction(*this));
1048                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
1049                      } else {                      } else {
1050                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1051                      }                      }
1052                      break;                      break;
1053                  default:                  default:
1054                      stringstream temp;                      stringstream temp;
1055                      temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();                      temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1056                      throw FinleyAdapterException(temp.str());                      throw FinleyAdapterException(temp.str());
                     break;  
1057              }              }
1058              break;              break;
1059          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
1060              switch(target.getFunctionSpace().getTypeCode()) {              switch(target.getFunctionSpace().getTypeCode()) {
1061                  case Nodes:                  case Nodes:
1062                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
                     break;  
1063                  case ReducedNodes:                  case ReducedNodes:
1064                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1065                          escript::Data temp=escript::Data(in);                          escript::Data in2=escript::Data(in);
1066                          temp.expand();                          in2.expand();
1067                          escriptDataC _in2 = temp.getDataC();                          Assemble_CopyNodalData(mesh->Nodes, target, in2);
                         Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
1068                      } else {                      } else {
1069                          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                          Assemble_CopyNodalData(mesh->Nodes, target, in);
1070                      }                      }
1071                      break;                      break;
1072                  case DegreesOfFreedom:                  case DegreesOfFreedom:
1073                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");                      throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1074                      break;                      break;
1075                  case ReducedDegreesOfFreedom:                  case ReducedDegreesOfFreedom:
1076                      Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                      Assemble_CopyNodalData(mesh->Nodes, target, in);
1077                      break;                      break;
1078                  case Elements:                  case Elements:
1079                  case ReducedElements:                  case ReducedElements:
1080                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1081                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );                          escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1082                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->Elements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
1083                      } else {                      } else {
1084                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->Elements, in, target);
1085                      }                      }
1086                      break;                      break;
1087                  case FaceElements:                  case FaceElements:
1088                  case ReducedFaceElements:                  case ReducedFaceElements:
1089                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1090                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this) );                          escript::Data in2=escript::Data(in, reducedContinuousFunction(*this) );
1091                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
1092                      } else {                      } else {
1093                          Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->FaceElements, in, target);
1094                      }                      }
1095                      break;                      break;
1096                  case Points:                  case Points:
1097                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1098                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this));                          escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1099                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->Points, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);  
1100                      } else {                      } else {
1101                          Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->Points, in, target);
1102                      }                      }
1103                      break;                      break;
1104                  case ContactElementsZero:                  case ContactElementsZero:
# Line 1155  void MeshAdapter::interpolateOnDomain(es Line 1106  void MeshAdapter::interpolateOnDomain(es
1106                  case ReducedContactElementsZero:                  case ReducedContactElementsZero:
1107                  case ReducedContactElementsOne:                  case ReducedContactElementsOne:
1108                      if (getMPISize()>1) {                      if (getMPISize()>1) {
1109                          escript::Data temp=escript::Data(in, reducedContinuousFunction(*this));                          escript::Data in2=escript::Data(in, reducedContinuousFunction(*this));
1110                          escriptDataC _in2 = temp.getDataC();                          Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in2, target);
                         Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
1111                      } else {                      } else {
1112                          Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                          Assemble_interpolate(mesh->Nodes, mesh->ContactElements, in, target);
1113                      }                      }
1114                      break;                      break;
1115                  default:                  default:
# Line 1186  void MeshAdapter::setToX(escript::Data& Line 1136  void MeshAdapter::setToX(escript::Data&
1136      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1137      if (argDomain!=*this)      if (argDomain!=*this)
1138          throw FinleyAdapterException("Error - Illegal domain of data point locations");          throw FinleyAdapterException("Error - Illegal domain of data point locations");
1139      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1140      // in case of appropriate function space we can do the job directly:      // in case of appropriate function space we can do the job directly:
1141      if (arg.getFunctionSpace().getTypeCode()==Nodes) {      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1142          escriptDataC _arg=arg.getDataC();          Assemble_NodeCoordinates(mesh->Nodes, arg);
         Finley_Assemble_NodeCoordinates(mesh->Nodes, &_arg);  
1143      } else {      } else {
1144          escript::Data tmp_data=Vector(0., continuousFunction(*this), true);          escript::Data tmp_data=Vector(0., continuousFunction(*this), true);
1145          escriptDataC _tmp_data=tmp_data.getDataC();          Assemble_NodeCoordinates(mesh->Nodes, tmp_data);
         Finley_Assemble_NodeCoordinates(mesh->Nodes, &_tmp_data);  
1146          // this is then interpolated onto arg:          // this is then interpolated onto arg:
1147          interpolateOnDomain(arg, tmp_data);          interpolateOnDomain(arg, tmp_data);
1148      }      }
# Line 1209  void MeshAdapter::setToNormal(escript::D Line 1157  void MeshAdapter::setToNormal(escript::D
1157      const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));      const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1158      if (normalDomain!=*this)      if (normalDomain!=*this)
1159          throw FinleyAdapterException("Error - Illegal domain of normal locations");          throw FinleyAdapterException("Error - Illegal domain of normal locations");
1160      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC _normal=normal.getDataC();  
1161      switch(normal.getFunctionSpace().getTypeCode()) {      switch(normal.getFunctionSpace().getTypeCode()) {
1162          case Nodes:          case Nodes:
1163              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
# Line 1226  void MeshAdapter::setToNormal(escript::D Line 1173  void MeshAdapter::setToNormal(escript::D
1173              break;              break;
1174          case FaceElements:          case FaceElements:
1175          case ReducedFaceElements:          case ReducedFaceElements:
1176              Finley_Assemble_setNormal(mesh->Nodes,mesh->FaceElements,&_normal);              Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
1177              break;              break;
1178          case Points:          case Points:
1179              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
# Line 1235  void MeshAdapter::setToNormal(escript::D Line 1182  void MeshAdapter::setToNormal(escript::D
1182          case ContactElementsZero:          case ContactElementsZero:
1183          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1184          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1185              Finley_Assemble_setNormal(mesh->Nodes,mesh->ContactElements,&_normal);              Assemble_getNormal(mesh->Nodes, mesh->ContactElements, normal);
1186              break;              break;
1187          case DegreesOfFreedom:          case DegreesOfFreedom:
1188              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
# Line 1255  void MeshAdapter::setToNormal(escript::D Line 1202  void MeshAdapter::setToNormal(escript::D
1202  //  //
1203  // interpolates data to other domain  // interpolates data to other domain
1204  //  //
1205  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target, const escript::Data& source) const
1206  {  {
1207      escript::const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();      throw FinleyAdapterException("Error - Finley does not allow interpolation across domains.");
     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());  
     if (targetDomain!=this)  
         throw FinleyAdapterException("Error - Illegal domain of interpolation target");  
   
     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");  
1208  }  }
1209    
1210  //  //
# Line 1275  void MeshAdapter::setToIntegrals(vector< Line 1217  void MeshAdapter::setToIntegrals(vector<
1217          throw FinleyAdapterException("Error - Illegal domain of integration kernel");          throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1218    
1219      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
1220      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC _temp;  
     escript::Data temp;  
     escriptDataC _arg=arg.getDataC();  
1221      switch(arg.getFunctionSpace().getTypeCode()) {      switch(arg.getFunctionSpace().getTypeCode()) {
1222          case Nodes:          case Nodes:
1223              temp=escript::Data(arg, escript::function(*this));              {
1224              _temp=temp.getDataC();                  escript::Data temp(arg, escript::function(*this));
1225              Finley_Assemble_integrate(mesh->Nodes, mesh->Elements, &_temp, &integrals[0]);                  Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1226                }
1227              break;              break;
1228          case ReducedNodes:          case ReducedNodes:
1229              temp=escript::Data( arg, escript::function(*this) );              {
1230              _temp=temp.getDataC();                  escript::Data temp(arg, escript::function(*this));
1231              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);                  Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1232                }
1233              break;              break;
1234          case Elements:          case Elements:
1235          case ReducedElements:          case ReducedElements:
1236              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);              Assemble_integrate(mesh->Nodes, mesh->Elements, arg, &integrals[0]);
1237              break;              break;
1238          case FaceElements:          case FaceElements:
1239          case ReducedFaceElements:          case ReducedFaceElements:
1240              Finley_Assemble_integrate(mesh->Nodes,mesh->FaceElements,&_arg,&integrals[0]);              Assemble_integrate(mesh->Nodes, mesh->FaceElements, arg, &integrals[0]);
1241              break;              break;
1242          case Points:          case Points:
1243              throw FinleyAdapterException("Error - Integral of data on points is not supported.");              throw FinleyAdapterException("Error - Integral of data on points is not supported.");
# Line 1305  void MeshAdapter::setToIntegrals(vector< Line 1246  void MeshAdapter::setToIntegrals(vector<
1246          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1247          case ContactElementsOne:          case ContactElementsOne:
1248          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1249              Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);              Assemble_integrate(mesh->Nodes, mesh->ContactElements, arg, &integrals[0]);
1250              break;              break;
1251          case DegreesOfFreedom:          case DegreesOfFreedom:
1252          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
1253              temp=escript::Data(arg, escript::function(*this));              {
1254              _temp=temp.getDataC();                  escript::Data temp(arg, escript::function(*this));
1255              Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);                  Assemble_integrate(mesh->Nodes, mesh->Elements, temp, &integrals[0]);
1256                }
1257              break;              break;
1258          default:          default:
1259              stringstream temp;              stringstream temp;
# Line 1335  void MeshAdapter::setToGradient(escript: Line 1277  void MeshAdapter::setToGradient(escript:
1277      if (gradDomain!=*this)      if (gradDomain!=*this)
1278          throw FinleyAdapterException("Error - Illegal domain of gradient");          throw FinleyAdapterException("Error - Illegal domain of gradient");
1279    
1280      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1281      escriptDataC _grad=grad.getDataC();      escript::Data nodeData;
     escriptDataC nodeDataC;  
     escript::Data temp;  
1282      if (getMPISize()>1) {      if (getMPISize()>1) {
1283          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
1284              temp=escript::Data(arg, continuousFunction(*this));              nodeData=escript::Data(arg, continuousFunction(*this));
             nodeDataC = temp.getDataC();  
1285          } else if(arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {          } else if(arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) {
1286              temp=escript::Data(arg, reducedContinuousFunction(*this));              nodeData=escript::Data(arg, reducedContinuousFunction(*this));
             nodeDataC = temp.getDataC();  
1287          } else {          } else {
1288              nodeDataC = arg.getDataC();              nodeData = arg;
1289          }          }
1290      } else {      } else {
1291          nodeDataC = arg.getDataC();          nodeData = arg;
1292      }      }
1293      switch(grad.getFunctionSpace().getTypeCode()) {      switch(grad.getFunctionSpace().getTypeCode()) {
1294          case Nodes:          case Nodes:
# Line 1361  void MeshAdapter::setToGradient(escript: Line 1299  void MeshAdapter::setToGradient(escript:
1299              break;              break;
1300          case Elements:          case Elements:
1301          case ReducedElements:          case ReducedElements:
1302              Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&_grad,&nodeDataC);              Assemble_gradient(mesh->Nodes, mesh->Elements, grad, nodeData);
1303              break;              break;
1304          case FaceElements:          case FaceElements:
1305          case ReducedFaceElements:          case ReducedFaceElements:
1306              Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&_grad,&nodeDataC);              Assemble_gradient(mesh->Nodes, mesh->FaceElements, grad, nodeData);
1307              break;              break;
1308          case Points:          case Points:
1309              throw FinleyAdapterException("Error - Gradient at points is not supported.");              throw FinleyAdapterException("Error - Gradient at points is not supported.");
# Line 1374  void MeshAdapter::setToGradient(escript: Line 1312  void MeshAdapter::setToGradient(escript:
1312          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1313          case ContactElementsOne:          case ContactElementsOne:
1314          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1315              Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&_grad,&nodeDataC);              Assemble_gradient(mesh->Nodes, mesh->ContactElements, grad, nodeData);
1316              break;              break;
1317          case DegreesOfFreedom:          case DegreesOfFreedom:
1318              throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");              throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
# Line 1396  void MeshAdapter::setToGradient(escript: Line 1334  void MeshAdapter::setToGradient(escript:
1334  //  //
1335  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1336  {  {
1337      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC tmp=size.getDataC();  
1338      switch(size.getFunctionSpace().getTypeCode()) {      switch(size.getFunctionSpace().getTypeCode()) {
1339          case Nodes:          case Nodes:
1340              throw FinleyAdapterException("Error - Size of nodes is not supported.");              throw FinleyAdapterException("Error - Size of nodes is not supported.");
# Line 1407  void MeshAdapter::setToSize(escript::Dat Line 1344  void MeshAdapter::setToSize(escript::Dat
1344              break;              break;
1345          case Elements:          case Elements:
1346          case ReducedElements:          case ReducedElements:
1347              Finley_Assemble_getSize(mesh->Nodes,mesh->Elements,&tmp);              Assemble_getSize(mesh->Nodes, mesh->Elements, size);
1348              break;              break;
1349          case FaceElements:          case FaceElements:
1350          case ReducedFaceElements:          case ReducedFaceElements:
1351              Finley_Assemble_getSize(mesh->Nodes,mesh->FaceElements,&tmp);              Assemble_getSize(mesh->Nodes, mesh->FaceElements, size);
1352              break;              break;
1353          case Points:          case Points:
1354              throw FinleyAdapterException("Error - Size of point elements is not supported.");              throw FinleyAdapterException("Error - Size of point elements is not supported.");
# Line 1420  void MeshAdapter::setToSize(escript::Dat Line 1357  void MeshAdapter::setToSize(escript::Dat
1357          case ContactElementsOne:          case ContactElementsOne:
1358          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1359          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1360              Finley_Assemble_getSize(mesh->Nodes,mesh->ContactElements,&tmp);              Assemble_getSize(mesh->Nodes,mesh->ContactElements,size);
1361              break;              break;
1362          case DegreesOfFreedom:          case DegreesOfFreedom:
1363              throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");              throw FinleyAdapterException("Error - Size of degrees of freedom is not supported.");
# Line 1442  void MeshAdapter::setToSize(escript::Dat Line 1379  void MeshAdapter::setToSize(escript::Dat
1379  //  //
1380  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1381  {  {
1382      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC tmp;  
1383      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1384      if (newDomain!=*this)      if (newDomain!=*this)
1385          throw FinleyAdapterException("Error - Illegal domain of new point locations");          throw FinleyAdapterException("Error - Illegal domain of new point locations");
1386      if (new_x.getFunctionSpace() == continuousFunction(*this)) {      if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1387          tmp = new_x.getDataC();          mesh->setCoordinates(new_x);
         Finley_Mesh_setCoordinates(mesh,&tmp);  
1388      } else {      } else {
1389          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
         //escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );  
         //tmp = new_x_inter.getDataC();  
         //Finley_Mesh_setCoordinates(mesh,&tmp);  
1390      }      }
1391      checkFinleyError();      checkFinleyError();
1392  }  }
# Line 1463  bool MeshAdapter::ownSample(int fs_code, Line 1395  bool MeshAdapter::ownSample(int fs_code,
1395  {  {
1396      if (getMPISize() > 1) {      if (getMPISize() > 1) {
1397  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1398          index_t myFirstNode=0, myLastNode=0, k=0;          int myFirstNode=0, myLastNode=0, k=0;
1399          index_t* globalNodeIndex=0;          int* globalNodeIndex=0;
1400          Finley_Mesh* mesh_p=m_finleyMesh.get();          Mesh* mesh_p=m_finleyMesh.get();
1401          /*          /*
1402           * this method is only used by saveDataCSV which would use the returned           * this method is only used by saveDataCSV which would use the returned
1403           * values for reduced nodes wrongly so this case is disabled for now           * values for reduced nodes wrongly so this case is disabled for now
1404          if (fs_code == FINLEY_REDUCED_NODES) {          if (fs_code == FINLEY_REDUCED_NODES) {
1405              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);              myFirstNode = NodeFile_getFirstReducedNode(mesh_p->Nodes);
1406              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);              myLastNode = NodeFile_getLastReducedNode(mesh_p->Nodes);
1407              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);              globalNodeIndex = NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1408          } else          } else
1409          */          */
1410          if (fs_code == FINLEY_NODES) {          if (fs_code == FINLEY_NODES) {
1411              myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);              myFirstNode = mesh_p->Nodes->getFirstNode();
1412              myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);              myLastNode = mesh_p->Nodes->getLastNode();
1413              globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);              globalNodeIndex = mesh_p->Nodes->borrowGlobalNodesIndex();
1414          } else {          } else {
1415              throw FinleyAdapterException("Unsupported function space type for ownSample()");              throw FinleyAdapterException("Unsupported function space type for ownSample()");
1416          }          }
1417    
1418          k=globalNodeIndex[id];          k=globalNodeIndex[id];
1419          return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );          return ((myFirstNode <= k) && (k < myLastNode));
1420  #endif  #endif
1421      }      }
1422      return true;      return true;
# Line 1523  escript::ASM_ptr MeshAdapter::newSystemM Line 1455  escript::ASM_ptr MeshAdapter::newSystemM
1455      }      }
1456    
1457      // generate matrix:      // generate matrix:
1458      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(      paso::SystemMatrixPattern* fsystemMatrixPattern=
1459              getFinley_Mesh(), reduceRowOrder, reduceColOrder);          getFinley_Mesh()->getPattern(reduceRowOrder, reduceColOrder);
1460      checkFinleyError();      checkFinleyError();
1461      Paso_SystemMatrix* fsystemMatrix;      Paso_SystemMatrix* fsystemMatrix;
1462      const int trilinos = 0;      const int trilinos = 0;
# Line 1537  escript::ASM_ptr MeshAdapter::newSystemM Line 1469  escript::ASM_ptr MeshAdapter::newSystemM
1469                  row_blocksize, column_blocksize, FALSE);                  row_blocksize, column_blocksize, FALSE);
1470      }      }
1471      checkPasoError();      checkPasoError();
1472      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      paso::SystemMatrixPattern_free(fsystemMatrixPattern);
1473      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1474      return escript::ASM_ptr(sma);      return escript::ASM_ptr(sma);
1475  }  }
# Line 1562  escript::ATP_ptr MeshAdapter::newTranspo Line 1494  escript::ATP_ptr MeshAdapter::newTranspo
1494      }      }
1495    
1496      // generate transport problem:      // generate transport problem:
1497      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(      paso::SystemMatrixPattern* fsystemMatrixPattern=
1498              getFinley_Mesh(), reduceOrder, reduceOrder);          getFinley_Mesh()->getPattern(reduceOrder, reduceOrder);
1499      checkFinleyError();      checkFinleyError();
1500      Paso_TransportProblem* transportProblem;      Paso_TransportProblem* transportProblem;
1501      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1502      checkPasoError();      checkPasoError();
1503      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      paso::SystemMatrixPattern_free(fsystemMatrixPattern);
1504      TransportProblemAdapter* tpa=new TransportProblemAdapter(      TransportProblemAdapter* tpa=new TransportProblemAdapter(
1505              transportProblem, blocksize, functionspace);              transportProblem, blocksize, functionspace);
1506      return escript::ATP_ptr(tpa);      return escript::ATP_ptr(tpa);
# Line 1914  bool MeshAdapter::operator!=(const escri Line 1846  bool MeshAdapter::operator!=(const escri
1846    
1847  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1848  {  {
1849      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1850      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1851                 package, symmetry, mesh->MPIInfo);                 package, symmetry, mesh->MPIInfo);
1852  }  }
1853    
1854  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1855  {  {
1856      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1857      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1858                  package, symmetry, mesh->MPIInfo);                  package, symmetry, mesh->MPIInfo);
1859  }  }
# Line 1944  escript::Data MeshAdapter::getSize() con Line 1876  escript::Data MeshAdapter::getSize() con
1876  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1877  {  {
1878      int *out = NULL;      int *out = NULL;
1879      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1880      switch (functionSpaceType) {      switch (functionSpaceType) {
1881          case Nodes:          case Nodes:
1882              out=mesh->Nodes->Id;              out=mesh->Nodes->Id;
# Line 1986  const int* MeshAdapter::borrowSampleRefe Line 1918  const int* MeshAdapter::borrowSampleRefe
1918  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1919  {  {
1920      int out=0;      int out=0;
1921      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1922      switch (functionSpaceType) {      switch (functionSpaceType) {
1923          case Nodes:          case Nodes:
1924              out=mesh->Nodes->Tag[sampleNo];              out=mesh->Nodes->Tag[sampleNo];
# Line 2029  int MeshAdapter::getTagFromSampleNo(int Line 1961  int MeshAdapter::getTagFromSampleNo(int
1961    
1962  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1963  {  {
1964      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     escriptDataC tmp=mask.getDataC();  
1965      switch(functionSpaceType) {      switch(functionSpaceType) {
1966          case Nodes:          case Nodes:
1967              Finley_NodeFile_setTags(mesh->Nodes,newTag,&tmp);              mesh->Nodes->setTags(newTag, mask);
1968              break;              break;
1969          case ReducedNodes:          case ReducedNodes:
1970              throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
             break;  
1971          case DegreesOfFreedom:          case DegreesOfFreedom:
1972              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
             break;  
1973          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
1974              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
             break;  
1975          case Elements:          case Elements:
1976          case ReducedElements:          case ReducedElements:
1977              Finley_ElementFile_setTags(mesh->Elements,newTag,&tmp);              mesh->Elements->setTags(newTag, mask);
1978              break;              break;
1979          case FaceElements:          case FaceElements:
1980          case ReducedFaceElements:          case ReducedFaceElements:
1981              Finley_ElementFile_setTags(mesh->FaceElements,newTag,&tmp);              mesh->FaceElements->setTags(newTag, mask);
1982              break;              break;
1983          case Points:          case Points:
1984              Finley_ElementFile_setTags(mesh->Points,newTag,&tmp);              mesh->Points->setTags(newTag, mask);
1985              break;              break;
1986          case ContactElementsZero:          case ContactElementsZero:
1987          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1988          case ContactElementsOne:          case ContactElementsOne:
1989          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1990              Finley_ElementFile_setTags(mesh->ContactElements,newTag,&tmp);              mesh->ContactElements->setTags(newTag, mask);
1991              break;              break;
1992          default:          default:
1993              stringstream temp;              stringstream temp;
# Line 2071  void MeshAdapter::setTags(const int func Line 1999  void MeshAdapter::setTags(const int func
1999    
2000  void MeshAdapter::setTagMap(const string& name, int tag)  void MeshAdapter::setTagMap(const string& name, int tag)
2001  {  {
2002      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2003      Finley_Mesh_addTagMap(mesh, name.c_str(), tag);      mesh->addTagMap(name.c_str(), tag);
2004      checkFinleyError();      checkFinleyError();
2005  }  }
2006    
2007  int MeshAdapter::getTag(const string& name) const  int MeshAdapter::getTag(const string& name) const
2008  {  {
2009      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2010      int tag = Finley_Mesh_getTag(mesh, name.c_str());      int tag = mesh->getTag(name.c_str());
2011      checkFinleyError();      checkFinleyError();
2012      return tag;      return tag;
2013  }  }
2014    
2015  bool MeshAdapter::isValidTagName(const string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2016  {  {
2017      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2018      return Finley_Mesh_isValidTagName(mesh, name.c_str());      return mesh->isValidTagName(name.c_str());
2019  }  }
2020    
2021  string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2022  {  {
2023      stringstream temp;      stringstream temp;
2024      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2025      Finley_TagMap* tag_map=mesh->TagMap;      TagMap::const_iterator it = mesh->tagMap.begin();
2026      while (tag_map) {      while (it != mesh->tagMap.end()) {
2027          temp << tag_map->name;          temp << it->first;
2028          tag_map=tag_map->next;          ++it;
2029          if (tag_map) temp << ", ";          if (it != mesh->tagMap.end())
2030                temp << ", ";
2031      }      }
2032      return temp.str();      return temp.str();
2033  }  }
2034    
2035  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2036  {  {
2037      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     dim_t numTags=0;  
2038      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2039          case Nodes:          case Nodes:
2040              numTags=mesh->Nodes->numTagsInUse;              return mesh->Nodes->tagsInUse.size();
             break;  
2041          case ReducedNodes:          case ReducedNodes:
2042              throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
             break;  
2043          case DegreesOfFreedom:          case DegreesOfFreedom:
2044              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
             break;  
2045          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
2046              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
             break;  
2047          case Elements:          case Elements:
2048          case ReducedElements:          case ReducedElements:
2049              numTags=mesh->Elements->numTagsInUse;              return mesh->Elements->tagsInUse.size();
             break;  
2050          case FaceElements:          case FaceElements:
2051          case ReducedFaceElements:          case ReducedFaceElements:
2052              numTags=mesh->FaceElements->numTagsInUse;              return mesh->FaceElements->tagsInUse.size();
             break;  
2053          case Points:          case Points:
2054              numTags=mesh->Points->numTagsInUse;              return mesh->Points->tagsInUse.size();
             break;  
2055          case ContactElementsZero:          case ContactElementsZero:
2056          case ReducedContactElementsZero:          case ReducedContactElementsZero:
2057          case ContactElementsOne:          case ContactElementsOne:
2058          case ReducedContactElementsOne:          case ReducedContactElementsOne:
2059              numTags=mesh->ContactElements->numTagsInUse;              return mesh->ContactElements->tagsInUse.size();
             break;  
2060          default:          default:
2061              stringstream temp;              stringstream ss;
2062              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              ss << "Finley does not know anything about function space type "
2063              throw FinleyAdapterException(temp.str());                   << functionSpaceCode;
2064                throw FinleyAdapterException(ss.str());
2065      }      }
2066      return numTags;      return 0;
2067  }  }
2068    
2069  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2070  {  {
2071      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
     index_t* tags=NULL;  
2072      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2073          case Nodes:          case Nodes:
2074              tags=mesh->Nodes->tagsInUse;              return &mesh->Nodes->tagsInUse[0];
             break;  
2075          case ReducedNodes:          case ReducedNodes:
2076              throw FinleyAdapterException("Error - ReducedNodes does not support tags");              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
             break;  
2077          case DegreesOfFreedom:          case DegreesOfFreedom:
2078              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
             break;  
2079          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
2080              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
             break;  
2081          case Elements:          case Elements:
2082          case ReducedElements:          case ReducedElements:
2083              tags=mesh->Elements->tagsInUse;              return &mesh->Elements->tagsInUse[0];
             break;  
2084          case FaceElements:          case FaceElements:
2085          case ReducedFaceElements:          case ReducedFaceElements:
2086              tags=mesh->FaceElements->tagsInUse;              return &mesh->FaceElements->tagsInUse[0];
             break;  
2087          case Points:          case Points:
2088              tags=mesh->Points->tagsInUse;              return &mesh->Points->tagsInUse[0];
             break;  
2089          case ContactElementsZero:          case ContactElementsZero:
2090          case ReducedContactElementsZero:          case ReducedContactElementsZero:
2091          case ContactElementsOne:          case ContactElementsOne:
2092          case ReducedContactElementsOne:          case ReducedContactElementsOne:
2093              tags=mesh->ContactElements->tagsInUse;              return &mesh->ContactElements->tagsInUse[0];
             break;  
2094          default:          default:
2095              stringstream temp;              stringstream temp;
2096              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;              temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2097              throw FinleyAdapterException(temp.str());              throw FinleyAdapterException(temp.str());
2098      }      }
2099      return tags;      return NULL;
2100  }  }
2101    
2102    
# Line 2209  bool MeshAdapter::canTag(int functionSpa Line 2121  bool MeshAdapter::canTag(int functionSpa
2121    
2122  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2123  {  {
2124      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2125      return Finley_Mesh_getStatus(mesh);      return mesh->getStatus();
2126  }  }
2127    
2128  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2129  {  {
2130      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2131      int order =-1;      int order =-1;
2132      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2133          case Nodes:          case Nodes:
# Line 2252  bool MeshAdapter::supportsContactElement Line 2164  bool MeshAdapter::supportsContactElement
2164      return true;      return true;
2165  }  }
2166    
2167  ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id,  escript::Data MeshAdapter::randomFill(const escript::DataTypes::ShapeType& shape,
2168          index_t order, index_t reducedOrder)         const escript::FunctionSpace& what, long seed, const boost::python::tuple& filter) const
2169  {  {
2170      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);      escript::Data towipe(0, shape, what, true);
2171        // since we just made this object, no sharing is possible and we don't need to check for
2172        // exlusive write
2173        escript::DataTypes::ValueType& dv=towipe.getExpandedVectorReference();
2174        const size_t dvsize=dv.size();
2175        esysUtils::randomFillArray(seed, &(dv[0]), dvsize);
2176        return towipe;  
2177  }  }
2178    
 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()  
 {  
     Finley_ReferenceElementSet_dealloc(m_refSet);  
 }  
2179    
2180  void MeshAdapter::addDiracPoints(const vector<double>& points,  void MeshAdapter::addDiracPoints(const vector<double>& points,
2181                                   const vector<int>& tags) const                                   const vector<int>& tags) const
# Line 2270  void MeshAdapter::addDiracPoints(const v Line 2184  void MeshAdapter::addDiracPoints(const v
2184      const int dim = getDim();      const int dim = getDim();
2185      int numPoints=points.size()/dim;      int numPoints=points.size()/dim;
2186      int numTags=tags.size();      int numTags=tags.size();
2187      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2188    
2189      if ( points.size() % dim != 0 ) {      if ( points.size() % dim != 0 ) {
2190          throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");          throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
# Line 2279  void MeshAdapter::addDiracPoints(const v Line 2193  void MeshAdapter::addDiracPoints(const v
2193      if ((numTags > 0) && (numPoints != numTags))      if ((numTags > 0) && (numPoints != numTags))
2194          throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");          throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2195    
2196      Finley_Mesh_addPoints(mesh, numPoints, &points[0], &tags[0]);      mesh->addPoints(numPoints, &points[0], &tags[0]);
2197      checkFinleyError();      checkFinleyError();
2198  }  }
2199    
# Line 2288  void MeshAdapter::addDiracPoints(const v Line 2202  void MeshAdapter::addDiracPoints(const v
2202  //       const int dim = getDim();  //       const int dim = getDim();
2203  //       int numPoints=boost::python::extract<int>(points.attr("__len__")());  //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2204  //       int numTags=boost::python::extract<int>(tags.attr("__len__")());  //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2205  //       Finley_Mesh* mesh=m_finleyMesh.get();  //       Mesh* mesh=m_finleyMesh.get();
2206  //  //
2207  //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )  //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2208  //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");  //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");

Legend:
Removed from v.4408  
changed lines
  Added in v.4800

  ViewVC Help
Powered by ViewVC 1.1.26