/[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 3490 by caltinay, Wed Mar 30 02:24:33 2011 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2010 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# 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 220  void MeshAdapter::dump(const string& fil Line 218  void MeshAdapter::dump(const string& fil
218        throw DataException(msgPrefix+"add_att(Name)");        throw DataException(msgPrefix+"add_att(Name)");
219     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
220        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
221     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
222        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
223     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
224        throw DataException(msgPrefix+"add_att(reduced_order)");        throw DataException(msgPrefix+"add_att(reduced_order)");
225     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
226        throw DataException(msgPrefix+"add_att(numNodes)");        throw DataException(msgPrefix+"add_att(numNodes)");
# 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 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) const
775  {  {
776       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
777       if (smat==0)
778       {
779        throw FinleyAdapterException("finley only supports its own system matrices.");
780       }
781     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
782     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
783     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 789  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792    
793     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
794    
795     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 );
796     checkFinleyError();     checkFinleyError();
797    
798     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
799     checkFinleyError();     checkFinleyError();
800    
801     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->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
802     checkFinleyError();     checkFinleyError();
803  }  }
804    
805  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
806                                          escript::Data& mat,                                          escript::Data& mat,
807                                          const escript::Data& D,                                          const escript::Data& D,
808                                          const escript::Data& d) const                                          const escript::Data& d,
809                                            const bool useHRZ) const
810  {  {
811     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
812     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
# Line 810  void  MeshAdapter::addPDEToLumpedSystem( Line 814  void  MeshAdapter::addPDEToLumpedSystem(
814    
815     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
816    
817     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
818     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
819      
820       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
821     checkFinleyError();     checkFinleyError();
822    
823  }  }
824    
825    
# Line 843  void MeshAdapter::addPDEToRHS( escript:: Line 849  void MeshAdapter::addPDEToRHS( escript::
849  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
850  //  //
851  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
852                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
853                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
854                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
855                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
856                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
857  {  {
858       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
859       if (tpa==0)
860       {
861        throw FinleyAdapterException("finley only supports its own transport problems.");
862       }
863    
864    
865     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
866     source.expand();     source.expand();
867     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 865  void MeshAdapter::addPDEToTransportProbl Line 878  void MeshAdapter::addPDEToTransportProbl
878     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
879    
880     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
881     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
882    
883     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 );
884     checkFinleyError();     checkFinleyError();
# Line 1042  void MeshAdapter::interpolateOnDomain(es Line 1055  void MeshAdapter::interpolateOnDomain(es
1055        case(Elements):        case(Elements):
1056        case(ReducedElements):        case(ReducedElements):
1057        if (getMPISize()>1) {        if (getMPISize()>1) {
1058           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1059           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1060           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1061        } else {        } else {
# Line 1052  void MeshAdapter::interpolateOnDomain(es Line 1065  void MeshAdapter::interpolateOnDomain(es
1065        case(FaceElements):        case(FaceElements):
1066        case(ReducedFaceElements):        case(ReducedFaceElements):
1067        if (getMPISize()>1) {        if (getMPISize()>1) {
1068           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1069           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1070           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1071        
# Line 1062  void MeshAdapter::interpolateOnDomain(es Line 1075  void MeshAdapter::interpolateOnDomain(es
1075        break;        break;
1076        case(Points):        case(Points):
1077        if (getMPISize()>1) {        if (getMPISize()>1) {
1078           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1079           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1080        } else {        } else {
1081           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1082        }        }
# Line 1073  void MeshAdapter::interpolateOnDomain(es Line 1086  void MeshAdapter::interpolateOnDomain(es
1086        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1087        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1088        if (getMPISize()>1) {        if (getMPISize()>1) {
1089           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1090           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1091           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1092        } else {        } else {
# Line 1111  void MeshAdapter::interpolateOnDomain(es Line 1124  void MeshAdapter::interpolateOnDomain(es
1124        case(Elements):        case(Elements):
1125        case(ReducedElements):        case(ReducedElements):
1126        if (getMPISize()>1) {        if (getMPISize()>1) {
1127           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1128           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1129           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1130        } else {        } else {
# Line 1121  void MeshAdapter::interpolateOnDomain(es Line 1134  void MeshAdapter::interpolateOnDomain(es
1134        case(FaceElements):        case(FaceElements):
1135        case(ReducedFaceElements):        case(ReducedFaceElements):
1136        if (getMPISize()>1) {        if (getMPISize()>1) {
1137           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1138           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1139           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1140        } else {        } else {
# Line 1130  void MeshAdapter::interpolateOnDomain(es Line 1143  void MeshAdapter::interpolateOnDomain(es
1143        break;        break;
1144        case(Points):        case(Points):
1145        if (getMPISize()>1) {        if (getMPISize()>1) {
1146           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1147           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1148           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1149        } else {        } else {
# Line 1142  void MeshAdapter::interpolateOnDomain(es Line 1155  void MeshAdapter::interpolateOnDomain(es
1155        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1156        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1157        if (getMPISize()>1) {        if (getMPISize()>1) {
1158           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1159           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1160           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1161        } else {        } else {
# Line 1179  void MeshAdapter::setToX(escript::Data& Line 1192  void MeshAdapter::setToX(escript::Data&
1192        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1193        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1194     } else {     } else {
1195        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1196        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1197        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1198        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1273  void MeshAdapter::setToIntegrals(vector< Line 1286  void MeshAdapter::setToIntegrals(vector<
1286     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1287     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1288     case(Nodes):     case(Nodes):
1289     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1290     _temp=temp.getDataC();     _temp=temp.getDataC();
1291     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1292     break;     break;
1293     case(ReducedNodes):     case(ReducedNodes):
1294     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1295     _temp=temp.getDataC();     _temp=temp.getDataC();
1296     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1297     break;     break;
# Line 1310  void MeshAdapter::setToIntegrals(vector< Line 1323  void MeshAdapter::setToIntegrals(vector<
1323     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1324     break;     break;
1325     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1326     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1327     _temp=temp.getDataC();     _temp=temp.getDataC();
1328     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1329     break;     break;
1330     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1331     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1332     _temp=temp.getDataC();     _temp=temp.getDataC();
1333     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1334     break;     break;
# Line 1347  void MeshAdapter::setToGradient(escript: Line 1360  void MeshAdapter::setToGradient(escript:
1360     escript::Data temp;     escript::Data temp;
1361     if (getMPISize()>1) {     if (getMPISize()>1) {
1362        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1363           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1364           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1365        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1366           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1367           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1368        } else {        } else {
1369           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1469  void MeshAdapter::setNewX(const escript: Line 1482  void MeshAdapter::setNewX(const escript:
1482     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1483     if (newDomain!=*this)     if (newDomain!=*this)
1484        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1485     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1486         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1487         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1488     } else {     } else {
1489         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1490         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1491         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1492     }     }
# Line 1536  void MeshAdapter::saveDX(const string& f Line 1549  void MeshAdapter::saveDX(const string& f
1549  //  //
1550  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
1551  {  {
1552     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1553     char **names;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1554     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);  
1555  }  }
1556    
1557  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1558  {  {
1559  #ifdef PASO_MPI  #ifdef ESYS_MPI
1560      index_t myFirstNode=0, myLastNode=0, k=0;      index_t myFirstNode=0, myLastNode=0, k=0;
1561      index_t* globalNodeIndex=0;      index_t* globalNodeIndex=0;
1562      Finley_Mesh* mesh_p=m_finleyMesh.get();      Finley_Mesh* mesh_p=m_finleyMesh.get();
# Line 1584  bool MeshAdapter::ownSample(int fs_code, Line 1583  bool MeshAdapter::ownSample(int fs_code,
1583  //  //
1584  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1585  //  //
1586  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1587                                                   const int row_blocksize,                                                   const int row_blocksize,
1588                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1589                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1631  SystemMatrixAdapter MeshAdapter::newSyst Line 1630  SystemMatrixAdapter MeshAdapter::newSyst
1630     }     }
1631     checkPasoError();     checkPasoError();
1632     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1633     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1634       return ASM_ptr(sma);
1635    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1636  }  }
1637    
1638  //  //
1639  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1640  //  //
1641  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1642                                                           const double theta,                                                           const bool useBackwardEuler,
1643                                                           const int blocksize,                                                           const int blocksize,
1644                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1645                                                           const int type) const                                                           const int type) const
# Line 1660  TransportProblemAdapter MeshAdapter::new Line 1661  TransportProblemAdapter MeshAdapter::new
1661    
1662     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1663     checkFinleyError();     checkFinleyError();
1664     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1665     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1666     checkPasoError();     checkPasoError();
1667     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1668     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1669       return ATP_ptr(tpa);
1670    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1671  }  }
1672    
1673  //  //
# Line 1727  MeshAdapter::commonFunctionSpace(const v Line 1730  MeshAdapter::commonFunctionSpace(const v
1730     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.
1731     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.
1732     */     */
1733      if (fs.size()==0)      if (fs.empty())
1734      {      {
1735          return false;          return false;
1736      }      }
# Line 2030  int MeshAdapter::getSystemMatrixTypeId(c Line 2033  int MeshAdapter::getSystemMatrixTypeId(c
2033  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
2034  {  {
2035     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2036     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);
2037     checkPasoError();     checkPasoError();
2038     return out;     return out;
2039  }  }
2040    
2041  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2042  {  {
2043     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2044  }  }
2045    
2046  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2047  {  {
2048     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2049  }  }
2050    
2051  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2052  {  {
2053     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2054  }  }
2055    
2056  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2360  AbstractDomain::StatusType MeshAdapter:: Line 2363  AbstractDomain::StatusType MeshAdapter::
2363    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2364  }  }
2365    
2366    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2367    {
2368      
2369      Finley_Mesh* mesh=m_finleyMesh.get();
2370      int order =-1;
2371      switch(functionSpaceCode) {
2372       case(Nodes):
2373       case(DegreesOfFreedom):
2374              order=mesh->approximationOrder;
2375              break;
2376       case(ReducedNodes):
2377       case(ReducedDegreesOfFreedom):
2378              order=mesh->reducedApproximationOrder;
2379              break;
2380       case(Elements):
2381       case(FaceElements):
2382       case(Points):
2383       case(ContactElementsZero):
2384       case(ContactElementsOne):
2385              order=mesh->integrationOrder;
2386              break;
2387       case(ReducedElements):
2388       case(ReducedFaceElements):
2389       case(ReducedContactElementsZero):
2390       case(ReducedContactElementsOne):
2391              order=mesh->reducedIntegrationOrder;
2392              break;
2393       default:
2394          stringstream temp;
2395          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2396          throw FinleyAdapterException(temp.str());
2397      }
2398      return order;
2399    }
2400    
2401    bool MeshAdapter::supportsContactElements() const
2402    {
2403      return true;
2404    }
2405    
2406    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2407    {
2408      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2409    }
2410    
2411    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2412    {
2413      Finley_ReferenceElementSet_dealloc(m_refSet);
2414    }
2415    
2416  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.2842  
changed lines
  Added in v.3490

  ViewVC Help
Powered by ViewVC 1.1.26