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

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

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

revision 1628 by phornby, Fri Jul 11 13:12:46 2008 UTC revision 2548 by jfenwick, Mon Jul 20 06:20:06 2009 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2009 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  #include "MeshAdapter.h"  #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
# Line 19  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21    #ifdef PASO_MPI
22    #include <mpi.h>
23    #include "paso/Paso_MPI.h"
24    #endif
25  extern "C" {  extern "C" {
26  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
27  }  }
 #include <vector>  
   
 #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  
28    
29  using namespace std;  using namespace std;
30  using namespace escript;  using namespace escript;
# Line 83  int MeshAdapter::getMPIRank() const Line 83  int MeshAdapter::getMPIRank() const
83  {  {
84     return m_finleyMesh.get()->MPIInfo->rank;     return m_finleyMesh.get()->MPIInfo->rank;
85  }  }
86    void MeshAdapter::MPIBarrier() const
87    {
88    #ifdef PASO_MPI
89       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
90    #endif
91       return;
92    }
93    bool MeshAdapter::onMasterProcessor() const
94    {
95       return m_finleyMesh.get()->MPIInfo->rank == 0;
96    }
97    
98    
99  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
# Line 98  void MeshAdapter::write(const std::strin Line 109  void MeshAdapter::write(const std::strin
109     TMPMEMFREE(fName);     TMPMEMFREE(fName);
110  }  }
111    
112  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
113  {  {
114     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
115  }  }
# Line 106  void MeshAdapter::Print_Mesh_Info(const Line 117  void MeshAdapter::Print_Mesh_Info(const
117  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const std::string& fileName) const
118  {  {
119  #ifdef USE_NETCDF  #ifdef USE_NETCDF
120     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
121     NcVar *ids, *data;     NcVar *ids;
122     int *int_ptr;     int *int_ptr;
123     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
124     Finley_TagMap* tag_map;     Finley_TagMap* tag_map;
# Line 123  void MeshAdapter::dump(const std::string Line 134  void MeshAdapter::dump(const std::string
134     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
135     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
136     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
137    #ifdef PASO_MPI
138       MPI_Status status;
139    #endif
140    
141     // I don't think the strdup is needed since Paso_MPI_appendRankToFileName  /* Incoming token indicates it's my turn to write */
142     // does it's own allocation.  #ifdef PASO_MPI
143     // char *newFileName = Paso_MPI_appendRankToFileName(strdup(fileName.c_str()), mpi_size, mpi_rank);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
144    #endif
145    
146     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),
147                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
148    
149     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
150     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
151       num_Tags = 0;
152     if (tag_map) {     if (tag_map) {
153        while (tag_map) {        while (tag_map) {
154           num_Tags++;           num_Tags++;
# Line 280  void MeshAdapter::dump(const std::string Line 296  void MeshAdapter::dump(const std::string
296        if (! (ids->put(int_ptr, mpi_size+1)) )        if (! (ids->put(int_ptr, mpi_size+1)) )
297           throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
298    
299          // Nodes nodeDistribution
300          if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
301             throw DataException("Error - MeshAdapter::dump: appending Nodes_NodeDistribution to netCDF file failed: " + *newFileName);
302          int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
303          if (! (ids->put(int_ptr, mpi_size+1)) )
304             throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);
305    
306     }     }
307    
308     // // // // // Elements // // // // //     // // // // // Elements // // // // //
309    
310     if (num_Elements>0) {     if (num_Elements>0) {
311    
       // Temp storage to gather node IDs  
       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);  
   
312        // Elements_Id        // Elements_Id
313        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
314           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);
# Line 318  void MeshAdapter::dump(const std::string Line 338  void MeshAdapter::dump(const std::string
338           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
339    
340        // Elements_Nodes        // Elements_Nodes
       for (int i=0; i<num_Elements; i++)  
          for (int j=0; j<num_Elements_numNodes; j++)  
             Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];  
341        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
342           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
343        if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
344           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
345    
       TMPMEMFREE(Elements_Nodes);  
   
346     }     }
347    
348     // // // // // Face_Elements // // // // //     // // // // // Face_Elements // // // // //
349    
350     if (num_FaceElements>0) {     if (num_FaceElements>0) {
351    
       // Temp storage to gather node IDs  
       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);  
   
352        // FaceElements_Id        // FaceElements_Id
353        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
354           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);
# Line 366  void MeshAdapter::dump(const std::string Line 378  void MeshAdapter::dump(const std::string
378           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
379    
380        // FaceElements_Nodes        // FaceElements_Nodes
       for (int i=0; i<num_FaceElements; i++)  
          for (int j=0; j<num_FaceElements_numNodes; j++)  
             FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];  
381        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
382           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
383        if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
384           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
385    
       TMPMEMFREE(FaceElements_Nodes);  
   
386     }     }
387    
388     // // // // // Contact_Elements // // // // //     // // // // // Contact_Elements // // // // //
389    
390     if (num_ContactElements>0) {     if (num_ContactElements>0) {
391    
       // Temp storage to gather node IDs  
       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);  
   
392        // ContactElements_Id        // ContactElements_Id
393        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
394           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);
# Line 414  void MeshAdapter::dump(const std::string Line 418  void MeshAdapter::dump(const std::string
418           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
419    
420        // ContactElements_Nodes        // ContactElements_Nodes
       for (int i=0; i<num_ContactElements; i++)  
          for (int j=0; j<num_ContactElements_numNodes; j++)  
             ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];  
421        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
422           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
423        if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
424           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
425    
       TMPMEMFREE(ContactElements_Nodes);  
   
426     }     }
427    
428     // // // // // Points // // // // //     // // // // // Points // // // // //
# Line 432  void MeshAdapter::dump(const std::string Line 431  void MeshAdapter::dump(const std::string
431    
432        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
433    
       // Temp storage to gather node IDs  
       int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
   
434        // Points_Id        // Points_Id
435        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
436           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
# Line 465  void MeshAdapter::dump(const std::string Line 461  void MeshAdapter::dump(const std::string
461    
462        // Points_Nodes        // Points_Nodes
463        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
       for (int i=0; i<num_Points; i++)  
          Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];  
464        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
465           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
466        if (! (ids->put(&(Points_Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
467           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
468    
       TMPMEMFREE(Points_Nodes);  
   
469     }     }
470    
471     // // // // // TagMap // // // // //     // // // // // TagMap // // // // //
# Line 520  void MeshAdapter::dump(const std::string Line 512  void MeshAdapter::dump(const std::string
512    
513     }     }
514    
515    /* Send token to next MPI process so he can take his turn */
516    #ifdef PASO_MPI
517       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
518    #endif
519    
520     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
521    
522  #else  #else
523     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
524  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
# Line 646  int MeshAdapter::getDiracDeltaFunctionCo Line 643  int MeshAdapter::getDiracDeltaFunctionCo
643  //  //
644  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
645  {  {
646     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
647       int numDim=Finley_Mesh_getDim(mesh);
648     checkFinleyError();     checkFinleyError();
649     return numDim;     return numDim;
650  }  }
651    
652  //  //
653    // Return the number of data points summed across all MPI processes
654    //
655    int MeshAdapter::getNumDataPointsGlobal() const
656    {
657       return Finley_NodeFile_getGlobalNumNodes(m_finleyMesh.get()->Nodes);
658    }
659    
660    //
661  // return the number of data points per sample and the number of samples  // return the number of data points per sample and the number of samples
662  // needed to represent data on a parts of the mesh.  // needed to represent data on a parts of the mesh.
663  //  //
# Line 827  void MeshAdapter::addPDEToTransportProbl Line 833  void MeshAdapter::addPDEToTransportProbl
833                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
834                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
835  {  {
836     DataArrayView::ShapeType shape;     DataTypes::ShapeType shape;
837     source.expand();     source.expand();
838     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
839     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
# Line 863  void MeshAdapter::addPDEToTransportProbl Line 869  void MeshAdapter::addPDEToTransportProbl
869  //  //
870  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
871  {  {
872     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
873     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
874     if (inDomain!=*this)       if (inDomain!=*this)  
875        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
876     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 875  void MeshAdapter::interpolateOnDomain(es Line 881  void MeshAdapter::interpolateOnDomain(es
881     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
882     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
883     case(Nodes):     case(Nodes):
884     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
885     case(Nodes):        case(Nodes):
886     case(ReducedNodes):        case(ReducedNodes):
887     case(DegreesOfFreedom):        case(DegreesOfFreedom):
888     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
889     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
    break;  
    case(Elements):  
    case(ReducedElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);  
    break;  
    case(FaceElements):  
    case(ReducedFaceElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);  
    break;  
    case(Points):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);  
    break;  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
890        break;        break;
891     }        case(Elements):
892     break;        case(ReducedElements):
893     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
894     switch(target.getFunctionSpace().getTypeCode()) {        break;
895     case(Nodes):        case(FaceElements):
896     case(ReducedNodes):        case(ReducedFaceElements):
897     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
898     case(ReducedDegreesOfFreedom):        break;
899     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
900     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
901     case(Elements):        break;
902     case(ReducedElements):        case(ContactElementsZero):
903     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
904     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
905     case(FaceElements):        break;
906     case(ReducedFaceElements):        case(ContactElementsOne):
907     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
908     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
909     case(Points):        break;
910     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
911     break;           stringstream temp;
912     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
913     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
914     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
915     break;        }
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
916        break;        break;
    }  
    break;  
    case(Elements):  
    if (target.getFunctionSpace().getTypeCode()==Elements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements possible.");  
    }  
    break;  
    case(ReducedElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");  
    }  
    break;  
    case(FaceElements):  
    if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");  
    }  
    break;  
    case(ReducedFaceElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");  
    }  
    break;  
    case(Points):  
    if (target.getFunctionSpace().getTypeCode()==Points) {  
       Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
    }  
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
    }  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");  
    }  
    break;  
    case(DegreesOfFreedom):        
    switch(target.getFunctionSpace().getTypeCode()) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
   
    case(Nodes):  
917     case(ReducedNodes):     case(ReducedNodes):
918     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
919        escript::Data temp=escript::Data(in);        case(Nodes):
920        temp.expand();        case(ReducedNodes):
921        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
922        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
923        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
924     }        break;
925     break;        case(Elements):
926     case(Elements):        case(ReducedElements):
    case(ReducedElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
    } else {  
927        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
928     }        break;
929     break;        case(FaceElements):
930     case(FaceElements):        case(ReducedFaceElements):
    case(ReducedFaceElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
   
    } else {  
931        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
932     }        break;
933     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
934        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
935     }        break;
936     break;        case(ContactElementsZero):
937     case(ContactElementsZero):        case(ReducedContactElementsZero):
    case(ContactElementsOne):  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
    } else {  
938        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
939        break;        break;
940     }        case(ContactElementsOne):
941     break;        case(ReducedContactElementsOne):
942     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
943     switch(target.getFunctionSpace().getTypeCode()) {        break;
944     case(Nodes):        default:
945     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
946     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
947     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
948     if (getMPISize()>1) {           break;
949        escript::Data temp=escript::Data(in);        }
950        temp.expand();        break;
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
    } else {  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    }  
    break;  
    case(DegreesOfFreedom):  
    throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
    break;  
    case(ReducedDegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
951     case(Elements):     case(Elements):
952          if (target.getFunctionSpace().getTypeCode()==Elements) {
953             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
954          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
955             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
956          } else {
957             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
958          }
959          break;
960     case(ReducedElements):     case(ReducedElements):
961     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
962        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
963        escriptDataC _in2 = temp.getDataC();        } else {
964        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
965     } else {        }
966        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
967     case(FaceElements):     case(FaceElements):
968          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
969             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
970          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
971             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
972          } else {
973             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
974          }
975          break;
976     case(ReducedFaceElements):     case(ReducedFaceElements):
977     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
978        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
979        escriptDataC _in2 = temp.getDataC();        } else {
980        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
981     } else {        }
982        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
983     case(Points):     case(Points):
984     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
985        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
986        escriptDataC _in2 = temp.getDataC();        } else {
987        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
988     } else {        }
989        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
990     case(ContactElementsZero):     case(ContactElementsZero):
991     case(ContactElementsOne):     case(ContactElementsOne):
992          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
993             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
994          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
995             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
996          } else {
997             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
998          }
999          break;
1000     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1001     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1002     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1003        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1004        escriptDataC _in2 = temp.getDataC();        } else {
1005        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1006     } else {        }
1007        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1008     }     case(DegreesOfFreedom):      
1009     break;        switch(target.getFunctionSpace().getTypeCode()) {
1010     default:        case(ReducedDegreesOfFreedom):
1011        stringstream temp;        case(DegreesOfFreedom):
1012        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1013        throw FinleyAdapterException(temp.str());        break;
1014      
1015          case(Nodes):
1016          case(ReducedNodes):
1017          if (getMPISize()>1) {
1018             escript::Data temp=escript::Data(in);
1019             temp.expand();
1020             escriptDataC _in2 = temp.getDataC();
1021             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1022          } else {
1023             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1024          }
1025          break;
1026          case(Elements):
1027          case(ReducedElements):
1028          if (getMPISize()>1) {
1029             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1030             escriptDataC _in2 = temp.getDataC();
1031             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1032          } else {
1033             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1034          }
1035          break;
1036          case(FaceElements):
1037          case(ReducedFaceElements):
1038          if (getMPISize()>1) {
1039             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1040             escriptDataC _in2 = temp.getDataC();
1041             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1042      
1043          } else {
1044             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1045          }
1046          break;
1047          case(Points):
1048          if (getMPISize()>1) {
1049             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1050             escriptDataC _in2 = temp.getDataC();
1051          } else {
1052             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1053          }
1054          break;
1055          case(ContactElementsZero):
1056          case(ContactElementsOne):
1057          case(ReducedContactElementsZero):
1058          case(ReducedContactElementsOne):
1059          if (getMPISize()>1) {
1060             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1061             escriptDataC _in2 = temp.getDataC();
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1063          } else {
1064             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1065          }
1066          break;
1067          default:
1068             stringstream temp;
1069             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1070             throw FinleyAdapterException(temp.str());
1071             break;
1072          }
1073          break;
1074       case(ReducedDegreesOfFreedom):
1075          switch(target.getFunctionSpace().getTypeCode()) {
1076          case(Nodes):
1077          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1078          break;
1079          case(ReducedNodes):
1080          if (getMPISize()>1) {
1081             escript::Data temp=escript::Data(in);
1082             temp.expand();
1083             escriptDataC _in2 = temp.getDataC();
1084             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1085          } else {
1086             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1087          }
1088          break;
1089          case(DegreesOfFreedom):
1090          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1091          break;
1092          case(ReducedDegreesOfFreedom):
1093          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1094          break;
1095          case(Elements):
1096          case(ReducedElements):
1097          if (getMPISize()>1) {
1098             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1099             escriptDataC _in2 = temp.getDataC();
1100             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1101          } else {
1102             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1103          }
1104          break;
1105          case(FaceElements):
1106          case(ReducedFaceElements):
1107          if (getMPISize()>1) {
1108             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1109             escriptDataC _in2 = temp.getDataC();
1110             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1111          } else {
1112             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1113          }
1114          break;
1115          case(Points):
1116          if (getMPISize()>1) {
1117             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1118             escriptDataC _in2 = temp.getDataC();
1119             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1120          } else {
1121             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1122          }
1123          break;
1124          case(ContactElementsZero):
1125          case(ContactElementsOne):
1126          case(ReducedContactElementsZero):
1127          case(ReducedContactElementsOne):
1128          if (getMPISize()>1) {
1129             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1130             escriptDataC _in2 = temp.getDataC();
1131             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1132          } else {
1133             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1134          }
1135          break;
1136          default:
1137             stringstream temp;
1138             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1139             throw FinleyAdapterException(temp.str());
1140             break;
1141          }
1142        break;        break;
    }  
    break;  
1143     default:     default:
1144        stringstream temp;        stringstream temp;
1145        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
# Line 1148  void MeshAdapter::interpolateOnDomain(es Line 1154  void MeshAdapter::interpolateOnDomain(es
1154  //  //
1155  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1156  {  {
1157     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1158     if (argDomain!=*this)     if (argDomain!=*this)
1159        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1160     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1171  void MeshAdapter::setToX(escript::Data& Line 1177  void MeshAdapter::setToX(escript::Data&
1177  //  //
1178  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1179  {  {
1180     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1181       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1182     if (normalDomain!=*this)     if (normalDomain!=*this)
1183        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1184     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1226  void MeshAdapter::setToNormal(escript::D Line 1233  void MeshAdapter::setToNormal(escript::D
1233  //  //
1234  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1235  {  {
1236     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1237     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1238       if (targetDomain!=this)
1239        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1240    
1241     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1238  void MeshAdapter::interpolateACross(escr Line 1246  void MeshAdapter::interpolateACross(escr
1246  //  //
1247  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1248  {  {
1249     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1250     if (argDomain!=*this)     if (argDomain!=*this)
1251        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1252    
# Line 1249  void MeshAdapter::setToIntegrals(std::ve Line 1257  void MeshAdapter::setToIntegrals(std::ve
1257     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1258     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1259     case(Nodes):     case(Nodes):
1260     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1261     _temp=temp.getDataC();     _temp=temp.getDataC();
1262     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1263     break;     break;
1264     case(ReducedNodes):     case(ReducedNodes):
1265     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1266     _temp=temp.getDataC();     _temp=temp.getDataC();
1267     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1268     break;     break;
# Line 1286  void MeshAdapter::setToIntegrals(std::ve Line 1294  void MeshAdapter::setToIntegrals(std::ve
1294     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1295     break;     break;
1296     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1297     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1298     _temp=temp.getDataC();     _temp=temp.getDataC();
1299     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1300     break;     break;
1301     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1302     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1303     _temp=temp.getDataC();     _temp=temp.getDataC();
1304     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1305     break;     break;
# Line 1310  void MeshAdapter::setToIntegrals(std::ve Line 1318  void MeshAdapter::setToIntegrals(std::ve
1318  //  //
1319  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1320  {  {
1321     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1322     if (argDomain!=*this)     if (argDomain!=*this)
1323        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1324     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1325     if (gradDomain!=*this)     if (gradDomain!=*this)
1326        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1327    
# Line 1435  void MeshAdapter::setToSize(escript::Dat Line 1443  void MeshAdapter::setToSize(escript::Dat
1443     checkFinleyError();     checkFinleyError();
1444  }  }
1445    
1446  // sets the location of nodes:  //
1447    // sets the location of nodes
1448    //
1449  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1450  {  {
1451     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1452     escriptDataC tmp;     escriptDataC tmp;
1453     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1454     if (newDomain!=*this)     if (newDomain!=*this)
1455        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1456     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1457     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1458           Finley_Mesh_setCoordinates(mesh,&tmp);
1459       } else {
1460           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );
1461           tmp = new_x_inter.getDataC();
1462           Finley_Mesh_setCoordinates(mesh,&tmp);
1463       }
1464     checkFinleyError();     checkFinleyError();
1465  }  }
1466    
1467  // saves a data array in openDX format:  //
1468  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1469    // corresponding pointers from python dictionary. Caller must free arrays.
1470    //
1471    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1472  {  {
1473     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1474     /* win32 refactor */     /* win32 refactor */
1475     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1476     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1477     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1478    
1479     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1480     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1481        std::string n=boost::python::extract<std::string>(keys[i]);        std::string n=boost::python::extract<std::string>(keys[i]);
1482        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1483        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1484           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1485        data[i]=d.getDataC();        data[i] = d.getDataC();
1486        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1487        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1488           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1489     }     }
1490     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1491    
1492    //
1493    // saves mesh and optionally data arrays in openDX format
1494    //
1495    void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1496    {
1497       int num_data;
1498       char **names;
1499       escriptDataC *data;
1500       escriptDataC **ptr_data;
1501    
1502       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1503       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1504     checkFinleyError();     checkFinleyError();
1505    
1506     /* win32 refactor */     /* win32 refactor */
1507     TMPMEMFREE(data);     TMPMEMFREE(data);
1508     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1509     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1510     {     {
1511        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1512     }     }
# Line 1493  void MeshAdapter::saveDX(const std::stri Line 1515  void MeshAdapter::saveDX(const std::stri
1515     return;     return;
1516  }  }
1517    
1518  // saves a data array in openVTK format:  //
1519  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1520    //
1521    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const
1522  {  {
1523     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     int num_data;
1524     const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1525     /* win32 refactor */     escriptDataC *data;
1526     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     escriptDataC **ptr_data;
    for(int i=0;i<num_data;i++)  
    {  
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
          throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]=&(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1527    
1528       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1529       Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1530     checkFinleyError();     checkFinleyError();
1531    
1532     /* win32 refactor */     /* win32 refactor */
1533     TMPMEMFREE(data);     TMPMEMFREE(data);
1534     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1535     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1536     {     {
1537        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1538     }     }
1539     TMPMEMFREE(names);     TMPMEMFREE(names);
   
    return;  
1540  }  }
1541                                                                                                                                                                        
1542                                                                                                                                                                        //
1543  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1544    //
1545  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1546                                                   const int row_blocksize,                                                   const int row_blocksize,
1547                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
# Line 1550  SystemMatrixAdapter MeshAdapter::newSyst Line 1552  SystemMatrixAdapter MeshAdapter::newSyst
1552     int reduceRowOrder=0;     int reduceRowOrder=0;
1553     int reduceColOrder=0;     int reduceColOrder=0;
1554     // is the domain right?     // is the domain right?
1555     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1556     if (row_domain!=*this)     if (row_domain!=*this)
1557        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1558     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1559     if (col_domain!=*this)     if (col_domain!=*this)
1560        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1561     // is the function space type right     // is the function space type right
# Line 1589  SystemMatrixAdapter MeshAdapter::newSyst Line 1591  SystemMatrixAdapter MeshAdapter::newSyst
1591     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1592     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1593  }  }
1594    
1595    //
1596  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1597    //
1598  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1599                                                           const double theta,                                                           const double theta,
1600                                                           const int blocksize,                                                           const int blocksize,
# Line 1598  TransportProblemAdapter MeshAdapter::new Line 1603  TransportProblemAdapter MeshAdapter::new
1603  {  {
1604     int reduceOrder=0;     int reduceOrder=0;
1605     // is the domain right?     // is the domain right?
1606     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1607     if (domain!=*this)     if (domain!=*this)
1608        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1609     // is the function space type right     // is the function space type right
# Line 1829  bool MeshAdapter::operator!=(const Abstr Line 1834  bool MeshAdapter::operator!=(const Abstr
1834     return !(operator==(other));     return !(operator==(other));
1835  }  }
1836    
1837  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1838    {
1839       Finley_Mesh* mesh=m_finleyMesh.get();
1840       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1841       checkPasoError();
1842       return out;
1843    }
1844    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1845  {  {
1846     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
1847       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1848     checkPasoError();     checkPasoError();
1849     return out;     return out;
1850  }  }
# Line 1848  escript::Data MeshAdapter::getNormal() c Line 1861  escript::Data MeshAdapter::getNormal() c
1861    
1862  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1863  {  {
1864     return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
1865  }  }
1866    
1867  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1868  {  {
1869     int *out = NULL;     int *out = NULL;
1870     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2047  std::string MeshAdapter::showTagNames() Line 2060  std::string MeshAdapter::showTagNames()
2060     return temp.str();     return temp.str();
2061  }  }
2062    
2063    int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2064    {
2065      Finley_Mesh* mesh=m_finleyMesh.get();
2066      dim_t numTags=0;
2067      switch(functionSpaceCode) {
2068       case(Nodes):
2069              numTags=mesh->Nodes->numTagsInUse;
2070              break;
2071       case(ReducedNodes):
2072              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2073              break;
2074       case(DegreesOfFreedom):
2075              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2076              break;
2077       case(ReducedDegreesOfFreedom):
2078              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2079              break;
2080       case(Elements):
2081       case(ReducedElements):
2082              numTags=mesh->Elements->numTagsInUse;
2083              break;
2084       case(FaceElements):
2085       case(ReducedFaceElements):
2086              numTags=mesh->FaceElements->numTagsInUse;
2087              break;
2088       case(Points):
2089              numTags=mesh->Points->numTagsInUse;
2090              break;
2091       case(ContactElementsZero):
2092       case(ReducedContactElementsZero):
2093       case(ContactElementsOne):
2094       case(ReducedContactElementsOne):
2095              numTags=mesh->ContactElements->numTagsInUse;
2096              break;
2097       default:
2098          stringstream temp;
2099          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2100          throw FinleyAdapterException(temp.str());
2101      }
2102      return numTags;
2103    }
2104    
2105    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2106    {
2107      Finley_Mesh* mesh=m_finleyMesh.get();
2108      index_t* tags=NULL;
2109      switch(functionSpaceCode) {
2110       case(Nodes):
2111              tags=mesh->Nodes->tagsInUse;
2112              break;
2113       case(ReducedNodes):
2114              throw FinleyAdapterException("Error - ReducedNodes does not support tags");
2115              break;
2116       case(DegreesOfFreedom):
2117              throw FinleyAdapterException("Error - DegreesOfFreedom does not support tags");
2118              break;
2119       case(ReducedDegreesOfFreedom):
2120              throw FinleyAdapterException("Error - ReducedDegreesOfFreedom does not support tags");
2121              break;
2122       case(Elements):
2123       case(ReducedElements):
2124              tags=mesh->Elements->tagsInUse;
2125              break;
2126       case(FaceElements):
2127       case(ReducedFaceElements):
2128              tags=mesh->FaceElements->tagsInUse;
2129              break;
2130       case(Points):
2131              tags=mesh->Points->tagsInUse;
2132              break;
2133       case(ContactElementsZero):
2134       case(ReducedContactElementsZero):
2135       case(ContactElementsOne):
2136       case(ReducedContactElementsOne):
2137              tags=mesh->ContactElements->tagsInUse;
2138              break;
2139       default:
2140          stringstream temp;
2141          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2142          throw FinleyAdapterException(temp.str());
2143      }
2144      return tags;
2145    }
2146    
2147    
2148    bool MeshAdapter::canTag(int functionSpaceCode) const
2149    {
2150      switch(functionSpaceCode) {
2151       case(Nodes):
2152       case(Elements):
2153       case(ReducedElements):
2154       case(FaceElements):
2155       case(ReducedFaceElements):
2156       case(Points):
2157       case(ContactElementsZero):
2158       case(ReducedContactElementsZero):
2159       case(ContactElementsOne):
2160       case(ReducedContactElementsOne):
2161              return true;
2162       case(ReducedNodes):
2163       case(DegreesOfFreedom):
2164       case(ReducedDegreesOfFreedom):
2165          return false;
2166       default:
2167        return false;
2168      }
2169    }
2170    
2171    AbstractDomain::StatusType MeshAdapter::getStatus() const
2172    {
2173      Finley_Mesh* mesh=m_finleyMesh.get();
2174      return Finley_Mesh_getStatus(mesh);
2175    }
2176    
2177    
2178    
2179  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1628  
changed lines
  Added in v.2548

  ViewVC Help
Powered by ViewVC 1.1.26