/[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 1801 by ksteube, Fri Sep 19 01:37:09 2008 UTC revision 2640 by jfenwick, Mon Aug 31 06:22:10 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 24  Line 23 
23  #include "paso/Paso_MPI.h"  #include "paso/Paso_MPI.h"
24  #endif  #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 87  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    #ifdef PASO_MPI
100      MPI_Comm
101    #else
102      unsigned int
103    #endif
104    MeshAdapter::getMPIComm() const
105    {
106    #ifdef PASO_MPI
107        return mesh->MPIInfo->comm;
108    #else
109        return 0;
110    #endif
111    }
112    
113    
114  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
# Line 102  void MeshAdapter::write(const std::strin Line 124  void MeshAdapter::write(const std::strin
124     TMPMEMFREE(fName);     TMPMEMFREE(fName);
125  }  }
126    
127  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
128  {  {
129     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
130  }  }
# Line 110  void MeshAdapter::Print_Mesh_Info(const Line 132  void MeshAdapter::Print_Mesh_Info(const
132  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const std::string& fileName) const
133  {  {
134  #ifdef USE_NETCDF  #ifdef USE_NETCDF
135     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
136     NcVar *ids;     NcVar *ids;
137     int *int_ptr;     int *int_ptr;
138     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
# Line 636  int MeshAdapter::getDiracDeltaFunctionCo Line 658  int MeshAdapter::getDiracDeltaFunctionCo
658  //  //
659  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
660  {  {
661     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
662       int numDim=Finley_Mesh_getDim(mesh);
663     checkFinleyError();     checkFinleyError();
664     return numDim;     return numDim;
665  }  }
# Line 861  void MeshAdapter::addPDEToTransportProbl Line 884  void MeshAdapter::addPDEToTransportProbl
884  //  //
885  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
886  {  {
887     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
888     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
889     if (inDomain!=*this)       if (inDomain!=*this)  
890        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
891     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 873  void MeshAdapter::interpolateOnDomain(es Line 896  void MeshAdapter::interpolateOnDomain(es
896     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
897     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
898     case(Nodes):     case(Nodes):
899     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
900     case(Nodes):        case(Nodes):
901     case(ReducedNodes):        case(ReducedNodes):
902     case(DegreesOfFreedom):        case(DegreesOfFreedom):
903     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
904     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());  
905        break;        break;
906     }        case(Elements):
907     break;        case(ReducedElements):
908     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
909     switch(target.getFunctionSpace().getTypeCode()) {        break;
910     case(Nodes):        case(FaceElements):
911     case(ReducedNodes):        case(ReducedFaceElements):
912     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
913     case(ReducedDegreesOfFreedom):        break;
914     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
915     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
916     case(Elements):        break;
917     case(ReducedElements):        case(ContactElementsZero):
918     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
919     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
920     case(FaceElements):        break;
921     case(ReducedFaceElements):        case(ContactElementsOne):
922     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
923     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
924     case(Points):        break;
925     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
926     break;           stringstream temp;
927     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
928     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
929     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
930     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());  
931        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):  
932     case(ReducedNodes):     case(ReducedNodes):
933     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
934        escript::Data temp=escript::Data(in);        case(Nodes):
935        temp.expand();        case(ReducedNodes):
936        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
937        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
938        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
939     }        break;
940     break;        case(Elements):
941     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 {  
942        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
943     }        break;
944     break;        case(FaceElements):
945     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 {  
946        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
947     }        break;
948     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
949        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
950     }        break;
951     break;        case(ContactElementsZero):
952     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 {  
953        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());  
954        break;        break;
955     }        case(ContactElementsOne):
956     break;        case(ReducedContactElementsOne):
957     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
958     switch(target.getFunctionSpace().getTypeCode()) {        break;
959     case(Nodes):        default:
960     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
961     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
962     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
963     if (getMPISize()>1) {           break;
964        escript::Data temp=escript::Data(in);        }
965        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;  
966     case(Elements):     case(Elements):
967          if (target.getFunctionSpace().getTypeCode()==Elements) {
968             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
969          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
970             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
971          } else {
972             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
973          }
974          break;
975     case(ReducedElements):     case(ReducedElements):
976     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
977        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
978        escriptDataC _in2 = temp.getDataC();        } else {
979        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
980     } else {        }
981        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
982     case(FaceElements):     case(FaceElements):
983          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
984             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
985          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
986             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
987          } else {
988             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
989          }
990          break;
991     case(ReducedFaceElements):     case(ReducedFaceElements):
992     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
993        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
994        escriptDataC _in2 = temp.getDataC();        } else {
995        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
996     } else {        }
997        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
998     case(Points):     case(Points):
999     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
1000        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
1001        escriptDataC _in2 = temp.getDataC();        } else {
1002        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
1003     } else {        }
1004        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
1005     case(ContactElementsZero):     case(ContactElementsZero):
1006     case(ContactElementsOne):     case(ContactElementsOne):
1007          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
1008             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1009          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1010             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
1011          } else {
1012             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
1013          }
1014          break;
1015     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1016     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1017     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1018        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1019        escriptDataC _in2 = temp.getDataC();        } else {
1020        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1021     } else {        }
1022        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1023     }     case(DegreesOfFreedom):      
1024     break;        switch(target.getFunctionSpace().getTypeCode()) {
1025     default:        case(ReducedDegreesOfFreedom):
1026        stringstream temp;        case(DegreesOfFreedom):
1027        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1028        throw FinleyAdapterException(temp.str());        break;
1029      
1030          case(Nodes):
1031          case(ReducedNodes):
1032          if (getMPISize()>1) {
1033             escript::Data temp=escript::Data(in);
1034             temp.expand();
1035             escriptDataC _in2 = temp.getDataC();
1036             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1037          } else {
1038             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1039          }
1040          break;
1041          case(Elements):
1042          case(ReducedElements):
1043          if (getMPISize()>1) {
1044             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1045             escriptDataC _in2 = temp.getDataC();
1046             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1047          } else {
1048             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1049          }
1050          break;
1051          case(FaceElements):
1052          case(ReducedFaceElements):
1053          if (getMPISize()>1) {
1054             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1055             escriptDataC _in2 = temp.getDataC();
1056             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1057      
1058          } else {
1059             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1060          }
1061          break;
1062          case(Points):
1063          if (getMPISize()>1) {
1064             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1065             escriptDataC _in2 = temp.getDataC();
1066          } else {
1067             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1068          }
1069          break;
1070          case(ContactElementsZero):
1071          case(ContactElementsOne):
1072          case(ReducedContactElementsZero):
1073          case(ReducedContactElementsOne):
1074          if (getMPISize()>1) {
1075             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1076             escriptDataC _in2 = temp.getDataC();
1077             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1078          } else {
1079             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1080          }
1081          break;
1082          default:
1083             stringstream temp;
1084             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1085             throw FinleyAdapterException(temp.str());
1086             break;
1087          }
1088          break;
1089       case(ReducedDegreesOfFreedom):
1090          switch(target.getFunctionSpace().getTypeCode()) {
1091          case(Nodes):
1092          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1093          break;
1094          case(ReducedNodes):
1095          if (getMPISize()>1) {
1096             escript::Data temp=escript::Data(in);
1097             temp.expand();
1098             escriptDataC _in2 = temp.getDataC();
1099             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1100          } else {
1101             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1102          }
1103          break;
1104          case(DegreesOfFreedom):
1105          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1106          break;
1107          case(ReducedDegreesOfFreedom):
1108          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1109          break;
1110          case(Elements):
1111          case(ReducedElements):
1112          if (getMPISize()>1) {
1113             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1114             escriptDataC _in2 = temp.getDataC();
1115             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1116          } else {
1117             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1118          }
1119          break;
1120          case(FaceElements):
1121          case(ReducedFaceElements):
1122          if (getMPISize()>1) {
1123             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1124             escriptDataC _in2 = temp.getDataC();
1125             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1126          } else {
1127             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1128          }
1129          break;
1130          case(Points):
1131          if (getMPISize()>1) {
1132             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1133             escriptDataC _in2 = temp.getDataC();
1134             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1135          } else {
1136             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1137          }
1138          break;
1139          case(ContactElementsZero):
1140          case(ContactElementsOne):
1141          case(ReducedContactElementsZero):
1142          case(ReducedContactElementsOne):
1143          if (getMPISize()>1) {
1144             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1145             escriptDataC _in2 = temp.getDataC();
1146             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1147          } else {
1148             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1149          }
1150          break;
1151          default:
1152             stringstream temp;
1153             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1154             throw FinleyAdapterException(temp.str());
1155             break;
1156          }
1157        break;        break;
    }  
    break;  
1158     default:     default:
1159        stringstream temp;        stringstream temp;
1160        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 1146  void MeshAdapter::interpolateOnDomain(es Line 1169  void MeshAdapter::interpolateOnDomain(es
1169  //  //
1170  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1171  {  {
1172     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1173     if (argDomain!=*this)     if (argDomain!=*this)
1174        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1175     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1169  void MeshAdapter::setToX(escript::Data& Line 1192  void MeshAdapter::setToX(escript::Data&
1192  //  //
1193  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1194  {  {
1195     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1196       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1197     if (normalDomain!=*this)     if (normalDomain!=*this)
1198        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1199     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1224  void MeshAdapter::setToNormal(escript::D Line 1248  void MeshAdapter::setToNormal(escript::D
1248  //  //
1249  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1250  {  {
1251     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1252     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1253       if (targetDomain!=this)
1254        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1255    
1256     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1236  void MeshAdapter::interpolateACross(escr Line 1261  void MeshAdapter::interpolateACross(escr
1261  //  //
1262  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1263  {  {
1264     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1265     if (argDomain!=*this)     if (argDomain!=*this)
1266        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1267    
# Line 1247  void MeshAdapter::setToIntegrals(std::ve Line 1272  void MeshAdapter::setToIntegrals(std::ve
1272     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1273     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1274     case(Nodes):     case(Nodes):
1275     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1276     _temp=temp.getDataC();     _temp=temp.getDataC();
1277     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1278     break;     break;
1279     case(ReducedNodes):     case(ReducedNodes):
1280     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1281     _temp=temp.getDataC();     _temp=temp.getDataC();
1282     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1283     break;     break;
# Line 1284  void MeshAdapter::setToIntegrals(std::ve Line 1309  void MeshAdapter::setToIntegrals(std::ve
1309     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1310     break;     break;
1311     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1312     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1313     _temp=temp.getDataC();     _temp=temp.getDataC();
1314     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1315     break;     break;
1316     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1317     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1318     _temp=temp.getDataC();     _temp=temp.getDataC();
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320     break;     break;
# Line 1308  void MeshAdapter::setToIntegrals(std::ve Line 1333  void MeshAdapter::setToIntegrals(std::ve
1333  //  //
1334  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1335  {  {
1336     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1337     if (argDomain!=*this)     if (argDomain!=*this)
1338        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1339     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1340     if (gradDomain!=*this)     if (gradDomain!=*this)
1341        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1342    
# Line 1433  void MeshAdapter::setToSize(escript::Dat Line 1458  void MeshAdapter::setToSize(escript::Dat
1458     checkFinleyError();     checkFinleyError();
1459  }  }
1460    
1461  // sets the location of nodes:  //
1462    // sets the location of nodes
1463    //
1464  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1465  {  {
1466     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1467     escriptDataC tmp;     escriptDataC tmp;
1468     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1469     if (newDomain!=*this)     if (newDomain!=*this)
1470        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1471     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1472     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1473           Finley_Mesh_setCoordinates(mesh,&tmp);
1474       } else {
1475           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );
1476           tmp = new_x_inter.getDataC();
1477           Finley_Mesh_setCoordinates(mesh,&tmp);
1478       }
1479     checkFinleyError();     checkFinleyError();
1480  }  }
1481    
1482  // saves a data array in openDX format:  //
1483  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
1484    // corresponding pointers from python dictionary. Caller must free arrays.
1485    //
1486    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1487  {  {
1488     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__")());  
1489     /* win32 refactor */     /* win32 refactor */
1490     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1491     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1492     {     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;  
1493    
1494     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1495     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1496        std::string n=boost::python::extract<std::string>(keys[i]);        std::string n=boost::python::extract<std::string>(keys[i]);
1497        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1498        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1499           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1500        data[i]=d.getDataC();        data[i] = d.getDataC();
1501        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1502        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1503           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());  
       }  
1504     }     }
1505     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1506    
1507    //
1508    // saves mesh and optionally data arrays in openDX format
1509    //
1510    void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1511    {
1512       int num_data;
1513       char **names;
1514       escriptDataC *data;
1515       escriptDataC **ptr_data;
1516    
1517       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1518       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1519     checkFinleyError();     checkFinleyError();
1520    
1521     /* win32 refactor */     /* win32 refactor */
1522     TMPMEMFREE(data);     TMPMEMFREE(data);
1523     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1524     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1525     {     {
1526        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1527     }     }
# Line 1491  void MeshAdapter::saveDX(const std::stri Line 1530  void MeshAdapter::saveDX(const std::stri
1530     return;     return;
1531  }  }
1532    
1533  // saves a data array in openVTK format:  //
1534  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1535    //
1536    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const
1537  {  {
1538     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     int num_data;
1539     const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1540     /* win32 refactor */     escriptDataC *data;
1541     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);  
1542    
1543       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1544       Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1545     checkFinleyError();     checkFinleyError();
1546    
1547     /* win32 refactor */     /* win32 refactor */
1548     TMPMEMFREE(data);     TMPMEMFREE(data);
1549     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1550     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1551     {     {
1552        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1553     }     }
1554     TMPMEMFREE(names);     TMPMEMFREE(names);
   
    return;  
1555  }  }
1556                                                                                                                                                                        
1557                                                                                                                                                                        //
1558  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1559    //
1560  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1561                                                   const int row_blocksize,                                                   const int row_blocksize,
1562                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
# Line 1548  SystemMatrixAdapter MeshAdapter::newSyst Line 1567  SystemMatrixAdapter MeshAdapter::newSyst
1567     int reduceRowOrder=0;     int reduceRowOrder=0;
1568     int reduceColOrder=0;     int reduceColOrder=0;
1569     // is the domain right?     // is the domain right?
1570     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1571     if (row_domain!=*this)     if (row_domain!=*this)
1572        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.");
1573     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1574     if (col_domain!=*this)     if (col_domain!=*this)
1575        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.");
1576     // is the function space type right     // is the function space type right
# Line 1581  SystemMatrixAdapter MeshAdapter::newSyst Line 1600  SystemMatrixAdapter MeshAdapter::newSyst
1600  #endif  #endif
1601     }     }
1602     else {     else {
1603        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1604     }     }
1605     checkPasoError();     checkPasoError();
1606     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1607     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1608  }  }
1609    
1610    //
1611  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1612    //
1613  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1614                                                           const double theta,                                                           const double theta,
1615                                                           const int blocksize,                                                           const int blocksize,
# Line 1596  TransportProblemAdapter MeshAdapter::new Line 1618  TransportProblemAdapter MeshAdapter::new
1618  {  {
1619     int reduceOrder=0;     int reduceOrder=0;
1620     // is the domain right?     // is the domain right?
1621     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1622     if (domain!=*this)     if (domain!=*this)
1623        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.");
1624     // is the function space type right     // is the function space type right
# Line 1654  bool MeshAdapter::isCellOriented(int fun Line 1676  bool MeshAdapter::isCellOriented(int fun
1676     return false;     return false;
1677  }  }
1678    
1679    bool
1680    MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1681    {
1682       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1683        class 1: DOF <-> Nodes
1684        class 2: ReducedDOF <-> ReducedNodes
1685        class 3: Points
1686        class 4: Elements
1687        class 5: ReducedElements
1688        class 6: FaceElements
1689        class 7: ReducedFaceElements
1690        class 8: ContactElementZero <-> ContactElementOne
1691        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1692    
1693       There is also a set of lines. Interpolation is possible down a line but not between lines.
1694       class 1 and 2 belong to all lines so aren't considered.
1695        line 0: class 3
1696        line 1: class 4,5
1697        line 2: class 6,7
1698        line 3: class 8,9
1699    
1700       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1701       eg hasnodes is true if we have at least one instance of Nodes.
1702       */
1703        if (fs.size()==0)
1704        {
1705        return false;
1706        }
1707        std::vector<int> hasclass(10);
1708        std::vector<int> hasline(4);    
1709        bool hasnodes=false;
1710        bool hasrednodes=false;
1711        bool hascez=false;
1712        bool hasrcez=false;
1713        for (int i=0;i<fs.size();++i)
1714        {
1715        switch(fs[i])
1716        {
1717        case(Nodes):   hasnodes=true;   // no break is deliberate
1718        case(DegreesOfFreedom):
1719            hasclass[1]=1;
1720            break;
1721        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1722        case(ReducedDegreesOfFreedom):
1723            hasclass[2]=1;
1724            break;
1725        case(Points):
1726            hasline[0]=1;
1727            hasclass[3]=1;
1728            break;
1729        case(Elements):
1730            hasclass[4]=1;
1731            hasline[1]=1;
1732            break;
1733        case(ReducedElements):
1734            hasclass[5]=1;
1735            hasline[1]=1;
1736            break;
1737        case(FaceElements):
1738            hasclass[6]=1;
1739            hasline[2]=1;
1740            break;
1741        case(ReducedFaceElements):
1742            hasclass[7]=1;
1743            hasline[2]=1;
1744            break;
1745        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1746        case(ContactElementsOne):
1747            hasclass[8]=1;
1748            hasline[3]=1;
1749            break;
1750        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1751        case(ReducedContactElementsOne):
1752            hasclass[9]=1;
1753            hasline[3]=1;
1754            break;
1755        default:
1756            return false;
1757        }
1758        }
1759        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1760        // fail if we have more than one leaf group
1761    
1762        if (totlines>1)
1763        {
1764        return false;   // there are at least two branches we can't interpolate between
1765        }
1766        else if (totlines==1)
1767        {
1768        if (hasline[0]==1)      // we have points
1769        {
1770            resultcode=Points;
1771        }
1772        else if (hasline[1]==1)
1773        {
1774            if (hasclass[5]==1)
1775            {
1776            resultcode=ReducedElements;
1777            }
1778            else
1779            {
1780            resultcode=Elements;
1781            }
1782        }
1783        else if (hasline[2]==1)
1784        {
1785            if (hasclass[7]==ReducedFaceElements)
1786            {
1787            resultcode=ReducedFaceElements;
1788            }
1789            else
1790            {
1791            resultcode=FaceElements;
1792            }
1793        }
1794        else    // so we must be in line3
1795        {
1796            if (hasclass[9]==1)
1797            {
1798            // need something from class 9
1799            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1800            }
1801            else
1802            {
1803            // something from class 8
1804            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1805            }
1806        }
1807        }
1808        else    // totlines==0
1809        {
1810        if (hasclass[2]==1)
1811        {
1812            // something from class 2
1813            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1814        }
1815        else
1816        {   // something from class 1
1817            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1818        }
1819        }
1820        return true;
1821    }
1822    
1823  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1824  {  {
1825     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1826     case(Nodes):     case(Nodes):
1827     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1828     case(Nodes):      case(Nodes):
1829     case(ReducedNodes):      case(ReducedNodes):
1830     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1831     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1832     case(Elements):      case(Elements):
1833     case(ReducedElements):      case(ReducedElements):
1834     case(FaceElements):      case(FaceElements):
1835     case(ReducedFaceElements):      case(ReducedFaceElements):
1836     case(Points):      case(Points):
1837     case(ContactElementsZero):      case(ContactElementsZero):
1838     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1839     case(ContactElementsOne):      case(ContactElementsOne):
1840     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1841     return true;      return true;
1842     default:      default:
1843        stringstream temp;            stringstream temp;
1844        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1845        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
1846     }     }
1847     break;     break;
1848     case(ReducedNodes):     case(ReducedNodes):
1849     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1850     case(ReducedNodes):      case(ReducedNodes):
1851     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1852     case(Elements):      case(Elements):
1853     case(ReducedElements):      case(ReducedElements):
1854     case(FaceElements):      case(FaceElements):
1855     case(ReducedFaceElements):      case(ReducedFaceElements):
1856     case(Points):      case(Points):
1857     case(ContactElementsZero):      case(ContactElementsZero):
1858     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1859     case(ContactElementsOne):      case(ContactElementsOne):
1860     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1861     return true;      return true;
1862     case(Nodes):      case(Nodes):
1863     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1864     return false;      return false;
1865     default:      default:
1866        stringstream temp;          stringstream temp;
1867        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1868        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1869     }     }
1870     break;     break;
1871     case(Elements):     case(Elements):
1872     if (functionSpaceType_target==Elements) {      if (functionSpaceType_target==Elements) {
1873        return true;        return true;
1874     } else if (functionSpaceType_target==ReducedElements) {      } else if (functionSpaceType_target==ReducedElements) {
1875        return true;        return true;
1876     } else {          } else {
1877        return false;            return false;
1878     }          }
1879     case(ReducedElements):     case(ReducedElements):
1880     if (functionSpaceType_target==ReducedElements) {      if (functionSpaceType_target==ReducedElements) {
1881        return true;        return true;
1882     } else {      } else {
1883        return false;            return false;
1884     }      }
1885     case(FaceElements):     case(FaceElements):
1886     if (functionSpaceType_target==FaceElements) {      if (functionSpaceType_target==FaceElements) {
1887        return true;              return true;
1888     } else if (functionSpaceType_target==ReducedFaceElements) {      } else if (functionSpaceType_target==ReducedFaceElements) {
1889        return true;              return true;
1890     } else {      } else {
1891        return false;              return false;
1892     }      }
1893     case(ReducedFaceElements):     case(ReducedFaceElements):
1894     if (functionSpaceType_target==ReducedFaceElements) {      if (functionSpaceType_target==ReducedFaceElements) {
1895        return true;              return true;
1896     } else {      } else {
1897        return false;          return false;
1898     }      }
1899     case(Points):     case(Points):
1900     if (functionSpaceType_target==Points) {      if (functionSpaceType_target==Points) {
1901        return true;              return true;
1902     } else {      } else {
1903        return false;              return false;
1904     }      }
1905     case(ContactElementsZero):     case(ContactElementsZero):
1906     case(ContactElementsOne):     case(ContactElementsOne):
1907     if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {      if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1908        return true;              return true;
1909     } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1910        return true;              return true;
1911     } else {      } else {
1912        return false;              return false;
1913     }      }
1914     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1915     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1916     if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {      if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1917        return true;              return true;
1918     } else {      } else {
1919        return false;              return false;
1920     }      }
1921     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1922     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1923     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1924     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1925     case(Nodes):      case(Nodes):
1926     case(ReducedNodes):      case(ReducedNodes):
1927     case(Elements):      case(Elements):
1928     case(ReducedElements):      case(ReducedElements):
1929     case(Points):      case(Points):
1930     case(FaceElements):      case(FaceElements):
1931     case(ReducedFaceElements):      case(ReducedFaceElements):
1932     case(ContactElementsZero):      case(ContactElementsZero):
1933     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1934     case(ContactElementsOne):      case(ContactElementsOne):
1935     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1936     return true;      return true;
1937     default:      default:
1938        stringstream temp;          stringstream temp;
1939        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1940        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1941     }      }
1942     break;      break;
1943     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1944     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1945     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1946     case(ReducedNodes):      case(ReducedNodes):
1947     case(Elements):      case(Elements):
1948     case(ReducedElements):      case(ReducedElements):
1949     case(FaceElements):      case(FaceElements):
1950     case(ReducedFaceElements):      case(ReducedFaceElements):
1951     case(Points):      case(Points):
1952     case(ContactElementsZero):      case(ContactElementsZero):
1953     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1954     case(ContactElementsOne):      case(ContactElementsOne):
1955     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1956     return true;      return true;
1957     case(Nodes):      case(Nodes):
1958     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1959     return false;      return false;
1960     default:      default:
1961        stringstream temp;          stringstream temp;
1962        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;          temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1963        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1964     }      }
1965     break;      break;
1966     default:     default:
1967        stringstream temp;        stringstream temp;
1968        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_source;
# Line 1827  bool MeshAdapter::operator!=(const Abstr Line 1993  bool MeshAdapter::operator!=(const Abstr
1993     return !(operator==(other));     return !(operator==(other));
1994  }  }
1995    
1996  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
1997    {
1998       Finley_Mesh* mesh=m_finleyMesh.get();
1999       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2000       checkPasoError();
2001       return out;
2002    }
2003    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2004  {  {
2005     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2006       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2007     checkPasoError();     checkPasoError();
2008     return out;     return out;
2009  }  }
# Line 1846  escript::Data MeshAdapter::getNormal() c Line 2020  escript::Data MeshAdapter::getNormal() c
2020    
2021  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2022  {  {
2023     return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
2024  }  }
2025    
2026  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2027  {  {
2028     int *out = NULL;     int *out = NULL;
2029     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2086  int MeshAdapter::getNumberOfTagsInUse(in Line 2260  int MeshAdapter::getNumberOfTagsInUse(in
2260    }    }
2261    return numTags;    return numTags;
2262  }  }
2263  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2264    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2265  {  {
2266    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2267    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2129  int* MeshAdapter::borrowListOfTagsInUse( Line 2304  int* MeshAdapter::borrowListOfTagsInUse(
2304  }  }
2305    
2306    
2307    bool MeshAdapter::canTag(int functionSpaceCode) const
2308    {
2309      switch(functionSpaceCode) {
2310       case(Nodes):
2311       case(Elements):
2312       case(ReducedElements):
2313       case(FaceElements):
2314       case(ReducedFaceElements):
2315       case(Points):
2316       case(ContactElementsZero):
2317       case(ReducedContactElementsZero):
2318       case(ContactElementsOne):
2319       case(ReducedContactElementsOne):
2320              return true;
2321       case(ReducedNodes):
2322       case(DegreesOfFreedom):
2323       case(ReducedDegreesOfFreedom):
2324          return false;
2325       default:
2326        return false;
2327      }
2328    }
2329    
2330    AbstractDomain::StatusType MeshAdapter::getStatus() const
2331    {
2332      Finley_Mesh* mesh=m_finleyMesh.get();
2333      return Finley_Mesh_getStatus(mesh);
2334    }
2335    
2336    
2337    
2338  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1801  
changed lines
  Added in v.2640

  ViewVC Help
Powered by ViewVC 1.1.26