/[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 2898 by caltinay, Mon Feb 1 01:23:46 2010 UTC revision 4346 by jfenwick, Tue Apr 2 04:46:45 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
 extern "C" {  
24  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
25  }  
26    #include <boost/python/import.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    using namespace paso;
31    
32  namespace finley {  namespace finley {
33    
# Line 85  int MeshAdapter::getMPIRank() const Line 85  int MeshAdapter::getMPIRank() const
85  }  }
86  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
87  {  {
88  #ifdef PASO_MPI  #ifdef ESYS_MPI
89     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
90  #endif  #endif
91     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 96  bool MeshAdapter::onMasterProcessor() co
96  }  }
97    
98    
99  #ifdef PASO_MPI  #ifdef ESYS_MPI
100    MPI_Comm    MPI_Comm
101  #else  #else
102    unsigned int    unsigned int
103  #endif  #endif
104  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
105  {  {
106  #ifdef PASO_MPI  #ifdef ESYS_MPI
107      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
108  #else  #else
109      return 0;      return 0;
# Line 149  void MeshAdapter::dump(const string& fil Line 149  void MeshAdapter::dump(const string& fil
149     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
150     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
151     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
152  #ifdef PASO_MPI  #ifdef ESYS_MPI
153     MPI_Status status;     MPI_Status status;
154  #endif  #endif
155    
156  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
157  #ifdef PASO_MPI  #ifdef ESYS_MPI
158     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
159  #endif  #endif
160    
161     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
162                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
163    
164     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
# Line 529  void MeshAdapter::dump(const string& fil Line 529  void MeshAdapter::dump(const string& fil
529     }     }
530    
531  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
532  #ifdef PASO_MPI  #ifdef ESYS_MPI
533     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
534  #endif  #endif
535    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 649  int MeshAdapter::getReducedSolutionCode(
649     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
650  }  }
651    
652  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
653  {  {
654     return Points;     return Points;
655  }  }
# Line 770  pair<int,int> MeshAdapter::getDataShape( Line 770  pair<int,int> MeshAdapter::getDataShape(
770  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
771  //  //
772  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
773                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
774                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
775                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
776                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
777                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
778  {  {
779       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
780       if (smat==0)
781       {
782        throw FinleyAdapterException("finley only supports Paso system matrices.");
783       }
784     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
785     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
786     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 786  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
793     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
794     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
795       escriptDataC _d_dirac=d_dirac.getDataC();
796       escriptDataC _y_dirac=y_dirac.getDataC();
797    
798     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
799    
800     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
801       checkFinleyError();
802    
803       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
804     checkFinleyError();     checkFinleyError();
805    
806     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
807     checkFinleyError();     checkFinleyError();
808    
809     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );      Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
810     checkFinleyError();     checkFinleyError();
811  }  }
812    
813  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
814                                          escript::Data& mat,                                          escript::Data& mat,
815                                          const escript::Data& D,                                          const escript::Data& D,
816                                          const escript::Data& d) const                                          const escript::Data& d,
817                                            const escript::Data& d_dirac,
818                                            const bool useHRZ) const
819  {  {
820     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
821     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
822     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
823       escriptDataC _d_dirac=d_dirac.getDataC();
824    
825     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
826    
827     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
828     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
829      
830       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
831       checkFinleyError();
832    
833       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
834     checkFinleyError();     checkFinleyError();
835    
836  }  }
837    
838    
839  //  //
840  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
841  //  //
842  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
843  {  {
844     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
845    
# Line 829  void MeshAdapter::addPDEToRHS( escript:: Line 848  void MeshAdapter::addPDEToRHS( escript::
848     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
849     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
850     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
851       escriptDataC _y_dirac=y_dirac.getDataC();
852    
853     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
854     checkFinleyError();     checkFinleyError();
# Line 838  void MeshAdapter::addPDEToRHS( escript:: Line 858  void MeshAdapter::addPDEToRHS( escript::
858    
859     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
860     checkFinleyError();     checkFinleyError();
861    
862       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
863       checkFinleyError();
864    
865  }  }
866  //  //
867  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
868  //  //
869  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
870                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
871                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
872                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
873                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
874                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
875                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
876  {  {
877       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
878       if (tpa==0)
879       {
880        throw FinleyAdapterException("finley only supports Paso transport problems.");
881       }
882    
883    
884     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
885     source.expand();     source.expand();
886     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 863  void MeshAdapter::addPDEToTransportProbl Line 895  void MeshAdapter::addPDEToTransportProbl
895     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
896     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
897     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
898       escriptDataC _d_dirac=d_dirac.getDataC();
899       escriptDataC _y_dirac=y_dirac.getDataC();
900    
901    
902     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
903     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
904    
905     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
906     checkFinleyError();     checkFinleyError();
# Line 878  void MeshAdapter::addPDEToTransportProbl Line 913  void MeshAdapter::addPDEToTransportProbl
913    
914     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
915     checkFinleyError();     checkFinleyError();
916    
917       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
918       checkFinleyError();
919    
920  }  }
921    
922  //  //
# Line 1042  void MeshAdapter::interpolateOnDomain(es Line 1081  void MeshAdapter::interpolateOnDomain(es
1081        case(Elements):        case(Elements):
1082        case(ReducedElements):        case(ReducedElements):
1083        if (getMPISize()>1) {        if (getMPISize()>1) {
1084           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1085           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1086           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1087        } else {        } else {
# Line 1052  void MeshAdapter::interpolateOnDomain(es Line 1091  void MeshAdapter::interpolateOnDomain(es
1091        case(FaceElements):        case(FaceElements):
1092        case(ReducedFaceElements):        case(ReducedFaceElements):
1093        if (getMPISize()>1) {        if (getMPISize()>1) {
1094           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1095           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1096           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1097        
# Line 1062  void MeshAdapter::interpolateOnDomain(es Line 1101  void MeshAdapter::interpolateOnDomain(es
1101        break;        break;
1102        case(Points):        case(Points):
1103        if (getMPISize()>1) {        if (getMPISize()>1) {
1104           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1105           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1106        } else {        } else {
1107           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1108        }        }
# Line 1073  void MeshAdapter::interpolateOnDomain(es Line 1112  void MeshAdapter::interpolateOnDomain(es
1112        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1113        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1114        if (getMPISize()>1) {        if (getMPISize()>1) {
1115           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1116           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1117           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1118        } else {        } else {
# Line 1111  void MeshAdapter::interpolateOnDomain(es Line 1150  void MeshAdapter::interpolateOnDomain(es
1150        case(Elements):        case(Elements):
1151        case(ReducedElements):        case(ReducedElements):
1152        if (getMPISize()>1) {        if (getMPISize()>1) {
1153           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1154           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1155           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1156        } else {        } else {
# Line 1121  void MeshAdapter::interpolateOnDomain(es Line 1160  void MeshAdapter::interpolateOnDomain(es
1160        case(FaceElements):        case(FaceElements):
1161        case(ReducedFaceElements):        case(ReducedFaceElements):
1162        if (getMPISize()>1) {        if (getMPISize()>1) {
1163           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1164           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1165           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1166        } else {        } else {
# Line 1130  void MeshAdapter::interpolateOnDomain(es Line 1169  void MeshAdapter::interpolateOnDomain(es
1169        break;        break;
1170        case(Points):        case(Points):
1171        if (getMPISize()>1) {        if (getMPISize()>1) {
1172           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1173           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1174           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1175        } else {        } else {
# Line 1142  void MeshAdapter::interpolateOnDomain(es Line 1181  void MeshAdapter::interpolateOnDomain(es
1181        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1182        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1183        if (getMPISize()>1) {        if (getMPISize()>1) {
1184           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1185           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1186           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1187        } else {        } else {
# Line 1179  void MeshAdapter::setToX(escript::Data& Line 1218  void MeshAdapter::setToX(escript::Data&
1218        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1219        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1220     } else {     } else {
1221        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1222        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1223        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1224        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1273  void MeshAdapter::setToIntegrals(vector< Line 1312  void MeshAdapter::setToIntegrals(vector<
1312     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1313     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1314     case(Nodes):     case(Nodes):
1315     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1316     _temp=temp.getDataC();     _temp=temp.getDataC();
1317     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1318     break;     break;
1319     case(ReducedNodes):     case(ReducedNodes):
1320     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1321     _temp=temp.getDataC();     _temp=temp.getDataC();
1322     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1323     break;     break;
# Line 1310  void MeshAdapter::setToIntegrals(vector< Line 1349  void MeshAdapter::setToIntegrals(vector<
1349     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1350     break;     break;
1351     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1352     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1353     _temp=temp.getDataC();     _temp=temp.getDataC();
1354     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1355     break;     break;
1356     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1357     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1358     _temp=temp.getDataC();     _temp=temp.getDataC();
1359     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1360     break;     break;
# Line 1347  void MeshAdapter::setToGradient(escript: Line 1386  void MeshAdapter::setToGradient(escript:
1386     escript::Data temp;     escript::Data temp;
1387     if (getMPISize()>1) {     if (getMPISize()>1) {
1388        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1389           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1390           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1391        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1392           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1393           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1394        } else {        } else {
1395           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1469  void MeshAdapter::setNewX(const escript: Line 1508  void MeshAdapter::setNewX(const escript:
1508     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1509     if (newDomain!=*this)     if (newDomain!=*this)
1510        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1511     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1512         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1513         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1514     } else {     } else {
1515         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1516    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1517         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1518         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1519     }     }
1520     checkFinleyError();     checkFinleyError();
1521  }  }
1522    
 //  
 // Helper for the save* methods. Extracts optional data variable names and the  
 // corresponding pointers from python dictionary. Caller must free arrays.  
 //  
 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const  
 {  
    numData = boost::python::extract<int>(arg.attr("__len__")());  
    /* win32 refactor */  
    names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;  
    data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;  
    dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0; i<numData; ++i) {  
       string n=boost::python::extract<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: Data must be defined on same Domain");  
       data[i] = d.getDataC();  
       dataPtr[i] = &(data[i]);  
       names[i] = TMPMEMALLOC(n.length()+1, char);  
       strcpy(names[i], n.c_str());  
    }  
 }  
   
 //  
 // saves mesh and optionally data arrays in openDX format  
 //  
 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
   
    return;  
 }  
   
 //  
 // saves mesh and optionally data arrays in VTK format  
 //  
 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveVTK(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data, metadata.c_str(), metadata_schema.c_str());  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
 }  
   
1523  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1524  {  {
1525  #ifdef PASO_MPI      if (getMPISize() > 1) {
1526      index_t myFirstNode=0, myLastNode=0, k=0;  #ifdef ESYS_MPI
1527      index_t* globalNodeIndex=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1528      Finley_Mesh* mesh_p=m_finleyMesh.get();          index_t* globalNodeIndex=0;
1529      if (fs_code == FINLEY_REDUCED_NODES)          Finley_Mesh* mesh_p=m_finleyMesh.get();
1530      {          /*
1531      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * this method is only used by saveDataCSV which would use the returned
1532      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1533      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES)
1534      }          {
1535      else              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1536      {              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1537      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1538      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          }
1539      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          else
1540      }          */
1541      k=globalNodeIndex[id];          if (fs_code == FINLEY_NODES)
1542      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );          {
1543                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1544                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1545                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1546            }
1547            else
1548            {
1549                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1550            }
1551    
1552            k=globalNodeIndex[id];
1553            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1554  #endif  #endif
1555        }
1556      return true;      return true;
1557  }  }
1558    
1559    
   
1560  //  //
1561  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1562  //  //
1563  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1564                                                   const int row_blocksize,                                                   const int row_blocksize,
1565                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1566                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1599  SystemMatrixAdapter MeshAdapter::newSyst Line 1575  SystemMatrixAdapter MeshAdapter::newSyst
1575        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.");
1576     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1577     if (col_domain!=*this)     if (col_domain!=*this)
1578        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1579     // is the function space type right     // is the function space type right
1580     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1581        reduceRowOrder=0;        reduceRowOrder=0;
# Line 1631  SystemMatrixAdapter MeshAdapter::newSyst Line 1607  SystemMatrixAdapter MeshAdapter::newSyst
1607     }     }
1608     checkPasoError();     checkPasoError();
1609     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1610     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1611       return ASM_ptr(sma);
1612    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1613  }  }
1614    
1615  //  //
1616  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1617  //  //
1618  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const double theta,  
1619                                                           const int blocksize,                                                           const int blocksize,
1620                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1621                                                           const int type) const                                                           const int type) const
# Line 1660  TransportProblemAdapter MeshAdapter::new Line 1637  TransportProblemAdapter MeshAdapter::new
1637    
1638     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1639     checkFinleyError();     checkFinleyError();
1640     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1641     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1642     checkPasoError();     checkPasoError();
1643     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1644     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1645       return ATP_ptr(tpa);
1646    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1647  }  }
1648    
1649  //  //
# Line 1727  MeshAdapter::commonFunctionSpace(const v Line 1706  MeshAdapter::commonFunctionSpace(const v
1706     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.     For classes with multiple members (eg class 2) we have vars to record if there is at least one instance.
1707     eg hasnodes is true if we have at least one instance of Nodes.     eg hasnodes is true if we have at least one instance of Nodes.
1708     */     */
1709      if (fs.size()==0)      if (fs.empty())
1710      {      {
1711          return false;          return false;
1712      }      }
# Line 2000  bool MeshAdapter::probeInterpolationOnDo Line 1979  bool MeshAdapter::probeInterpolationOnDo
1979     return false;     return false;
1980  }  }
1981    
1982    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1983    {
1984        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1985        {
1986            return 1;
1987        }
1988        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1989        {
1990            return -1;
1991        }
1992        return 0;
1993    }  
1994      
1995  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1996  {  {
1997     return false;     return false;
# Line 2023  bool MeshAdapter::operator!=(const Abstr Line 2015  bool MeshAdapter::operator!=(const Abstr
2015  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
2016  {  {
2017     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2018     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2019     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2020  }  }
2021    
2022  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
2023  {  {
2024     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2025     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2026     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2027  }  }
2028    
2029  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2030  {  {
2031     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2032  }  }
2033    
2034  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2035  {  {
2036     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2037  }  }
2038    
2039  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2040  {  {
2041     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2042  }  }
2043    
2044  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2213  void MeshAdapter::setTagMap(const string Line 2204  void MeshAdapter::setTagMap(const string
2204  {  {
2205     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2206     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2207     checkPasoError();     checkFinleyError();
2208     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2209  }  }
2210    
# Line 2222  int MeshAdapter::getTag(const string& na Line 2213  int MeshAdapter::getTag(const string& na
2213     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2214     int tag=0;     int tag=0;
2215     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2216     checkPasoError();     checkFinleyError();
2217     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2218     return tag;     return tag;
2219  }  }
# Line 2395  int MeshAdapter::getApproximationOrder(c Line 2386  int MeshAdapter::getApproximationOrder(c
2386    return order;    return order;
2387  }  }
2388    
2389  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  bool MeshAdapter::supportsContactElements() const
2390    {
2391      return true;
2392    }
2393    
2394    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2395  {  {
2396    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2397  }  }
# Line 2405  ReferenceElementSetWrapper::~ReferenceEl Line 2401  ReferenceElementSetWrapper::~ReferenceEl
2401    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2402  }  }
2403    
2404    // points will be flattened
2405    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2406    {
2407          const int dim = getDim();
2408          int numPoints=points.size()/dim;
2409          int numTags=tags.size();
2410          Finley_Mesh* mesh=m_finleyMesh.get();
2411          
2412          if ( points.size() % dim != 0 )
2413          {
2414        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2415          }
2416          
2417          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2418         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2419          
2420          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2421          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2422          
2423          for (int i=0;i<numPoints;++i) {
2424           points_ptr[ i * dim     ] = points[i * dim ];
2425           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2426           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2427               tags_ptr[i]=tags[i];
2428          }
2429          
2430          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2431          checkFinleyError();
2432          
2433          TMPMEMFREE(points_ptr);
2434          TMPMEMFREE(tags_ptr);
2435    }
2436    
2437    
2438    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2439    // {
2440    //       const int dim = getDim();
2441    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2442    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2443    //       Finley_Mesh* mesh=m_finleyMesh.get();
2444    //      
2445    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2446    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2447    //      
2448    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2449    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2450    //      
2451    //       for (int i=0;i<numPoints;++i) {
2452    //     int tag_id=-1;
2453    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2454    //     if  ( numComps !=   dim ) {
2455    //                stringstream temp;            
2456    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2457    //                throw FinleyAdapterException(temp.str());
2458    //     }
2459    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2460    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2461    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2462    //    
2463    //     if ( numTags > 0) {
2464    //            boost::python::extract<string> ex_str(tags[i]);
2465    //        if  ( ex_str.check() ) {
2466    //            tag_id=getTag( ex_str());
2467    //        } else {
2468    //             boost::python::extract<int> ex_int(tags[i]);
2469    //             if ( ex_int.check() ) {
2470    //                 tag_id=ex_int();
2471    //             } else {
2472    //              stringstream temp;          
2473    //                  temp << "Error - unable to extract tag for point " << i;
2474    //              throw FinleyAdapterException(temp.str());
2475    //            }
2476    //        }
2477    //     }      
2478    //            tags_ptr[i]=tag_id;
2479    //       }
2480    //      
2481    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2482    //       checkPasoError();
2483    //      
2484    //       TMPMEMFREE(points_ptr);
2485    //       TMPMEMFREE(tags_ptr);
2486    // }
2487    
2488    /*
2489    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2490    {  
2491        boost::python::list points =  boost::python::list();
2492        boost::python::list tags =  boost::python::list();
2493        points.append(point);
2494        tags.append(tag);
2495        addDiracPoints(points, tags);
2496    }
2497    */
2498    
2499    /*
2500    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2501    {
2502            boost::python::list points =   boost::python::list();
2503            boost::python::list tags =   boost::python::list();
2504            points.append(point);
2505            tags.append(tag);
2506            addDiracPoints(points, tags);
2507    }
2508    */
2509  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2898  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26