/[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 2842 by caltinay, Thu Jan 14 05:42:02 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-2009 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 220  void MeshAdapter::dump(const string& fil Line 222  void MeshAdapter::dump(const string& fil
222        throw DataException(msgPrefix+"add_att(Name)");        throw DataException(msgPrefix+"add_att(Name)");
223     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
224        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
225     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
226        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
227     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
228        throw DataException(msgPrefix+"add_att(reduced_order)");        throw DataException(msgPrefix+"add_att(reduced_order)");
229     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
230        throw DataException(msgPrefix+"add_att(numNodes)");        throw DataException(msgPrefix+"add_att(numNodes)");
# 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     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
831      
832       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  }  }
839    
840    
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 829  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 838  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 863  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_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     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 878  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 1042  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 1052  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 1062  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 1073  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 1111  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 1121  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 1130  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 1142  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 1179  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 1273  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 1310  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 1347  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 1469  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());  
1521     }     }
 }  
   
 //  
 // 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);  
1522     checkFinleyError();     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  }  }
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 1599  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 1631  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 double theta,  
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 1660  TransportProblemAdapter MeshAdapter::new Line 1639  TransportProblemAdapter MeshAdapter::new
1639    
1640     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1641     checkFinleyError();     checkFinleyError();
1642     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1643     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1644     checkPasoError();     checkPasoError();
1645     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1646     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1647       return ATP_ptr(tpa);
1648    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1649  }  }
1650    
1651  //  //
# Line 1727  MeshAdapter::commonFunctionSpace(const v Line 1708  MeshAdapter::commonFunctionSpace(const v
1708     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.
1709     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.
1710     */     */
1711      if (fs.size()==0)      if (fs.empty())
1712      {      {
1713          return false;          return false;
1714      }      }
# Line 2000  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 2023  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_FCTransportProblem_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 2213  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 2222  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 2360  AbstractDomain::StatusType MeshAdapter:: Line 2353  AbstractDomain::StatusType MeshAdapter::
2353    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2354  }  }
2355    
2356    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2357    {
2358      
2359      Finley_Mesh* mesh=m_finleyMesh.get();
2360      int order =-1;
2361      switch(functionSpaceCode) {
2362       case(Nodes):
2363       case(DegreesOfFreedom):
2364              order=mesh->approximationOrder;
2365              break;
2366       case(ReducedNodes):
2367       case(ReducedDegreesOfFreedom):
2368              order=mesh->reducedApproximationOrder;
2369              break;
2370       case(Elements):
2371       case(FaceElements):
2372       case(Points):
2373       case(ContactElementsZero):
2374       case(ContactElementsOne):
2375              order=mesh->integrationOrder;
2376              break;
2377       case(ReducedElements):
2378       case(ReducedFaceElements):
2379       case(ReducedContactElementsZero):
2380       case(ReducedContactElementsOne):
2381              order=mesh->reducedIntegrationOrder;
2382              break;
2383       default:
2384          stringstream temp;
2385          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2386          throw FinleyAdapterException(temp.str());
2387      }
2388      return order;
2389    }
2390    
2391    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);
2399    }
2400    
2401    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2402    {
2403      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.2842  
changed lines
  Added in v.4286

  ViewVC Help
Powered by ViewVC 1.1.26