/[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 3522 by gross, Tue May 24 00:57:58 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 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();     checkFinleyError();
800    
801     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 );
802     checkFinleyError();     checkFinleyError();
803    
804     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 );
805       checkFinleyError();
806    
807        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 1727  MeshAdapter::commonFunctionSpace(const v Line 1754  MeshAdapter::commonFunctionSpace(const v
1754     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.
1755     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.
1756     */     */
1757      if (fs.size()==0)      if (fs.empty())
1758      {      {
1759          return false;          return false;
1760      }      }
# 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 2360  AbstractDomain::StatusType MeshAdapter:: Line 2387  AbstractDomain::StatusType MeshAdapter::
2387    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2388  }  }
2389    
2390    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2391    {
2392      
2393      Finley_Mesh* mesh=m_finleyMesh.get();
2394      int order =-1;
2395      switch(functionSpaceCode) {
2396       case(Nodes):
2397       case(DegreesOfFreedom):
2398              order=mesh->approximationOrder;
2399              break;
2400       case(ReducedNodes):
2401       case(ReducedDegreesOfFreedom):
2402              order=mesh->reducedApproximationOrder;
2403              break;
2404       case(Elements):
2405       case(FaceElements):
2406       case(Points):
2407       case(ContactElementsZero):
2408       case(ContactElementsOne):
2409              order=mesh->integrationOrder;
2410              break;
2411       case(ReducedElements):
2412       case(ReducedFaceElements):
2413       case(ReducedContactElementsZero):
2414       case(ReducedContactElementsOne):
2415              order=mesh->reducedIntegrationOrder;
2416              break;
2417       default:
2418          stringstream temp;
2419          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2420          throw FinleyAdapterException(temp.str());
2421      }
2422      return order;
2423    }
2424    
2425    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);
2433    }
2434    
2435    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2436    {
2437      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) and ( 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.2842  
changed lines
  Added in v.3522

  ViewVC Help
Powered by ViewVC 1.1.26