/[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 1802 by jfenwick, Tue Sep 23 01:03:29 2008 UTC revision 2635 by jfenwick, Thu Aug 27 04:54:41 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  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
# Line 102  void MeshAdapter::write(const std::strin Line 109  void MeshAdapter::write(const std::strin
109     TMPMEMFREE(fName);     TMPMEMFREE(fName);
110  }  }
111    
112  void MeshAdapter::Print_Mesh_Info(const bool full=false) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
113  {  {
114     Finley_PrintMesh_Info(m_finleyMesh.get(), full);     Finley_PrintMesh_Info(m_finleyMesh.get(), full);
115  }  }
# Line 110  void MeshAdapter::Print_Mesh_Info(const Line 117  void MeshAdapter::Print_Mesh_Info(const
117  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const std::string& fileName) const
118  {  {
119  #ifdef USE_NETCDF  #ifdef USE_NETCDF
120     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
121     NcVar *ids;     NcVar *ids;
122     int *int_ptr;     int *int_ptr;
123     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
# Line 636  int MeshAdapter::getDiracDeltaFunctionCo Line 643  int MeshAdapter::getDiracDeltaFunctionCo
643  //  //
644  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
645  {  {
646     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
647       int numDim=Finley_Mesh_getDim(mesh);
648     checkFinleyError();     checkFinleyError();
649     return numDim;     return numDim;
650  }  }
# Line 861  void MeshAdapter::addPDEToTransportProbl Line 869  void MeshAdapter::addPDEToTransportProbl
869  //  //
870  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
871  {  {
872     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
873     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
874     if (inDomain!=*this)       if (inDomain!=*this)  
875        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
876     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 873  void MeshAdapter::interpolateOnDomain(es Line 881  void MeshAdapter::interpolateOnDomain(es
881     escriptDataC _in=in.getDataC();     escriptDataC _in=in.getDataC();
882     switch(in.getFunctionSpace().getTypeCode()) {     switch(in.getFunctionSpace().getTypeCode()) {
883     case(Nodes):     case(Nodes):
884     switch(target.getFunctionSpace().getTypeCode()) {        switch(target.getFunctionSpace().getTypeCode()) {
885     case(Nodes):        case(Nodes):
886     case(ReducedNodes):        case(ReducedNodes):
887     case(DegreesOfFreedom):        case(DegreesOfFreedom):
888     case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
889     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
    break;  
    case(Elements):  
    case(ReducedElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);  
    break;  
    case(FaceElements):  
    case(ReducedFaceElements):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);  
    break;  
    case(Points):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);  
    break;  
    case(ContactElementsZero):  
    case(ReducedContactElementsZero):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
890        break;        break;
891     }        case(Elements):
892     break;        case(ReducedElements):
893     case(ReducedNodes):        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
894     switch(target.getFunctionSpace().getTypeCode()) {        break;
895     case(Nodes):        case(FaceElements):
896     case(ReducedNodes):        case(ReducedFaceElements):
897     case(DegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
898     case(ReducedDegreesOfFreedom):        break;
899     Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        case(Points):
900     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
901     case(Elements):        break;
902     case(ReducedElements):        case(ContactElementsZero):
903     Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        case(ReducedContactElementsZero):
904     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
905     case(FaceElements):        break;
906     case(ReducedFaceElements):        case(ContactElementsOne):
907     Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        case(ReducedContactElementsOne):
908     break;        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
909     case(Points):        break;
910     Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        default:
911     break;           stringstream temp;
912     case(ContactElementsZero):           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
913     case(ReducedContactElementsZero):           throw FinleyAdapterException(temp.str());
914     Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);           break;
915     break;        }
    case(ContactElementsOne):  
    case(ReducedContactElementsOne):  
    Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
916        break;        break;
    }  
    break;  
    case(Elements):  
    if (target.getFunctionSpace().getTypeCode()==Elements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements possible.");  
    }  
    break;  
    case(ReducedElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedElements) {  
       Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");  
    }  
    break;  
    case(FaceElements):  
    if (target.getFunctionSpace().getTypeCode()==FaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");  
    }  
    break;  
    case(ReducedFaceElements):  
    if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {  
       Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");  
    }  
    break;  
    case(Points):  
    if (target.getFunctionSpace().getTypeCode()==Points) {  
       Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on points possible.");  
    }  
    break;  
    case(ContactElementsZero):  
    case(ContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");  
    }  
    break;  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {  
       Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);  
    } else {  
       throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");  
    }  
    break;  
    case(DegreesOfFreedom):        
    switch(target.getFunctionSpace().getTypeCode()) {  
    case(ReducedDegreesOfFreedom):  
    case(DegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
   
    case(Nodes):  
917     case(ReducedNodes):     case(ReducedNodes):
918     if (getMPISize()>1) {        switch(target.getFunctionSpace().getTypeCode()) {
919        escript::Data temp=escript::Data(in);        case(Nodes):
920        temp.expand();        case(ReducedNodes):
921        escriptDataC _in2 = temp.getDataC();        case(DegreesOfFreedom):
922        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);        case(ReducedDegreesOfFreedom):
    } else {  
923        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
924     }        break;
925     break;        case(Elements):
926     case(Elements):        case(ReducedElements):
    case(ReducedElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);  
    } else {  
927        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
928     }        break;
929     break;        case(FaceElements):
930     case(FaceElements):        case(ReducedFaceElements):
    case(ReducedFaceElements):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);  
   
    } else {  
931        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
932     }        break;
933     break;        case(Points):
    case(Points):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
    } else {  
934        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
935     }        break;
936     break;        case(ContactElementsZero):
937     case(ContactElementsZero):        case(ReducedContactElementsZero):
    case(ContactElementsOne):  
    case(ReducedContactElementsZero):  
    case(ReducedContactElementsOne):  
    if (getMPISize()>1) {  
       escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );  
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);  
    } else {  
938        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
    }  
    break;  
    default:  
       stringstream temp;  
       temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();  
       throw FinleyAdapterException(temp.str());  
939        break;        break;
940     }        case(ContactElementsOne):
941     break;        case(ReducedContactElementsOne):
942     case(ReducedDegreesOfFreedom):        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
943     switch(target.getFunctionSpace().getTypeCode()) {        break;
944     case(Nodes):        default:
945     throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");           stringstream temp;
946     break;           temp << "Error - Interpolation on Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
947     case(ReducedNodes):           throw FinleyAdapterException(temp.str());
948     if (getMPISize()>1) {           break;
949        escript::Data temp=escript::Data(in);        }
950        temp.expand();        break;
       escriptDataC _in2 = temp.getDataC();  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);  
    } else {  
       Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    }  
    break;  
    case(DegreesOfFreedom):  
    throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");  
    break;  
    case(ReducedDegreesOfFreedom):  
    Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);  
    break;  
951     case(Elements):     case(Elements):
952          if (target.getFunctionSpace().getTypeCode()==Elements) {
953             Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
954          } else if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
955             Finley_Assemble_AverageElementData(mesh->Elements,&_target,&_in);
956          } else {
957             throw FinleyAdapterException("Error - No interpolation with data on elements possible.");
958          }
959          break;
960     case(ReducedElements):     case(ReducedElements):
961     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedElements) {
962        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Elements,&_target,&_in);
963        escriptDataC _in2 = temp.getDataC();        } else {
964        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on elements with reduced integration order possible.");
965     } else {        }
966        Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);        break;
    }  
    break;  
967     case(FaceElements):     case(FaceElements):
968          if (target.getFunctionSpace().getTypeCode()==FaceElements) {
969             Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
970          } else if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
971             Finley_Assemble_AverageElementData(mesh->FaceElements,&_target,&_in);
972          } else {
973             throw FinleyAdapterException("Error - No interpolation with data on face elements possible.");
974          }
975          break;
976     case(ReducedFaceElements):     case(ReducedFaceElements):
977     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedFaceElements) {
978        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->FaceElements,&_target,&_in);
979        escriptDataC _in2 = temp.getDataC();        } else {
980        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on face elements with reduced integration order possible.");
981     } else {        }
982        Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);        break;
    }  
    break;  
983     case(Points):     case(Points):
984     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==Points) {
985        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->Points,&_target,&_in);
986        escriptDataC _in2 = temp.getDataC();        } else {
987        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on points possible.");
988     } else {        }
989        Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);        break;
    }  
    break;  
990     case(ContactElementsZero):     case(ContactElementsZero):
991     case(ContactElementsOne):     case(ContactElementsOne):
992          if (target.getFunctionSpace().getTypeCode()==ContactElementsZero || target.getFunctionSpace().getTypeCode()==ContactElementsOne) {
993             Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
994          } else if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
995             Finley_Assemble_AverageElementData(mesh->ContactElements,&_target,&_in);
996          } else {
997             throw FinleyAdapterException("Error - No interpolation with data on contact elements possible.");
998          }
999          break;
1000     case(ReducedContactElementsZero):     case(ReducedContactElementsZero):
1001     case(ReducedContactElementsOne):     case(ReducedContactElementsOne):
1002     if (getMPISize()>1) {        if (target.getFunctionSpace().getTypeCode()==ReducedContactElementsZero || target.getFunctionSpace().getTypeCode()==ReducedContactElementsOne) {
1003        escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           Finley_Assemble_CopyElementData(mesh->ContactElements,&_target,&_in);
1004        escriptDataC _in2 = temp.getDataC();        } else {
1005        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           throw FinleyAdapterException("Error - No interpolation with data on contact elements with reduced integration order possible.");
1006     } else {        }
1007        Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);        break;
1008     }     case(DegreesOfFreedom):      
1009     break;        switch(target.getFunctionSpace().getTypeCode()) {
1010     default:        case(ReducedDegreesOfFreedom):
1011        stringstream temp;        case(DegreesOfFreedom):
1012        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();        Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1013        throw FinleyAdapterException(temp.str());        break;
1014      
1015          case(Nodes):
1016          case(ReducedNodes):
1017          if (getMPISize()>1) {
1018             escript::Data temp=escript::Data(in);
1019             temp.expand();
1020             escriptDataC _in2 = temp.getDataC();
1021             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1022          } else {
1023             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1024          }
1025          break;
1026          case(Elements):
1027          case(ReducedElements):
1028          if (getMPISize()>1) {
1029             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1030             escriptDataC _in2 = temp.getDataC();
1031             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1032          } else {
1033             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1034          }
1035          break;
1036          case(FaceElements):
1037          case(ReducedFaceElements):
1038          if (getMPISize()>1) {
1039             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1040             escriptDataC _in2 = temp.getDataC();
1041             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1042      
1043          } else {
1044             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1045          }
1046          break;
1047          case(Points):
1048          if (getMPISize()>1) {
1049             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1050             escriptDataC _in2 = temp.getDataC();
1051          } else {
1052             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1053          }
1054          break;
1055          case(ContactElementsZero):
1056          case(ContactElementsOne):
1057          case(ReducedContactElementsZero):
1058          case(ReducedContactElementsOne):
1059          if (getMPISize()>1) {
1060             escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );
1061             escriptDataC _in2 = temp.getDataC();
1062             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1063          } else {
1064             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1065          }
1066          break;
1067          default:
1068             stringstream temp;
1069             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1070             throw FinleyAdapterException(temp.str());
1071             break;
1072          }
1073          break;
1074       case(ReducedDegreesOfFreedom):
1075          switch(target.getFunctionSpace().getTypeCode()) {
1076          case(Nodes):
1077          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to mesh nodes.");
1078          break;
1079          case(ReducedNodes):
1080          if (getMPISize()>1) {
1081             escript::Data temp=escript::Data(in);
1082             temp.expand();
1083             escriptDataC _in2 = temp.getDataC();
1084             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in2);
1085          } else {
1086             Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1087          }
1088          break;
1089          case(DegreesOfFreedom):
1090          throw FinleyAdapterException("Error - Finley does not support interpolation from reduced degrees of freedom to degrees of freedom");
1091          break;
1092          case(ReducedDegreesOfFreedom):
1093          Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
1094          break;
1095          case(Elements):
1096          case(ReducedElements):
1097          if (getMPISize()>1) {
1098             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1099             escriptDataC _in2 = temp.getDataC();
1100             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1101          } else {
1102             Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
1103          }
1104          break;
1105          case(FaceElements):
1106          case(ReducedFaceElements):
1107          if (getMPISize()>1) {
1108             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1109             escriptDataC _in2 = temp.getDataC();
1110             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1111          } else {
1112             Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in,&_target);
1113          }
1114          break;
1115          case(Points):
1116          if (getMPISize()>1) {
1117             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1118             escriptDataC _in2 = temp.getDataC();
1119             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1120          } else {
1121             Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1122          }
1123          break;
1124          case(ContactElementsZero):
1125          case(ContactElementsOne):
1126          case(ReducedContactElementsZero):
1127          case(ReducedContactElementsOne):
1128          if (getMPISize()>1) {
1129             escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );
1130             escriptDataC _in2 = temp.getDataC();
1131             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1132          } else {
1133             Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
1134          }
1135          break;
1136          default:
1137             stringstream temp;
1138             temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
1139             throw FinleyAdapterException(temp.str());
1140             break;
1141          }
1142        break;        break;
    }  
    break;  
1143     default:     default:
1144        stringstream temp;        stringstream temp;
1145        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type %d" << in.getFunctionSpace().getTypeCode();
# Line 1146  void MeshAdapter::interpolateOnDomain(es Line 1154  void MeshAdapter::interpolateOnDomain(es
1154  //  //
1155  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1156  {  {
1157     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1158     if (argDomain!=*this)     if (argDomain!=*this)
1159        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1160     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1169  void MeshAdapter::setToX(escript::Data& Line 1177  void MeshAdapter::setToX(escript::Data&
1177  //  //
1178  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1179  {  {
1180     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1181       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1182     if (normalDomain!=*this)     if (normalDomain!=*this)
1183        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1184     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1224  void MeshAdapter::setToNormal(escript::D Line 1233  void MeshAdapter::setToNormal(escript::D
1233  //  //
1234  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1235  {  {
1236     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1237     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1238       if (targetDomain!=this)
1239        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1240    
1241     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1236  void MeshAdapter::interpolateACross(escr Line 1246  void MeshAdapter::interpolateACross(escr
1246  //  //
1247  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1248  {  {
1249     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1250     if (argDomain!=*this)     if (argDomain!=*this)
1251        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1252    
# Line 1247  void MeshAdapter::setToIntegrals(std::ve Line 1257  void MeshAdapter::setToIntegrals(std::ve
1257     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1258     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1259     case(Nodes):     case(Nodes):
1260     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1261     _temp=temp.getDataC();     _temp=temp.getDataC();
1262     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1263     break;     break;
1264     case(ReducedNodes):     case(ReducedNodes):
1265     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1266     _temp=temp.getDataC();     _temp=temp.getDataC();
1267     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1268     break;     break;
# Line 1284  void MeshAdapter::setToIntegrals(std::ve Line 1294  void MeshAdapter::setToIntegrals(std::ve
1294     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1295     break;     break;
1296     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1297     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1298     _temp=temp.getDataC();     _temp=temp.getDataC();
1299     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1300     break;     break;
1301     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1302     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1303     _temp=temp.getDataC();     _temp=temp.getDataC();
1304     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1305     break;     break;
# Line 1308  void MeshAdapter::setToIntegrals(std::ve Line 1318  void MeshAdapter::setToIntegrals(std::ve
1318  //  //
1319  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1320  {  {
1321     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1322     if (argDomain!=*this)     if (argDomain!=*this)
1323        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1324     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1325     if (gradDomain!=*this)     if (gradDomain!=*this)
1326        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1327    
# Line 1433  void MeshAdapter::setToSize(escript::Dat Line 1443  void MeshAdapter::setToSize(escript::Dat
1443     checkFinleyError();     checkFinleyError();
1444  }  }
1445    
1446  // sets the location of nodes:  //
1447    // sets the location of nodes
1448    //
1449  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1450  {  {
1451     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1452     escriptDataC tmp;     escriptDataC tmp;
1453     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1454     if (newDomain!=*this)     if (newDomain!=*this)
1455        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1456     tmp = new_x.getDataC();     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {
1457     Finley_Mesh_setCoordinates(mesh,&tmp);         tmp = new_x.getDataC();
1458           Finley_Mesh_setCoordinates(mesh,&tmp);
1459       } else {
1460           escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );
1461           tmp = new_x_inter.getDataC();
1462           Finley_Mesh_setCoordinates(mesh,&tmp);
1463       }
1464     checkFinleyError();     checkFinleyError();
1465  }  }
1466    
1467  // saves a data array in openDX format:  //
1468  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  // Helper for the save* methods. Extracts optional data variable names and the
1469    // corresponding pointers from python dictionary. Caller must free arrays.
1470    //
1471    void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const
1472  {  {
1473     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     numData = boost::python::extract<int>(arg.attr("__len__")());
    const int num_data=boost::python::extract<int>(arg.attr("__len__")());  
1474     /* win32 refactor */     /* win32 refactor */
1475     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;
1476     for(int i=0;i<num_data;i++)     data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;
1477     {     dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
1478    
1479     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1480     for (int i=0;i<num_data;++i) {     for (int i=0; i<numData; ++i) {
1481        std::string n=boost::python::extract<std::string>(keys[i]);        std::string n=boost::python::extract<std::string>(keys[i]);
1482        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1483        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1484           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error: Data must be defined on same Domain");
1485        data[i]=d.getDataC();        data[i] = d.getDataC();
1486        ptr_data[i]= &(data[i]);        dataPtr[i] = &(data[i]);
1487        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1488           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1489     }     }
1490     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  }
1491    
1492    //
1493    // saves mesh and optionally data arrays in openDX format
1494    //
1495    void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1496    {
1497       int num_data;
1498       char **names;
1499       escriptDataC *data;
1500       escriptDataC **ptr_data;
1501    
1502       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1503       Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);
1504     checkFinleyError();     checkFinleyError();
1505    
1506     /* win32 refactor */     /* win32 refactor */
1507     TMPMEMFREE(data);     TMPMEMFREE(data);
1508     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1509     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1510     {     {
1511        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1512     }     }
# Line 1491  void MeshAdapter::saveDX(const std::stri Line 1515  void MeshAdapter::saveDX(const std::stri
1515     return;     return;
1516  }  }
1517    
1518  // saves a data array in openVTK format:  //
1519  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  // saves mesh and optionally data arrays in VTK format
1520    //
1521    void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg,  const std::string& metadata, const std::string& metadata_schema) const
1522  {  {
1523     unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;     int num_data;
1524     const int num_data=boost::python::extract<int>(arg.attr("__len__")());     char **names;
1525     /* win32 refactor */     escriptDataC *data;
1526     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     escriptDataC **ptr_data;
    for(int i=0;i<num_data;i++)  
    {  
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
   
    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;  
    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0;i<num_data;++i) {  
       std::string n=boost::python::extract<std::string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)  
          throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");  
       data[i]=d.getDataC();  
       ptr_data[i]=&(data[i]);  
       if (n.length()>MAX_namelength-1) {  
          strncpy(names[i],n.c_str(),MAX_namelength-1);  
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
    }  
    Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);  
1527    
1528       extractArgsFromDict(arg, num_data, names, data, ptr_data);
1529       Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());
1530     checkFinleyError();     checkFinleyError();
1531    
1532     /* win32 refactor */     /* win32 refactor */
1533     TMPMEMFREE(data);     TMPMEMFREE(data);
1534     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1535     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1536     {     {
1537        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1538     }     }
1539     TMPMEMFREE(names);     TMPMEMFREE(names);
   
    return;  
1540  }  }
1541                                                                                                                                                                        
1542                                                                                                                                                                        //
1543  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros:  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1544    //
1545  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  SystemMatrixAdapter MeshAdapter::newSystemMatrix(
1546                                                   const int row_blocksize,                                                   const int row_blocksize,
1547                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
# Line 1548  SystemMatrixAdapter MeshAdapter::newSyst Line 1552  SystemMatrixAdapter MeshAdapter::newSyst
1552     int reduceRowOrder=0;     int reduceRowOrder=0;
1553     int reduceColOrder=0;     int reduceColOrder=0;
1554     // is the domain right?     // is the domain right?
1555     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1556     if (row_domain!=*this)     if (row_domain!=*this)
1557        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1558     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1559     if (col_domain!=*this)     if (col_domain!=*this)
1560        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1561     // is the function space type right     // is the function space type right
# Line 1581  SystemMatrixAdapter MeshAdapter::newSyst Line 1585  SystemMatrixAdapter MeshAdapter::newSyst
1585  #endif  #endif
1586     }     }
1587     else {     else {
1588        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize,FALSE);
1589     }     }
1590     checkPasoError();     checkPasoError();
1591     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1592     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1593  }  }
1594    
1595    //
1596  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1597    //
1598  TransportProblemAdapter MeshAdapter::newTransportProblem(  TransportProblemAdapter MeshAdapter::newTransportProblem(
1599                                                           const double theta,                                                           const double theta,
1600                                                           const int blocksize,                                                           const int blocksize,
# Line 1596  TransportProblemAdapter MeshAdapter::new Line 1603  TransportProblemAdapter MeshAdapter::new
1603  {  {
1604     int reduceOrder=0;     int reduceOrder=0;
1605     // is the domain right?     // is the domain right?
1606     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1607     if (domain!=*this)     if (domain!=*this)
1608        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1609     // is the function space type right     // is the function space type right
# Line 1654  bool MeshAdapter::isCellOriented(int fun Line 1661  bool MeshAdapter::isCellOriented(int fun
1661     return false;     return false;
1662  }  }
1663    
1664    bool
1665    MeshAdapter::commonFunctionSpace(const std::vector<int>& fs, int& resultcode) const
1666    {
1667       /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1668        class 1: DOF <-> Nodes
1669        class 2: ReducedDOF <-> ReducedNodes
1670        class 3: Points
1671        class 4: Elements
1672        class 5: ReducedElements
1673        class 6: FaceElements
1674        class 7: ReducedFaceElements
1675        class 8: ContactElementZero <-> ContactElementOne
1676        class 9: ReducedContactElementZero <-> ReducedContactElementOne
1677    
1678       There is also a set of lines. Interpolation is possible down a line but not between lines.
1679       class 1 and 2 belong to all lines so aren't considered.
1680        line 0: class 3
1681        line 1: class 4,5
1682        line 2: class 6,7
1683        line 3: class 8,9
1684    
1685       For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1686       eg hasnodes is true if we have at least one instance of Nodes.
1687       */
1688        if (fs.size()==0)
1689        {
1690        return false;
1691        }
1692        std::vector<int> hasclass(10);
1693        std::vector<int> hasline(4);    
1694        bool hasnodes=false;
1695        bool hasrednodes=false;
1696        bool hascez=false;
1697        bool hasrcez=false;
1698        for (int i=0;i<fs.size();++i)
1699        {
1700        switch(fs[i])
1701        {
1702        case(Nodes):   hasnodes=true;   // no break is deliberate
1703        case(DegreesOfFreedom):
1704            hasclass[1]=1;
1705            break;
1706        case(ReducedNodes):    hasrednodes=true;    // no break is deliberate
1707        case(ReducedDegreesOfFreedom):
1708            hasclass[2]=1;
1709            break;
1710        case(Points):
1711            hasline[0]=1;
1712            hasclass[3]=1;
1713            break;
1714        case(Elements):
1715            hasclass[4]=1;
1716            hasline[1]=1;
1717            break;
1718        case(ReducedElements):
1719            hasclass[5]=1;
1720            hasline[1]=1;
1721            break;
1722        case(FaceElements):
1723            hasclass[6]=1;
1724            hasline[2]=1;
1725            break;
1726        case(ReducedFaceElements):
1727            hasclass[7]=1;
1728            hasline[2]=1;
1729            break;
1730        case(ContactElementsZero):  hascez=true;    // no break is deliberate
1731        case(ContactElementsOne):
1732            hasclass[8]=1;
1733            hasline[3]=1;
1734            break;
1735        case(ReducedContactElementsZero):   hasrcez=true;   // no break is deliberate
1736        case(ReducedContactElementsOne):
1737            hasclass[9]=1;
1738            hasline[3]=1;
1739            break;
1740        default:
1741            return false;
1742        }
1743        }
1744        int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1745        // fail if we have more than one leaf group
1746    
1747        if (totlines>1)
1748        {
1749        return false;   // there are at least two branches we can't interpolate between
1750        }
1751        else if (totlines==1)
1752        {
1753        if (hasline[0]==1)      // we have points
1754        {
1755            resultcode=Points;
1756        }
1757        else if (hasline[1]==1)
1758        {
1759            if (hasclass[5]==1)
1760            {
1761            resultcode=ReducedElements;
1762            }
1763            else
1764            {
1765            resultcode=Elements;
1766            }
1767        }
1768        else if (hasline[2]==1)
1769        {
1770            if (hasclass[7]==ReducedFaceElements)
1771            {
1772            resultcode=ReducedFaceElements;
1773            }
1774            else
1775            {
1776            resultcode=FaceElements;
1777            }
1778        }
1779        else    // so we must be in line3
1780        {
1781            if (hasclass[9]==1)
1782            {
1783            // need something from class 9
1784            resultcode=(hasrcez?ReducedContactElementsZero:ReducedContactElementsOne);
1785            }
1786            else
1787            {
1788            // something from class 8
1789            resultcode=(hascez?ContactElementsZero:ContactElementsOne);
1790            }
1791        }
1792        }
1793        else    // totlines==0
1794        {
1795        if (hasclass[2]==1)
1796        {
1797            // something from class 2
1798            resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1799        }
1800        else
1801        {   // something from class 1
1802            resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1803        }
1804        }
1805        return true;
1806    }
1807    
1808  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1809  {  {
1810     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1811     case(Nodes):     case(Nodes):
1812     switch(functionSpaceType_target) {      switch(functionSpaceType_target) {
1813     case(Nodes):      case(Nodes):
1814     case(ReducedNodes):      case(ReducedNodes):
1815     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1816     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1817     case(Elements):      case(Elements):
1818     case(ReducedElements):      case(ReducedElements):
1819     case(FaceElements):      case(FaceElements):
1820     case(ReducedFaceElements):      case(ReducedFaceElements):
1821     case(Points):      case(Points):
1822     case(ContactElementsZero):      case(ContactElementsZero):
1823     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1824     case(ContactElementsOne):      case(ContactElementsOne):
1825     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1826     return true;      return true;
1827     default:      default:
1828        stringstream temp;            stringstream temp;
1829        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;
1830        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());  
1831     }     }
1832     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):  
1833     case(ReducedNodes):     case(ReducedNodes):
1834     case(Elements):      switch(functionSpaceType_target) {
1835     case(ReducedElements):      case(ReducedNodes):
1836     case(Points):      case(ReducedDegreesOfFreedom):
1837     case(FaceElements):      case(Elements):
1838     case(ReducedFaceElements):      case(ReducedElements):
1839     case(ContactElementsZero):      case(FaceElements):
1840     case(ReducedContactElementsZero):      case(ReducedFaceElements):
1841     case(ContactElementsOne):      case(Points):
1842     case(ReducedContactElementsOne):      case(ContactElementsZero):
1843     return true;      case(ReducedContactElementsZero):
1844     default:      case(ContactElementsOne):
1845        stringstream temp;      case(ReducedContactElementsOne):
1846        temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;      return true;
1847        throw FinleyAdapterException(temp.str());      case(Nodes):
1848     }      case(DegreesOfFreedom):
1849     break;      return false;
1850        default:
1851            stringstream temp;
1852            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1853            throw FinleyAdapterException(temp.str());
1854       }
1855       break;
1856       case(Elements):
1857        if (functionSpaceType_target==Elements) {
1858          return true;
1859        } else if (functionSpaceType_target==ReducedElements) {
1860          return true;
1861            } else {
1862              return false;
1863            }
1864       case(ReducedElements):
1865        if (functionSpaceType_target==ReducedElements) {
1866          return true;
1867        } else {
1868              return false;
1869        }
1870       case(FaceElements):
1871        if (functionSpaceType_target==FaceElements) {
1872                return true;
1873        } else if (functionSpaceType_target==ReducedFaceElements) {
1874                return true;
1875        } else {
1876                return false;
1877        }
1878       case(ReducedFaceElements):
1879        if (functionSpaceType_target==ReducedFaceElements) {
1880                return true;
1881        } else {
1882            return false;
1883        }
1884       case(Points):
1885        if (functionSpaceType_target==Points) {
1886                return true;
1887        } else {
1888                return false;
1889        }
1890       case(ContactElementsZero):
1891       case(ContactElementsOne):
1892        if (functionSpaceType_target==ContactElementsZero || functionSpaceType_target==ContactElementsOne) {
1893                return true;
1894        } else if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1895                return true;
1896        } else {
1897                return false;
1898        }
1899       case(ReducedContactElementsZero):
1900       case(ReducedContactElementsOne):
1901        if (functionSpaceType_target==ReducedContactElementsZero || functionSpaceType_target==ReducedContactElementsOne) {
1902                return true;
1903        } else {
1904                return false;
1905        }
1906       case(DegreesOfFreedom):
1907        switch(functionSpaceType_target) {
1908        case(ReducedDegreesOfFreedom):
1909        case(DegreesOfFreedom):
1910        case(Nodes):
1911        case(ReducedNodes):
1912        case(Elements):
1913        case(ReducedElements):
1914        case(Points):
1915        case(FaceElements):
1916        case(ReducedFaceElements):
1917        case(ContactElementsZero):
1918        case(ReducedContactElementsZero):
1919        case(ContactElementsOne):
1920        case(ReducedContactElementsOne):
1921        return true;
1922        default:
1923            stringstream temp;
1924            temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << functionSpaceType_target;
1925            throw FinleyAdapterException(temp.str());
1926        }
1927        break;
1928     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1929     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1930     case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1931     case(ReducedNodes):      case(ReducedNodes):
1932     case(Elements):      case(Elements):
1933     case(ReducedElements):      case(ReducedElements):
1934     case(FaceElements):      case(FaceElements):
1935     case(ReducedFaceElements):      case(ReducedFaceElements):
1936     case(Points):      case(Points):
1937     case(ContactElementsZero):      case(ContactElementsZero):
1938     case(ReducedContactElementsZero):      case(ReducedContactElementsZero):
1939     case(ContactElementsOne):      case(ContactElementsOne):
1940     case(ReducedContactElementsOne):      case(ReducedContactElementsOne):
1941     return true;      return true;
1942     case(Nodes):      case(Nodes):
1943     case(DegreesOfFreedom):      case(DegreesOfFreedom):
1944     return false;      return false;
1945     default:      default:
1946        stringstream temp;          stringstream temp;
1947        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;
1948        throw FinleyAdapterException(temp.str());          throw FinleyAdapterException(temp.str());
1949     }      }
1950     break;      break;
1951     default:     default:
1952        stringstream temp;        stringstream temp;
1953        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 1978  bool MeshAdapter::operator!=(const Abstr
1978     return !(operator==(other));     return !(operator==(other));
1979  }  }
1980    
1981  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
1982  {  {
1983     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     Finley_Mesh* mesh=m_finleyMesh.get();
1984       int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1985       checkPasoError();
1986       return out;
1987    }
1988    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1989    {
1990       Finley_Mesh* mesh=m_finleyMesh.get();
1991       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
1992     checkPasoError();     checkPasoError();
1993     return out;     return out;
1994  }  }
# Line 1846  escript::Data MeshAdapter::getNormal() c Line 2005  escript::Data MeshAdapter::getNormal() c
2005    
2006  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2007  {  {
2008     return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
2009  }  }
2010    
2011  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
2012  {  {
2013     int *out = NULL;     int *out = NULL;
2014     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 2086  int MeshAdapter::getNumberOfTagsInUse(in Line 2245  int MeshAdapter::getNumberOfTagsInUse(in
2245    }    }
2246    return numTags;    return numTags;
2247  }  }
2248  int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  
2249    const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2250  {  {
2251    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
2252    index_t* tags=NULL;    index_t* tags=NULL;
# Line 2152  bool MeshAdapter::canTag(int functionSpa Line 2312  bool MeshAdapter::canTag(int functionSpa
2312    }    }
2313  }  }
2314    
2315    AbstractDomain::StatusType MeshAdapter::getStatus() const
2316    {
2317      Finley_Mesh* mesh=m_finleyMesh.get();
2318      return Finley_Mesh_getStatus(mesh);
2319    }
2320    
2321    
2322    
2323  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1802  
changed lines
  Added in v.2635

  ViewVC Help
Powered by ViewVC 1.1.26