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

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

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

revision 1990 by ksteube, Fri Nov 7 04:19:07 2008 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2008 by University of Queensland  * Copyright (c) 2003-2009 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 23  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 99  bool MeshAdapter::onMasterProcessor() co Line 96  bool MeshAdapter::onMasterProcessor() co
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 m_finleyMesh->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 {
115     return m_finleyMesh.get();     return m_finleyMesh.get();
116  }  }
# Line 112  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 228  void MeshAdapter::dump(const std::string Line 240  void MeshAdapter::dump(const std::string
240        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_FaceElements_numNodes to NetCDF file failed: " + *newFileName);
241     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )     if (!dataFile.add_att("num_ContactElements_numNodes",num_ContactElements_numNodes) )
242        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_ContactElements_numNodes to NetCDF file failed: " + *newFileName);
243     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Elements_TypeId", mesh->Elements->referenceElementSet->referenceElement->Type->TypeId) )
244        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending Elements_TypeId to NetCDF file failed: " + *newFileName);
245     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("FaceElements_TypeId", mesh->FaceElements->referenceElementSet->referenceElement->Type->TypeId) )
246        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending FaceElements_TypeId to NetCDF file failed: " + *newFileName);
247     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("ContactElements_TypeId", mesh->ContactElements->referenceElementSet->referenceElement->Type->TypeId) )
248        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending ContactElements_TypeId to NetCDF file failed: " + *newFileName);
249     if (!dataFile.add_att("Points_TypeId", mesh->Points->ReferenceElement->Type->TypeId) )     if (!dataFile.add_att("Points_TypeId", mesh->Points->referenceElementSet->referenceElement->Type->TypeId) )
250        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending Points_TypeId to NetCDF file failed: " + *newFileName);
251     if (!dataFile.add_att("num_Tags", num_Tags) )     if (!dataFile.add_att("num_Tags", num_Tags) )
252        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);        throw DataException("Error - MeshAdapter::dump: appending num_Tags to NetCDF file failed: " + *newFileName);
# Line 646  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 680  pair<int,int> MeshAdapter::getDataShape( Line 693  pair<int,int> MeshAdapter::getDataShape(
693     case(Elements):     case(Elements):
694     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
695        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
696        numDataPointsPerSample=mesh->Elements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
697     }     }
698     break;     break;
699     case(ReducedElements):     case(ReducedElements):
700     if (mesh->Elements!=NULL) {     if (mesh->Elements!=NULL) {
701        numSamples=mesh->Elements->numElements;        numSamples=mesh->Elements->numElements;
702        numDataPointsPerSample=mesh->Elements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->Elements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
703     }     }
704     break;     break;
705     case(FaceElements):     case(FaceElements):
706     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
707        numDataPointsPerSample=mesh->FaceElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
708        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
709     }     }
710     break;     break;
711     case(ReducedFaceElements):     case(ReducedFaceElements):
712     if (mesh->FaceElements!=NULL) {     if (mesh->FaceElements!=NULL) {
713        numDataPointsPerSample=mesh->FaceElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->FaceElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
714        numSamples=mesh->FaceElements->numElements;        numSamples=mesh->FaceElements->numElements;
715     }     }
716     break;     break;
# Line 709  pair<int,int> MeshAdapter::getDataShape( Line 722  pair<int,int> MeshAdapter::getDataShape(
722     break;     break;
723     case(ContactElementsZero):     case(ContactElementsZero):
724     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
725        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
726        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
727     }     }
728     break;     break;
729     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
730     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
731        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
732        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
733     }     }
734     break;     break;
735     case(ContactElementsOne):     case(ContactElementsOne):
736     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
737        numDataPointsPerSample=mesh->ContactElements->ReferenceElement->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElement->Parametrization->numQuadNodes;
738        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
739     }     }
740     break;     break;
741     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
742     if (mesh->ContactElements!=NULL) {     if (mesh->ContactElements!=NULL) {
743        numDataPointsPerSample=mesh->ContactElements->ReferenceElementReducedOrder->numQuadNodes;        numDataPointsPerSample=mesh->ContactElements->referenceElementSet->referenceElementReducedQuadrature->Parametrization->numQuadNodes;
744        numSamples=mesh->ContactElements->numElements;        numSamples=mesh->ContactElements->numElements;
745     }     }
746     break;     break;
# Line 883  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 1235  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 1444  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();
# Line 1452  void MeshAdapter::setNewX(const escript: Line 1468  void MeshAdapter::setNewX(const escript:
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 1502  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);
1555    }
1556    
1557     return;  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558    {
1559    #ifdef PASO_MPI
1560        index_t myFirstNode=0, myLastNode=0, k=0;
1561        index_t* globalNodeIndex=0;
1562        Finley_Mesh* mesh_p=m_finleyMesh.get();
1563        if (fs_code == FINLEY_REDUCED_NODES)
1564        {
1565        myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1566        myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1567        globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1568        }
1569        else
1570        {
1571        myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1572        myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1573        globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1574        }
1575        k=globalNodeIndex[id];
1576        return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1577    #endif
1578        return true;
1579  }  }
1580                                                                                                                                                                        
1581                                                                                                                                                                        
1582  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  
1583    //
1584    // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585    //
1586  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1587                                                   const int row_blocksize,                                                   const int row_blocksize,
1588                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
# Line 1592  SystemMatrixAdapter MeshAdapter::newSyst Line 1626  SystemMatrixAdapter MeshAdapter::newSyst
1626  #endif  #endif
1627     }     }
1628     else {     else {
1629        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1630     }     }
1631     checkPasoError();     checkPasoError();
1632     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1634  }  }
1635    
1636    //
1637  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1638    //
1639  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1640                                                           const double theta,                                                           const double theta,
1641                                                           const int blocksize,                                                           const int blocksize,
# Line 1665  bool MeshAdapter::isCellOriented(int fun Line 1702  bool MeshAdapter::isCellOriented(int fun
1702     return false;     return false;
1703  }  }
1704    
1705    bool
1706    MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1707    {
1708       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1709        class 1: DOF <-> Nodes
1710        class 2: ReducedDOF <-> ReducedNodes
1711        class 3: Points
1712        class 4: Elements
1713        class 5: ReducedElements
1714        class 6: FaceElements
1715        class 7: ReducedFaceElements
1716        class 8: ContactElementZero <-> ContactElementOne
1717        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1718    
1719       There is also a set of lines. Interpolation is possible down a line but not between lines.
1720       class 1 and 2 belong to all lines so aren't considered.
1721        line 0: class 3
1722        line 1: class 4,5
1723        line 2: class 6,7
1724        line 3: class 8,9
1725    
1726       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1727       eg hasnodes is true if we have at least one instance of Nodes.
1728       */
1729        if (fs.size()==0)
1730        {
1731        return false;
1732        }
1733        std::vector<int> hasclass(10);
1734        std::vector<int> hasline(4);    
1735        bool hasnodes=false;
1736        bool hasrednodes=false;
1737        bool hascez=false;
1738        bool hasrcez=false;
1739        for (int i=0;i<fs.size();++i)
1740        {
1741        switch(fs[i])
1742        {
1743        case(Nodes):   hasnodes=true;   // no break is deliberate
1744        case(DegreesOfFreedom):
1745            hasclass[1]=1;
1746            break;
1747        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1748        case(ReducedDegreesOfFreedom):
1749            hasclass[2]=1;
1750            break;
1751        case(Points):
1752            hasline[0]=1;
1753            hasclass[3]=1;
1754            break;
1755        case(Elements):
1756            hasclass[4]=1;
1757            hasline[1]=1;
1758            break;
1759        case(ReducedElements):
1760            hasclass[5]=1;
1761            hasline[1]=1;
1762            break;
1763        case(FaceElements):
1764            hasclass[6]=1;
1765            hasline[2]=1;
1766            break;
1767        case(ReducedFaceElements):
1768            hasclass[7]=1;
1769            hasline[2]=1;
1770            break;
1771        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1772        case(ContactElementsOne):
1773            hasclass[8]=1;
1774            hasline[3]=1;
1775            break;
1776        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1777        case(ReducedContactElementsOne):
1778            hasclass[9]=1;
1779            hasline[3]=1;
1780            break;
1781        default:
1782            return false;
1783        }
1784        }
1785        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1786        // fail if we have more than one leaf group
1787    
1788        if (totlines>1)
1789        {
1790        return false;   // there are at least two branches we can't interpolate between
1791        }
1792        else if (totlines==1)
1793        {
1794        if (hasline[0]==1)      // we have points
1795        {
1796            resultcode=Points;
1797        }
1798        else if (hasline[1]==1)
1799        {
1800            if (hasclass[5]==1)
1801            {
1802            resultcode=ReducedElements;
1803            }
1804            else
1805            {
1806            resultcode=Elements;
1807            }
1808        }
1809        else if (hasline[2]==1)
1810        {
1811            if (hasclass[7]==1)
1812            {
1813            resultcode=ReducedFaceElements;
1814            }
1815            else
1816            {
1817            resultcode=FaceElements;
1818            }
1819        }
1820        else    // so we must be in line3
1821        {
1822            if (hasclass[9]==1)
1823            {
1824            // need something from class 9
1825            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1826            }
1827            else
1828            {
1829            // something from class 8
1830            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1831            }
1832        }
1833        }
1834        else    // totlines==0
1835        {
1836        if (hasclass[2]==1)
1837        {
1838            // something from class 2
1839            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1840        }
1841        else
1842        {   // something from class 1
1843            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1844        }
1845        }
1846        return true;
1847    }
1848    
1849  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1850  {  {
1851     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1852     case(Nodes):     case(Nodes):
1853     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1854     case(Nodes):      case(Nodes):
1855     case(ReducedNodes):      case(ReducedNodes):
1856     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1857     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1858     case(Elements):      case(Elements):
1859     case(ReducedElements):      case(ReducedElements):
1860     case(FaceElements):      case(FaceElements):
1861     case(ReducedFaceElements):      case(ReducedFaceElements):
1862     case(Points):      case(Points):
1863     case(ContactElementsZero):      case(ContactElementsZero):
1864     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1865     case(ContactElementsOne):      case(ContactElementsOne):
1866     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1867     return true;      return true;
1868     default:      default:
1869        stringstream temp;            stringstream temp;
1870        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;
1871        throw FinleyAdapterException(temp.str());            throw FinleyAdapterException(temp.str());
    }  
    break;  
    case(ReducedNodes):  
    switch(functionSpaceType_target) {  
    case(ReducedNodes):  
    case(ReducedDegreesOfFreedom):  
    case(Elements):  
    case(ReducedElements):  
    case(FaceElements):  
    case(ReducedFaceElements):  
    case(Points):  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    return true;  
    case(Nodes):  
    case(DegreesOfFreedom):  
    return false;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;  
       throw FinleyAdapterException(temp.str());  
1872     }     }
1873     break;     break;
    case(Elements):  
    if (functionSpaceType_target==Elements) {  
       return true;  
    } else if (functionSpaceType_target==ReducedElements) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ReducedElements):  
    if (functionSpaceType_target==ReducedElements) {  
       return true;  
    } else {  
       return false;  
    }  
    case(FaceElements):  
    if (functionSpaceType_target==FaceElements) {  
       return true;  
    } else if (functionSpaceType_target==ReducedFaceElements) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ReducedFaceElements):  
    if (functionSpaceType_target==ReducedFaceElements) {  
       return true;  
    } else {  
       return false;  
    }  
    case(Points):  
    if (functionSpaceType_target==Points) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {  
       return true;  
    } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {  
       return true;  
    } else {  
       return false;  
    }  
    case(DegreesOfFreedom):  
    switch(functionSpaceType_target) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    case(Nodes):  
1874     case(ReducedNodes):     case(ReducedNodes):
1875     case(Elements):      switch(functionSpaceType_target) {
1876     case(ReducedElements):      case(ReducedNodes):
1877     case(Points):      case(ReducedDegreesOfFreedom):
1878     case(FaceElements):      case(Elements):
1879     case(ReducedFaceElements):      case(ReducedElements):
1880     case(ContactElementsZero):      case(FaceElements):
1881     case(ReducedContactElementsZero):      case(ReducedFaceElements):
1882     case(ContactElementsOne):      case(Points):
1883     case(ReducedContactElementsOne):      case(ContactElementsZero):
1884     return true;      case(ReducedContactElementsZero):
1885     default:      case(ContactElementsOne):
1886        stringstream temp;      case(ReducedContactElementsOne):
1887        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      return true;
1888        throw FinleyAdapterException(temp.str());      case(Nodes):
1889     }      case(DegreesOfFreedom):
1890     break;      return false;
1891        default:
1892            stringstream temp;
1893            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1894            throw FinleyAdapterException(temp.str());
1895       }
1896       break;
1897       case(Elements):
1898        if (functionSpaceType_target==Elements) {
1899          return true;
1900        } else if (functionSpaceType_target==ReducedElements) {
1901          return true;
1902            } else {
1903              return false;
1904            }
1905       case(ReducedElements):
1906        if (functionSpaceType_target==ReducedElements) {
1907          return true;
1908        } else {
1909              return false;
1910        }
1911       case(FaceElements):
1912        if (functionSpaceType_target==FaceElements) {
1913                return true;
1914        } else if (functionSpaceType_target==ReducedFaceElements) {
1915                return true;
1916        } else {
1917                return false;
1918        }
1919       case(ReducedFaceElements):
1920        if (functionSpaceType_target==ReducedFaceElements) {
1921                return true;
1922        } else {
1923            return false;
1924        }
1925       case(Points):
1926        if (functionSpaceType_target==Points) {
1927                return true;
1928        } else {
1929                return false;
1930        }
1931       case(ContactElementsZero):
1932       case(ContactElementsOne):
1933        if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1934                return true;
1935        } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1936                return true;
1937        } else {
1938                return false;
1939        }
1940       case(ReducedContactElementsZero):
1941       case(ReducedContactElementsOne):
1942        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1943                return true;
1944        } else {
1945                return false;
1946        }
1947       case(DegreesOfFreedom):
1948        switch(functionSpaceType_target) {
1949        case(ReducedDegreesOfFreedom):
1950        case(DegreesOfFreedom):
1951        case(Nodes):
1952        case(ReducedNodes):
1953        case(Elements):
1954        case(ReducedElements):
1955        case(Points):
1956        case(FaceElements):
1957        case(ReducedFaceElements):
1958        case(ContactElementsZero):
1959        case(ReducedContactElementsZero):
1960        case(ContactElementsOne):
1961        case(ReducedContactElementsOne):
1962        return true;
1963        default:
1964            stringstream temp;
1965            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1966            throw FinleyAdapterException(temp.str());
1967        }
1968        break;
1969     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1970     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1971     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1972     case(ReducedNodes):      case(ReducedNodes):
1973     case(Elements):      case(Elements):
1974     case(ReducedElements):      case(ReducedElements):
1975     case(FaceElements):      case(FaceElements):
1976     case(ReducedFaceElements):      case(ReducedFaceElements):
1977     case(Points):      case(Points):
1978     case(ContactElementsZero):      case(ContactElementsZero):
1979     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1980     case(ContactElementsOne):      case(ContactElementsOne):
1981     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1982     return true;      return true;
1983     case(Nodes):      case(Nodes):
1984     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1985     return false;      return false;
1986     default:      default:
1987        stringstream temp;          stringstream temp;
1988        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;
1989        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1990     }      }
1991     break;      break;
1992     default:     default:
1993        stringstream temp;        stringstream temp;
1994        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 1840  bool MeshAdapter::operator!=(const Abstr Line 2021  bool MeshAdapter::operator!=(const Abstr
2021    
2022  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2023  {  {
2024     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2025       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2026     checkPasoError();     checkPasoError();
2027     return out;     return out;
2028  }  }
2029  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2030  {  {
2031     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
2032       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2033     checkPasoError();     checkPasoError();
2034     return out;     return out;
2035  }  }
# Line 1866  escript::Data MeshAdapter::getSize() con Line 2049  escript::Data MeshAdapter::getSize() con
2049     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
2050  }  }
2051    
2052  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2053  {  {
2054     int *out = NULL;     int *out = NULL;
2055     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2103  int MeshAdapter::getNumberOfTagsInUse(in Line 2286  int MeshAdapter::getNumberOfTagsInUse(in
2286    }    }
2287    return numTags;    return numTags;
2288  }  }
2289  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2290    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2291  {  {
2292    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2293    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2169  bool MeshAdapter::canTag(int functionSpa Line 2353  bool MeshAdapter::canTag(int functionSpa
2353    }    }
2354  }  }
2355    
2356    AbstractDomain::StatusType MeshAdapter::getStatus() const
2357    {
2358      Finley_Mesh* mesh=m_finleyMesh.get();
2359      return Finley_Mesh_getStatus(mesh);
2360    }
2361    
2362    
2363    
2364  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1990  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26