/[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 2940 by caltinay, Fri Feb 19 00:38:45 2010 UTC revision 3546 by jfenwick, Thu Jul 28 01:32:20 2011 UTC
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
21  extern "C" {  extern "C" {
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23  }  }
24    
25    #include <boost/python/import.hpp>
26    
27  using namespace std;  using namespace std;
28  using namespace escript;  using namespace escript;
29    
# Line 85  int MeshAdapter::getMPIRank() const Line 83  int MeshAdapter::getMPIRank() const
83  }  }
84  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
85  {  {
86  #ifdef PASO_MPI  #ifdef ESYS_MPI
87     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
88  #endif  #endif
89     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 94  bool MeshAdapter::onMasterProcessor() co
94  }  }
95    
96    
97  #ifdef PASO_MPI  #ifdef ESYS_MPI
98    MPI_Comm    MPI_Comm
99  #else  #else
100    unsigned int    unsigned int
101  #endif  #endif
102  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
103  {  {
104  #ifdef PASO_MPI  #ifdef ESYS_MPI
105      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
106  #else  #else
107      return 0;      return 0;
# Line 149  void MeshAdapter::dump(const string& fil Line 147  void MeshAdapter::dump(const string& fil
147     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
148     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
149     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
150  #ifdef PASO_MPI  #ifdef ESYS_MPI
151     MPI_Status status;     MPI_Status status;
152  #endif  #endif
153    
154  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
155  #ifdef PASO_MPI  #ifdef ESYS_MPI
156     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);
157  #endif  #endif
158    
159     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
160                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
161    
162     /* 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 527  void MeshAdapter::dump(const string& fil
527     }     }
528    
529  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
530  #ifdef PASO_MPI  #ifdef ESYS_MPI
531     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);
532  #endif  #endif
533    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 647  int MeshAdapter::getReducedSolutionCode(
647     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
648  }  }
649    
650  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
651  {  {
652     return Points;     return Points;
653  }  }
# Line 770  pair<int,int> MeshAdapter::getDataShape( Line 768  pair<int,int> MeshAdapter::getDataShape(
768  // 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:
769  //  //
770  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
771                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
772                                   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,
773                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
774                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
775                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
776  {  {
777       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
778       if (smat==0)
779       {
780        throw FinleyAdapterException("finley only supports its own system matrices.");
781       }
782     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
783     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
784     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 786  void MeshAdapter::addPDEToSystem( Line 790  void MeshAdapter::addPDEToSystem(
790     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
791     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
792     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
793       escriptDataC _d_dirac=d_dirac.getDataC();
794       escriptDataC _y_dirac=y_dirac.getDataC();
795    
796     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
797    
798     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 );
799       checkFinleyError();
800    
801       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
802     checkFinleyError();     checkFinleyError();
803    
804     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 );
805     checkFinleyError();     checkFinleyError();
806    
807     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 );
808     checkFinleyError();     checkFinleyError();
809  }  }
810    
811  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
812                                          escript::Data& mat,                                          escript::Data& mat,
813                                          const escript::Data& D,                                          const escript::Data& D,
814                                          const escript::Data& d) const                                          const escript::Data& d,
815                                            const escript::Data& d_dirac,
816                                            const bool useHRZ) const
817  {  {
818     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
819     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
820     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
821       escriptDataC _d_dirac=d_dirac.getDataC();
822    
823     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
824    
825     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
826     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
827      
828       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
829       checkFinleyError();
830    
831       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
832     checkFinleyError();     checkFinleyError();
833    
834  }  }
835    
836    
837  //  //
838  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
839  //  //
840  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
841  {  {
842     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
843    
# Line 829  void MeshAdapter::addPDEToRHS( escript:: Line 846  void MeshAdapter::addPDEToRHS( escript::
846     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
847     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
848     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
849       escriptDataC _y_dirac=y_dirac.getDataC();
850    
851     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 );
852     checkFinleyError();     checkFinleyError();
# Line 838  void MeshAdapter::addPDEToRHS( escript:: Line 856  void MeshAdapter::addPDEToRHS( escript::
856    
857     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 );
858     checkFinleyError();     checkFinleyError();
859    
860       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
861       checkFinleyError();
862    
863  }  }
864  //  //
865  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
866  //  //
867  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
868                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
869                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
870                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
871                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
872                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
873                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
874  {  {
875       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
876       if (tpa==0)
877       {
878        throw FinleyAdapterException("finley only supports its own transport problems.");
879       }
880    
881    
882     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
883     source.expand();     source.expand();
884     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 863  void MeshAdapter::addPDEToTransportProbl Line 893  void MeshAdapter::addPDEToTransportProbl
893     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
894     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
895     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
896       escriptDataC _d_dirac=d_dirac.getDataC();
897       escriptDataC _y_dirac=y_dirac.getDataC();
898    
899    
900     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
901     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
902    
903     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 );
904     checkFinleyError();     checkFinleyError();
# Line 878  void MeshAdapter::addPDEToTransportProbl Line 911  void MeshAdapter::addPDEToTransportProbl
911    
912     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 );
913     checkFinleyError();     checkFinleyError();
914    
915       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
916       checkFinleyError();
917    
918  }  }
919    
920  //  //
# Line 1042  void MeshAdapter::interpolateOnDomain(es Line 1079  void MeshAdapter::interpolateOnDomain(es
1079        case(Elements):        case(Elements):
1080        case(ReducedElements):        case(ReducedElements):
1081        if (getMPISize()>1) {        if (getMPISize()>1) {
1082           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1083           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1084           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1085        } else {        } else {
# Line 1052  void MeshAdapter::interpolateOnDomain(es Line 1089  void MeshAdapter::interpolateOnDomain(es
1089        case(FaceElements):        case(FaceElements):
1090        case(ReducedFaceElements):        case(ReducedFaceElements):
1091        if (getMPISize()>1) {        if (getMPISize()>1) {
1092           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1093           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1094           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1095        
# Line 1062  void MeshAdapter::interpolateOnDomain(es Line 1099  void MeshAdapter::interpolateOnDomain(es
1099        break;        break;
1100        case(Points):        case(Points):
1101        if (getMPISize()>1) {        if (getMPISize()>1) {
1102           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1103           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1104        } else {        } else {
1105           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1106        }        }
# Line 1073  void MeshAdapter::interpolateOnDomain(es Line 1110  void MeshAdapter::interpolateOnDomain(es
1110        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1111        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1112        if (getMPISize()>1) {        if (getMPISize()>1) {
1113           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1114           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1115           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1116        } else {        } else {
# Line 1111  void MeshAdapter::interpolateOnDomain(es Line 1148  void MeshAdapter::interpolateOnDomain(es
1148        case(Elements):        case(Elements):
1149        case(ReducedElements):        case(ReducedElements):
1150        if (getMPISize()>1) {        if (getMPISize()>1) {
1151           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1152           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1153           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1154        } else {        } else {
# Line 1121  void MeshAdapter::interpolateOnDomain(es Line 1158  void MeshAdapter::interpolateOnDomain(es
1158        case(FaceElements):        case(FaceElements):
1159        case(ReducedFaceElements):        case(ReducedFaceElements):
1160        if (getMPISize()>1) {        if (getMPISize()>1) {
1161           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1162           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1163           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1164        } else {        } else {
# Line 1130  void MeshAdapter::interpolateOnDomain(es Line 1167  void MeshAdapter::interpolateOnDomain(es
1167        break;        break;
1168        case(Points):        case(Points):
1169        if (getMPISize()>1) {        if (getMPISize()>1) {
1170           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1171           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1172           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1173        } else {        } else {
# Line 1142  void MeshAdapter::interpolateOnDomain(es Line 1179  void MeshAdapter::interpolateOnDomain(es
1179        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1180        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1181        if (getMPISize()>1) {        if (getMPISize()>1) {
1182           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1183           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1184           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1185        } else {        } else {
# Line 1179  void MeshAdapter::setToX(escript::Data& Line 1216  void MeshAdapter::setToX(escript::Data&
1216        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1217        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1218     } else {     } else {
1219        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1220        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1222        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1273  void MeshAdapter::setToIntegrals(vector< Line 1310  void MeshAdapter::setToIntegrals(vector<
1310     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1311     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1312     case(Nodes):     case(Nodes):
1313     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1314     _temp=temp.getDataC();     _temp=temp.getDataC();
1315     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1316     break;     break;
1317     case(ReducedNodes):     case(ReducedNodes):
1318     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1319     _temp=temp.getDataC();     _temp=temp.getDataC();
1320     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1321     break;     break;
# Line 1310  void MeshAdapter::setToIntegrals(vector< Line 1347  void MeshAdapter::setToIntegrals(vector<
1347     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1348     break;     break;
1349     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1350     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1351     _temp=temp.getDataC();     _temp=temp.getDataC();
1352     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1353     break;     break;
1354     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1355     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1356     _temp=temp.getDataC();     _temp=temp.getDataC();
1357     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1358     break;     break;
# Line 1347  void MeshAdapter::setToGradient(escript: Line 1384  void MeshAdapter::setToGradient(escript:
1384     escript::Data temp;     escript::Data temp;
1385     if (getMPISize()>1) {     if (getMPISize()>1) {
1386        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1387           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1388           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1389        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1390           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1391           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1392        } else {        } else {
1393           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1469  void MeshAdapter::setNewX(const escript: Line 1506  void MeshAdapter::setNewX(const escript:
1506     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1507     if (newDomain!=*this)     if (newDomain!=*this)
1508        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1509     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1510         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1511         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1512     } else {     } else {
1513         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1514         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1515         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1516     }     }
# Line 1536  void MeshAdapter::saveDX(const string& f Line 1573  void MeshAdapter::saveDX(const string& f
1573  //  //
1574  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const
1575  {  {
1576     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1577     char **names;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1578     escriptDataC *data;                metadata, metadata_schema, arg);
    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);  
1579  }  }
1580    
1581  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1582  {  {
1583  #ifdef PASO_MPI  #ifdef ESYS_MPI
1584      index_t myFirstNode=0, myLastNode=0, k=0;      index_t myFirstNode=0, myLastNode=0, k=0;
1585      index_t* globalNodeIndex=0;      index_t* globalNodeIndex=0;
1586      Finley_Mesh* mesh_p=m_finleyMesh.get();      Finley_Mesh* mesh_p=m_finleyMesh.get();
# Line 1584  bool MeshAdapter::ownSample(int fs_code, Line 1607  bool MeshAdapter::ownSample(int fs_code,
1607  //  //
1608  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1609  //  //
1610  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1611                                                   const int row_blocksize,                                                   const int row_blocksize,
1612                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1613                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1631  SystemMatrixAdapter MeshAdapter::newSyst Line 1654  SystemMatrixAdapter MeshAdapter::newSyst
1654     }     }
1655     checkPasoError();     checkPasoError();
1656     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1657     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1658       return ASM_ptr(sma);
1659    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1660  }  }
1661    
1662  //  //
1663  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1664  //  //
1665  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1666                                                           const double theta,                                                           const bool useBackwardEuler,
1667                                                           const int blocksize,                                                           const int blocksize,
1668                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1669                                                           const int type) const                                                           const int type) const
# Line 1660  TransportProblemAdapter MeshAdapter::new Line 1685  TransportProblemAdapter MeshAdapter::new
1685    
1686     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1687     checkFinleyError();     checkFinleyError();
1688     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1689     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1690     checkPasoError();     checkPasoError();
1691     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1692     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1693       return ATP_ptr(tpa);
1694    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1695  }  }
1696    
1697  //  //
# Line 2030  int MeshAdapter::getSystemMatrixTypeId(c Line 2057  int MeshAdapter::getSystemMatrixTypeId(c
2057  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
2058  {  {
2059     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2060     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);
2061     checkPasoError();     checkPasoError();
2062     return out;     return out;
2063  }  }
2064    
2065  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2066  {  {
2067     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2068  }  }
2069    
2070  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2071  {  {
2072     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2073  }  }
2074    
2075  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2076  {  {
2077     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2078  }  }
2079    
2080  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2395  int MeshAdapter::getApproximationOrder(c Line 2422  int MeshAdapter::getApproximationOrder(c
2422    return order;    return order;
2423  }  }
2424    
2425  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  bool MeshAdapter::supportsContactElements() const
2426    {
2427      return true;
2428    }
2429    
2430    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2431  {  {
2432    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2433  }  }
# Line 2405  ReferenceElementSetWrapper::~ReferenceEl Line 2437  ReferenceElementSetWrapper::~ReferenceEl
2437    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
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    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2491    {  
2492        boost::python::list points =  boost::python::list();
2493        boost::python::list tags =  boost::python::list();
2494        points.append(point);
2495        tags.append(tag);
2496        addDiracPoints(points, tags);
2497    }
2498    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2499    {
2500            boost::python::list points =   boost::python::list();
2501            boost::python::list tags =   boost::python::list();
2502            points.append(point);
2503            tags.append(tag);
2504            addDiracPoints(points, tags);
2505    }  
2506  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2940  
changed lines
  Added in v.3546

  ViewVC Help
Powered by ViewVC 1.1.26