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

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

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

revision 2989 by gross, Thu Mar 18 01:45:55 2010 UTC revision 4286 by caltinay, Thu Mar 7 04:28:11 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  
24  extern "C" {  extern "C" {
25  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
26  }  }
27    
28    #include <boost/python/import.hpp>
29    
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    using namespace paso;
33    
34  namespace finley {  namespace finley {
35    
# Line 85  int MeshAdapter::getMPIRank() const Line 87  int MeshAdapter::getMPIRank() const
87  }  }
88  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
89  {  {
90  #ifdef PASO_MPI  #ifdef ESYS_MPI
91     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
92  #endif  #endif
93     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 98  bool MeshAdapter::onMasterProcessor() co
98  }  }
99    
100    
101  #ifdef PASO_MPI  #ifdef ESYS_MPI
102    MPI_Comm    MPI_Comm
103  #else  #else
104    unsigned int    unsigned int
105  #endif  #endif
106  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
107  {  {
108  #ifdef PASO_MPI  #ifdef ESYS_MPI
109      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
110  #else  #else
111      return 0;      return 0;
# Line 149  void MeshAdapter::dump(const string& fil Line 151  void MeshAdapter::dump(const string& fil
151     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
152     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
153     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
154  #ifdef PASO_MPI  #ifdef ESYS_MPI
155     MPI_Status status;     MPI_Status status;
156  #endif  #endif
157    
158  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
159  #ifdef PASO_MPI  #ifdef ESYS_MPI
160     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);
161  #endif  #endif
162    
163     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
164                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
165    
166     /* 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 531  void MeshAdapter::dump(const string& fil
531     }     }
532    
533  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
534  #ifdef PASO_MPI  #ifdef ESYS_MPI
535     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);
536  #endif  #endif
537    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 651  int MeshAdapter::getReducedSolutionCode(
651     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
652  }  }
653    
654  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
655  {  {
656     return Points;     return Points;
657  }  }
# Line 770  pair<int,int> MeshAdapter::getDataShape( Line 772  pair<int,int> MeshAdapter::getDataShape(
772  // 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:
773  //  //
774  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
775                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
776                                   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,
777                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
778                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
779                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
780  {  {
781       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782       if (smat==0)
783       {
784        throw FinleyAdapterException("finley only supports Paso system matrices.");
785       }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
788     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 786  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
795     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
796     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
797       escriptDataC _d_dirac=d_dirac.getDataC();
798       escriptDataC _y_dirac=y_dirac.getDataC();
799    
800     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
801    
802     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 );
803       checkFinleyError();
804    
805       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
806     checkFinleyError();     checkFinleyError();
807    
808     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 );
809     checkFinleyError();     checkFinleyError();
810    
811     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 );
812     checkFinleyError();     checkFinleyError();
813  }  }
814    
815  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
816                                          escript::Data& mat,                                          escript::Data& mat,
817                                          const escript::Data& D,                                          const escript::Data& D,
818                                          const escript::Data& d) const                                          const escript::Data& d,
819                                            const escript::Data& d_dirac,
820                                            const bool useHRZ) const
821  {  {
822     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
823     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
824     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
825       escriptDataC _d_dirac=d_dirac.getDataC();
826    
827     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
828    
829     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
830     checkFinleyError();     checkFinleyError();
831        
832     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
833       checkFinleyError();
834    
835       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
836     checkFinleyError();     checkFinleyError();
837    
838  }  }
# Line 822  void  MeshAdapter::addPDEToLumpedSystem( Line 841  void  MeshAdapter::addPDEToLumpedSystem(
841  //  //
842  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
843  //  //
844  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
845  {  {
846     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
847    
# Line 831  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
851     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
852     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
853       escriptDataC _y_dirac=y_dirac.getDataC();
854    
855     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 );
856     checkFinleyError();     checkFinleyError();
# Line 840  void MeshAdapter::addPDEToRHS( escript:: Line 860  void MeshAdapter::addPDEToRHS( escript::
860    
861     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 );
862     checkFinleyError();     checkFinleyError();
863    
864       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
865       checkFinleyError();
866    
867  }  }
868  //  //
869  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
870  //  //
871  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
872                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
873                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
874                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
875                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
876                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
877                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
878  {  {
879       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
880       if (tpa==0)
881       {
882        throw FinleyAdapterException("finley only supports Paso transport problems.");
883       }
884    
885    
886     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
887     source.expand();     source.expand();
888     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 865  void MeshAdapter::addPDEToTransportProbl Line 897  void MeshAdapter::addPDEToTransportProbl
897     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
898     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
899     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
900       escriptDataC _d_dirac=d_dirac.getDataC();
901       escriptDataC _y_dirac=y_dirac.getDataC();
902    
903    
904     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
905     Paso_TransportProblem* _tp = tp.getPaso_TransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
906    
907     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 );
908     checkFinleyError();     checkFinleyError();
# Line 880  void MeshAdapter::addPDEToTransportProbl Line 915  void MeshAdapter::addPDEToTransportProbl
915    
916     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 );
917     checkFinleyError();     checkFinleyError();
918    
919       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
920       checkFinleyError();
921    
922  }  }
923    
924  //  //
# Line 1044  void MeshAdapter::interpolateOnDomain(es Line 1083  void MeshAdapter::interpolateOnDomain(es
1083        case(Elements):        case(Elements):
1084        case(ReducedElements):        case(ReducedElements):
1085        if (getMPISize()>1) {        if (getMPISize()>1) {
1086           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1088           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1089        } else {        } else {
# Line 1054  void MeshAdapter::interpolateOnDomain(es Line 1093  void MeshAdapter::interpolateOnDomain(es
1093        case(FaceElements):        case(FaceElements):
1094        case(ReducedFaceElements):        case(ReducedFaceElements):
1095        if (getMPISize()>1) {        if (getMPISize()>1) {
1096           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1097           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1098           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1099        
# Line 1064  void MeshAdapter::interpolateOnDomain(es Line 1103  void MeshAdapter::interpolateOnDomain(es
1103        break;        break;
1104        case(Points):        case(Points):
1105        if (getMPISize()>1) {        if (getMPISize()>1) {
1106           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1107           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1108        } else {        } else {
1109           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1110        }        }
# Line 1075  void MeshAdapter::interpolateOnDomain(es Line 1114  void MeshAdapter::interpolateOnDomain(es
1114        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1115        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1116        if (getMPISize()>1) {        if (getMPISize()>1) {
1117           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1118           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1119           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120        } else {        } else {
# Line 1113  void MeshAdapter::interpolateOnDomain(es Line 1152  void MeshAdapter::interpolateOnDomain(es
1152        case(Elements):        case(Elements):
1153        case(ReducedElements):        case(ReducedElements):
1154        if (getMPISize()>1) {        if (getMPISize()>1) {
1155           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1157           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1158        } else {        } else {
# Line 1123  void MeshAdapter::interpolateOnDomain(es Line 1162  void MeshAdapter::interpolateOnDomain(es
1162        case(FaceElements):        case(FaceElements):
1163        case(ReducedFaceElements):        case(ReducedFaceElements):
1164        if (getMPISize()>1) {        if (getMPISize()>1) {
1165           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1166           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1167           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1168        } else {        } else {
# Line 1132  void MeshAdapter::interpolateOnDomain(es Line 1171  void MeshAdapter::interpolateOnDomain(es
1171        break;        break;
1172        case(Points):        case(Points):
1173        if (getMPISize()>1) {        if (getMPISize()>1) {
1174           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1175           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1176           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1177        } else {        } else {
# Line 1144  void MeshAdapter::interpolateOnDomain(es Line 1183  void MeshAdapter::interpolateOnDomain(es
1183        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1184        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1185        if (getMPISize()>1) {        if (getMPISize()>1) {
1186           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1187           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1188           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1189        } else {        } else {
# Line 1181  void MeshAdapter::setToX(escript::Data& Line 1220  void MeshAdapter::setToX(escript::Data&
1220        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1222     } else {     } else {
1223        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1224        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1225        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1226        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1275  void MeshAdapter::setToIntegrals(vector< Line 1314  void MeshAdapter::setToIntegrals(vector<
1314     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1315     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1316     case(Nodes):     case(Nodes):
1317     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1318     _temp=temp.getDataC();     _temp=temp.getDataC();
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320     break;     break;
1321     case(ReducedNodes):     case(ReducedNodes):
1322     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1323     _temp=temp.getDataC();     _temp=temp.getDataC();
1324     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1325     break;     break;
# Line 1312  void MeshAdapter::setToIntegrals(vector< Line 1351  void MeshAdapter::setToIntegrals(vector<
1351     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1352     break;     break;
1353     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1354     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1355     _temp=temp.getDataC();     _temp=temp.getDataC();
1356     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1357     break;     break;
1358     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1359     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1360     _temp=temp.getDataC();     _temp=temp.getDataC();
1361     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1362     break;     break;
# Line 1349  void MeshAdapter::setToGradient(escript: Line 1388  void MeshAdapter::setToGradient(escript:
1388     escript::Data temp;     escript::Data temp;
1389     if (getMPISize()>1) {     if (getMPISize()>1) {
1390        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1391           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1392           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1393        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1394           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1395           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1396        } else {        } else {
1397           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1471  void MeshAdapter::setNewX(const escript: Line 1510  void MeshAdapter::setNewX(const escript:
1510     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1511     if (newDomain!=*this)     if (newDomain!=*this)
1512        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1513     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1514         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1515         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1516     } else {     } else {
1517         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.");    
1518    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1520         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
    }  
    checkFinleyError();  
 }  
   
 //  
 // 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]);  
1521     }     }
    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());  
1522     checkFinleyError();     checkFinleyError();
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1523  }  }
1524    
1525  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1526  {  {
1527  #ifdef PASO_MPI      if (getMPISize() > 1) {
1528      index_t myFirstNode=0, myLastNode=0, k=0;  #ifdef ESYS_MPI
1529      index_t* globalNodeIndex=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1530      Finley_Mesh* mesh_p=m_finleyMesh.get();          index_t* globalNodeIndex=0;
1531      if (fs_code == FINLEY_REDUCED_NODES)          Finley_Mesh* mesh_p=m_finleyMesh.get();
1532      {          /*
1533      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * this method is only used by saveDataCSV which would use the returned
1534      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1535      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES)
1536      }          {
1537      else              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1538      {              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1539      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1540      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          }
1541      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          else
1542      }          */
1543      k=globalNodeIndex[id];          if (fs_code == FINLEY_NODES)
1544      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );          {
1545                myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1546                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1547                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1548            }
1549            else
1550            {
1551                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1552            }
1553    
1554            k=globalNodeIndex[id];
1555            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1556  #endif  #endif
1557        }
1558      return true;      return true;
1559  }  }
1560    
1561    
   
1562  //  //
1563  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1564  //  //
1565  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1566                                                   const int row_blocksize,                                                   const int row_blocksize,
1567                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1568                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1601  SystemMatrixAdapter MeshAdapter::newSyst Line 1577  SystemMatrixAdapter MeshAdapter::newSyst
1577        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.");
1578     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1579     if (col_domain!=*this)     if (col_domain!=*this)
1580        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.");
1581     // is the function space type right     // is the function space type right
1582     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1583        reduceRowOrder=0;        reduceRowOrder=0;
# Line 1633  SystemMatrixAdapter MeshAdapter::newSyst Line 1609  SystemMatrixAdapter MeshAdapter::newSyst
1609     }     }
1610     checkPasoError();     checkPasoError();
1611     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1612     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1613       return ASM_ptr(sma);
1614    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1615  }  }
1616    
1617  //  //
1618  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1619  //  //
1620  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const bool useBackwardEuler,  
1621                                                           const int blocksize,                                                           const int blocksize,
1622                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1623                                                           const int type) const                                                           const int type) const
# Line 1663  TransportProblemAdapter MeshAdapter::new Line 1640  TransportProblemAdapter MeshAdapter::new
1640     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1641     checkFinleyError();     checkFinleyError();
1642     Paso_TransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1643     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1644     checkPasoError();     checkPasoError();
1645     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1646     return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1647       return ATP_ptr(tpa);
1648    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1649  }  }
1650    
1651  //  //
# Line 2002  bool MeshAdapter::probeInterpolationOnDo Line 1981  bool MeshAdapter::probeInterpolationOnDo
1981     return false;     return false;
1982  }  }
1983    
1984    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1985    {
1986        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1987        {
1988            return 1;
1989        }
1990        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1991        {
1992            return -1;
1993        }
1994        return 0;
1995    }  
1996      
1997  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
1998  {  {
1999     return false;     return false;
# Line 2025  bool MeshAdapter::operator!=(const Abstr Line 2017  bool MeshAdapter::operator!=(const Abstr
2017  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
2018  {  {
2019     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2020     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2021     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2022  }  }
2023    
2024  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
2025  {  {
2026     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2027     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2028     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2029  }  }
2030    
2031  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2032  {  {
2033     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2034  }  }
2035    
2036  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2037  {  {
2038     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2039  }  }
2040    
2041  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2042  {  {
2043     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2044  }  }
2045    
2046  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2215  void MeshAdapter::setTagMap(const string Line 2206  void MeshAdapter::setTagMap(const string
2206  {  {
2207     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2208     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2209     checkPasoError();     checkFinleyError();
2210     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2211  }  }
2212    
# Line 2224  int MeshAdapter::getTag(const string& na Line 2215  int MeshAdapter::getTag(const string& na
2215     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2216     int tag=0;     int tag=0;
2217     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2218     checkPasoError();     checkFinleyError();
2219     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2220     return tag;     return tag;
2221  }  }
# Line 2397  int MeshAdapter::getApproximationOrder(c Line 2388  int MeshAdapter::getApproximationOrder(c
2388    return order;    return order;
2389  }  }
2390    
2391  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  bool MeshAdapter::supportsContactElements() const
2392    {
2393      return true;
2394    }
2395    
2396    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2397  {  {
2398    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2399  }  }
# Line 2407  ReferenceElementSetWrapper::~ReferenceEl Line 2403  ReferenceElementSetWrapper::~ReferenceEl
2403    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2404  }  }
2405    
2406    // points will be flattened
2407    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2408    {
2409          const int dim = getDim();
2410          int numPoints=points.size()/dim;
2411          int numTags=tags.size();
2412          Finley_Mesh* mesh=m_finleyMesh.get();
2413          
2414          if ( points.size() % dim != 0 )
2415          {
2416        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2417          }
2418          
2419          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2420         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2421          
2422          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2423          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2424          
2425          for (int i=0;i<numPoints;++i) {
2426           points_ptr[ i * dim     ] = points[i * dim ];
2427           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2428           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2429               tags_ptr[i]=tags[i];
2430          }
2431          
2432          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2433          checkFinleyError();
2434          
2435          TMPMEMFREE(points_ptr);
2436          TMPMEMFREE(tags_ptr);
2437    }
2438    
2439    
2440    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2441    // {
2442    //       const int dim = getDim();
2443    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2444    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2445    //       Finley_Mesh* mesh=m_finleyMesh.get();
2446    //      
2447    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2448    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2449    //      
2450    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2451    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2452    //      
2453    //       for (int i=0;i<numPoints;++i) {
2454    //     int tag_id=-1;
2455    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2456    //     if  ( numComps !=   dim ) {
2457    //                stringstream temp;            
2458    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2459    //                throw FinleyAdapterException(temp.str());
2460    //     }
2461    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2462    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2463    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2464    //    
2465    //     if ( numTags > 0) {
2466    //            boost::python::extract<string> ex_str(tags[i]);
2467    //        if  ( ex_str.check() ) {
2468    //            tag_id=getTag( ex_str());
2469    //        } else {
2470    //             boost::python::extract<int> ex_int(tags[i]);
2471    //             if ( ex_int.check() ) {
2472    //                 tag_id=ex_int();
2473    //             } else {
2474    //              stringstream temp;          
2475    //                  temp << "Error - unable to extract tag for point " << i;
2476    //              throw FinleyAdapterException(temp.str());
2477    //            }
2478    //        }
2479    //     }      
2480    //            tags_ptr[i]=tag_id;
2481    //       }
2482    //      
2483    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2484    //       checkPasoError();
2485    //      
2486    //       TMPMEMFREE(points_ptr);
2487    //       TMPMEMFREE(tags_ptr);
2488    // }
2489    
2490    /*
2491    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2492    {  
2493        boost::python::list points =  boost::python::list();
2494        boost::python::list tags =  boost::python::list();
2495        points.append(point);
2496        tags.append(tag);
2497        addDiracPoints(points, tags);
2498    }
2499    */
2500    
2501    /*
2502    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2503    {
2504            boost::python::list points =   boost::python::list();
2505            boost::python::list tags =   boost::python::list();
2506            points.append(point);
2507            tags.append(tag);
2508            addDiracPoints(points, tags);
2509    }
2510    */
2511  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2989  
changed lines
  Added in v.4286

  ViewVC Help
Powered by ViewVC 1.1.26