/[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 3981 by jfenwick, Fri Sep 21 02:47:54 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2009 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * http://www.uq.edu.au
 * http://www.uq.edu.au/esscc  
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
8  * Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
9  * http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
10  *  *
11  *******************************************************/  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12    * Development since 2012 by School of Earth Sciences
13    *
14    *****************************************************************************/
15    
16    #include <pasowrap/PasoException.h>
17    #include <pasowrap/TransportProblemAdapter.h>
18  #include "MeshAdapter.h"  #include "MeshAdapter.h"
19  #include "escript/Data.h"  #include "escript/Data.h"
20  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
21  #ifdef USE_NETCDF  #ifdef USE_NETCDF
22  #include <netcdfcpp.h>  #include <netcdfcpp.h>
23  #endif  #endif
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #include "paso/Paso_MPI.h"  
 #endif  
24  extern "C" {  extern "C" {
25  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
26  }  }
27    
28    #include <boost/python/import.hpp>
29    
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    using namespace paso;
33    
34  namespace finley {  namespace finley {
35    
# Line 85  int MeshAdapter::getMPIRank() const Line 87  int MeshAdapter::getMPIRank() const
87  }  }
88  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
89  {  {
90  #ifdef PASO_MPI  #ifdef ESYS_MPI
91     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
92  #endif  #endif
93     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 98  bool MeshAdapter::onMasterProcessor() co
98  }  }
99    
100    
101  #ifdef PASO_MPI  #ifdef ESYS_MPI
102    MPI_Comm    MPI_Comm
103  #else  #else
104    unsigned int    unsigned int
105  #endif  #endif
106  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
107  {  {
108  #ifdef PASO_MPI  #ifdef ESYS_MPI
109      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
110  #else  #else
111      return 0;      return 0;
# Line 149  void MeshAdapter::dump(const string& fil Line 151  void MeshAdapter::dump(const string& fil
151     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
152     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
153     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
154  #ifdef PASO_MPI  #ifdef ESYS_MPI
155     MPI_Status status;     MPI_Status status;
156  #endif  #endif
157    
158  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
159  #ifdef PASO_MPI  #ifdef ESYS_MPI
160     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);     if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
161  #endif  #endif
162    
163     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
164                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
165    
166     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
# Line 220  void MeshAdapter::dump(const string& fil Line 222  void MeshAdapter::dump(const string& fil
222        throw DataException(msgPrefix+"add_att(Name)");        throw DataException(msgPrefix+"add_att(Name)");
223     if (!dataFile.add_att("numDim",numDim) )     if (!dataFile.add_att("numDim",numDim) )
224        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
225     if (!dataFile.add_att("order",mesh->order) )     if (!dataFile.add_att("order",mesh->integrationOrder) )
226        throw DataException(msgPrefix+"add_att(order)");        throw DataException(msgPrefix+"add_att(order)");
227     if (!dataFile.add_att("reduced_order",mesh->reduced_order) )     if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )
228        throw DataException(msgPrefix+"add_att(reduced_order)");        throw DataException(msgPrefix+"add_att(reduced_order)");
229     if (!dataFile.add_att("numNodes",numNodes) )     if (!dataFile.add_att("numNodes",numNodes) )
230        throw DataException(msgPrefix+"add_att(numNodes)");        throw DataException(msgPrefix+"add_att(numNodes)");
# Line 529  void MeshAdapter::dump(const string& fil Line 531  void MeshAdapter::dump(const string& fil
531     }     }
532    
533  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
534  #ifdef PASO_MPI  #ifdef ESYS_MPI
535     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);     if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
536  #endif  #endif
537    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 651  int MeshAdapter::getReducedSolutionCode(
651     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
652  }  }
653    
654  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
655  {  {
656     return Points;     return Points;
657  }  }
# Line 770  pair<int,int> MeshAdapter::getDataShape( Line 772  pair<int,int> MeshAdapter::getDataShape(
772  // adds linear PDE of second order into a given stiffness matrix and right hand side:  // adds linear PDE of second order into a given stiffness matrix and right hand side:
773  //  //
774  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(
775                                   SystemMatrixAdapter& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
776                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                   const escript::Data& A, const escript::Data& B, const escript::Data& C,const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
777                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
778                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
779                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
780  {  {
781       SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
782       if (smat==0)
783       {
784        throw FinleyAdapterException("finley only supports Paso system matrices.");
785       }
786     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
787     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
788     escriptDataC _B=B.getDataC();     escriptDataC _B=B.getDataC();
# Line 786  void MeshAdapter::addPDEToSystem( Line 794  void MeshAdapter::addPDEToSystem(
794     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
795     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
796     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
797       escriptDataC _d_dirac=d_dirac.getDataC();
798       escriptDataC _y_dirac=y_dirac.getDataC();
799    
800     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
801    
802     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
803       checkFinleyError();
804    
805       Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );
806     checkFinleyError();     checkFinleyError();
807    
808     Finley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, 0, 0, 0, &_d, 0, &_y );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );
809     checkFinleyError();     checkFinleyError();
810    
811     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_contact, 0, &_y_contact );      Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
812     checkFinleyError();     checkFinleyError();
813  }  }
814    
815  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
816                                          escript::Data& mat,                                          escript::Data& mat,
817                                          const escript::Data& D,                                          const escript::Data& D,
818                                          const escript::Data& d) const                                          const escript::Data& d,
819                                            const escript::Data& d_dirac,
820                                            const bool useHRZ) const
821  {  {
822     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
823     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
824     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
825       escriptDataC _d_dirac=d_dirac.getDataC();
826    
827     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
828    
829     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
830     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     checkFinleyError();
831      
832       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
833       checkFinleyError();
834    
835       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
836     checkFinleyError();     checkFinleyError();
837    
838  }  }
839    
840    
841  //  //
842  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
843  //  //
844  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact,  const escript::Data& y_dirac) const
845  {  {
846     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
847    
# Line 829  void MeshAdapter::addPDEToRHS( escript:: Line 850  void MeshAdapter::addPDEToRHS( escript::
850     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
851     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
852     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
853       escriptDataC _y_dirac=y_dirac.getDataC();
854    
855     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements, 0, &_rhs, 0, 0, 0, 0, &_X, &_Y );
856     checkFinleyError();     checkFinleyError();
# Line 838  void MeshAdapter::addPDEToRHS( escript:: Line 860  void MeshAdapter::addPDEToRHS( escript::
860    
861     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, 0, &_rhs , 0, 0, 0, 0, 0, &_y_contact );
862     checkFinleyError();     checkFinleyError();
863    
864       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
865       checkFinleyError();
866    
867  }  }
868  //  //
869  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
870  //  //
871  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
872                                             TransportProblemAdapter& tp, escript::Data& source, const escript::Data& M,                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,
873                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
874                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
875                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
876                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
877                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
878  {  {
879       TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
880       if (tpa==0)
881       {
882        throw FinleyAdapterException("finley only supports Paso transport problems.");
883       }
884    
885    
886     DataTypes::ShapeType shape;     DataTypes::ShapeType shape;
887     source.expand();     source.expand();
888     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
# Line 863  void MeshAdapter::addPDEToTransportProbl Line 897  void MeshAdapter::addPDEToTransportProbl
897     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
898     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
899     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
900       escriptDataC _d_dirac=d_dirac.getDataC();
901       escriptDataC _y_dirac=y_dirac.getDataC();
902    
903    
904     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
905     Paso_FCTransportProblem* _tp = tp.getPaso_FCTransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
906    
907     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &_source, 0, 0, 0, &_M, 0, 0 );
908     checkFinleyError();     checkFinleyError();
# Line 878  void MeshAdapter::addPDEToTransportProbl Line 915  void MeshAdapter::addPDEToTransportProbl
915    
916     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, _tp->transport_matrix, &_source , 0, 0, 0, &_d_contact, 0, &_y_contact );
917     checkFinleyError();     checkFinleyError();
918    
919       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
920       checkFinleyError();
921    
922  }  }
923    
924  //  //
# Line 1042  void MeshAdapter::interpolateOnDomain(es Line 1083  void MeshAdapter::interpolateOnDomain(es
1083        case(Elements):        case(Elements):
1084        case(ReducedElements):        case(ReducedElements):
1085        if (getMPISize()>1) {        if (getMPISize()>1) {
1086           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1087           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1088           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1089        } else {        } else {
# Line 1052  void MeshAdapter::interpolateOnDomain(es Line 1093  void MeshAdapter::interpolateOnDomain(es
1093        case(FaceElements):        case(FaceElements):
1094        case(ReducedFaceElements):        case(ReducedFaceElements):
1095        if (getMPISize()>1) {        if (getMPISize()>1) {
1096           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1097           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1098           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1099        
# Line 1062  void MeshAdapter::interpolateOnDomain(es Line 1103  void MeshAdapter::interpolateOnDomain(es
1103        break;        break;
1104        case(Points):        case(Points):
1105        if (getMPISize()>1) {        if (getMPISize()>1) {
1106           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1107           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1108        } else {        } else {
1109           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1110        }        }
# Line 1073  void MeshAdapter::interpolateOnDomain(es Line 1114  void MeshAdapter::interpolateOnDomain(es
1114        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1115        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1116        if (getMPISize()>1) {        if (getMPISize()>1) {
1117           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1118           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1119           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1120        } else {        } else {
# Line 1111  void MeshAdapter::interpolateOnDomain(es Line 1152  void MeshAdapter::interpolateOnDomain(es
1152        case(Elements):        case(Elements):
1153        case(ReducedElements):        case(ReducedElements):
1154        if (getMPISize()>1) {        if (getMPISize()>1) {
1155           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1156           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1157           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1158        } else {        } else {
# Line 1121  void MeshAdapter::interpolateOnDomain(es Line 1162  void MeshAdapter::interpolateOnDomain(es
1162        case(FaceElements):        case(FaceElements):
1163        case(ReducedFaceElements):        case(ReducedFaceElements):
1164        if (getMPISize()>1) {        if (getMPISize()>1) {
1165           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1166           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1167           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1168        } else {        } else {
# Line 1130  void MeshAdapter::interpolateOnDomain(es Line 1171  void MeshAdapter::interpolateOnDomain(es
1171        break;        break;
1172        case(Points):        case(Points):
1173        if (getMPISize()>1) {        if (getMPISize()>1) {
1174           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1175           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1176           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1177        } else {        } else {
# Line 1142  void MeshAdapter::interpolateOnDomain(es Line 1183  void MeshAdapter::interpolateOnDomain(es
1183        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1184        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1185        if (getMPISize()>1) {        if (getMPISize()>1) {
1186           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1187           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1188           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1189        } else {        } else {
# Line 1179  void MeshAdapter::setToX(escript::Data& Line 1220  void MeshAdapter::setToX(escript::Data&
1220        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1221        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1222     } else {     } else {
1223        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1224        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1225        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1226        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1273  void MeshAdapter::setToIntegrals(vector< Line 1314  void MeshAdapter::setToIntegrals(vector<
1314     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1315     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1316     case(Nodes):     case(Nodes):
1317     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1318     _temp=temp.getDataC();     _temp=temp.getDataC();
1319     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1320     break;     break;
1321     case(ReducedNodes):     case(ReducedNodes):
1322     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1323     _temp=temp.getDataC();     _temp=temp.getDataC();
1324     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1325     break;     break;
# Line 1310  void MeshAdapter::setToIntegrals(vector< Line 1351  void MeshAdapter::setToIntegrals(vector<
1351     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1352     break;     break;
1353     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1354     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1355     _temp=temp.getDataC();     _temp=temp.getDataC();
1356     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1357     break;     break;
1358     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1359     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1360     _temp=temp.getDataC();     _temp=temp.getDataC();
1361     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1362     break;     break;
# Line 1347  void MeshAdapter::setToGradient(escript: Line 1388  void MeshAdapter::setToGradient(escript:
1388     escript::Data temp;     escript::Data temp;
1389     if (getMPISize()>1) {     if (getMPISize()>1) {
1390        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1391           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1392           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1393        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1394           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1395           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1396        } else {        } else {
1397           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1469  void MeshAdapter::setNewX(const escript: Line 1510  void MeshAdapter::setNewX(const escript:
1510     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1511     if (newDomain!=*this)     if (newDomain!=*this)
1512        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1513     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1514         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1515         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1516     } else {     } else {
1517         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1518    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1519         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1520         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
1521     }     }
1522     checkFinleyError();     checkFinleyError();
1523  }  }
# Line 1536  void MeshAdapter::saveDX(const string& f Line 1578  void MeshAdapter::saveDX(const string& f
1578  //  //
1579  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
1580  {  {
1581     int num_data;      boost::python::object pySaveVTK = boost::python::import("esys.weipa").attr("_saveVTK");
1582     char **names;      pySaveVTK(filename, const_cast<MeshAdapter*>(this)->getPtr(),
1583     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);  
1584  }  }
1585    
1586  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1587  {  {
1588  #ifdef PASO_MPI  #ifdef ESYS_MPI
1589      index_t myFirstNode=0, myLastNode=0, k=0;      index_t myFirstNode=0, myLastNode=0, k=0;
1590      index_t* globalNodeIndex=0;      index_t* globalNodeIndex=0;
1591      Finley_Mesh* mesh_p=m_finleyMesh.get();      Finley_Mesh* mesh_p=m_finleyMesh.get();
# Line 1584  bool MeshAdapter::ownSample(int fs_code, Line 1612  bool MeshAdapter::ownSample(int fs_code,
1612  //  //
1613  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1614  //  //
1615  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1616                                                   const int row_blocksize,                                                   const int row_blocksize,
1617                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1618                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1631  SystemMatrixAdapter MeshAdapter::newSyst Line 1659  SystemMatrixAdapter MeshAdapter::newSyst
1659     }     }
1660     checkPasoError();     checkPasoError();
1661     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1662     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1663       return ASM_ptr(sma);
1664    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1665  }  }
1666    
1667  //  //
1668  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1669  //  //
1670  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const double theta,  
1671                                                           const int blocksize,                                                           const int blocksize,
1672                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1673                                                           const int type) const                                                           const int type) const
# Line 1660  TransportProblemAdapter MeshAdapter::new Line 1689  TransportProblemAdapter MeshAdapter::new
1689    
1690     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1691     checkFinleyError();     checkFinleyError();
1692     Paso_FCTransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1693     transportProblem=Paso_FCTransportProblem_alloc(theta,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1694     checkPasoError();     checkPasoError();
1695     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1696     return TransportProblemAdapter(transportProblem,theta,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1697       return ATP_ptr(tpa);
1698    //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1699  }  }
1700    
1701  //  //
# Line 1727  MeshAdapter::commonFunctionSpace(const v Line 1758  MeshAdapter::commonFunctionSpace(const v
1758     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.
1759     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.
1760     */     */
1761      if (fs.size()==0)      if (fs.empty())
1762      {      {
1763          return false;          return false;
1764      }      }
# Line 2023  bool MeshAdapter::operator!=(const Abstr Line 2054  bool MeshAdapter::operator!=(const Abstr
2054  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
2055  {  {
2056     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2057     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2058     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2059  }  }
2060    
2061  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
2062  {  {
2063     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2064     int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2065     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2066  }  }
2067    
2068  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2069  {  {
2070     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2071  }  }
2072    
2073  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2074  {  {
2075     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2076  }  }
2077    
2078  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2079  {  {
2080     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2081  }  }
2082    
2083  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2213  void MeshAdapter::setTagMap(const string Line 2243  void MeshAdapter::setTagMap(const string
2243  {  {
2244     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2245     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2246     checkPasoError();     checkFinleyError();
2247     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2248  }  }
2249    
# Line 2222  int MeshAdapter::getTag(const string& na Line 2252  int MeshAdapter::getTag(const string& na
2252     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2253     int tag=0;     int tag=0;
2254     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2255     checkPasoError();     checkFinleyError();
2256     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2257     return tag;     return tag;
2258  }  }
# Line 2360  AbstractDomain::StatusType MeshAdapter:: Line 2390  AbstractDomain::StatusType MeshAdapter::
2390    return Finley_Mesh_getStatus(mesh);    return Finley_Mesh_getStatus(mesh);
2391  }  }
2392    
2393    int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2394    {
2395      
2396      Finley_Mesh* mesh=m_finleyMesh.get();
2397      int order =-1;
2398      switch(functionSpaceCode) {
2399       case(Nodes):
2400       case(DegreesOfFreedom):
2401              order=mesh->approximationOrder;
2402              break;
2403       case(ReducedNodes):
2404       case(ReducedDegreesOfFreedom):
2405              order=mesh->reducedApproximationOrder;
2406              break;
2407       case(Elements):
2408       case(FaceElements):
2409       case(Points):
2410       case(ContactElementsZero):
2411       case(ContactElementsOne):
2412              order=mesh->integrationOrder;
2413              break;
2414       case(ReducedElements):
2415       case(ReducedFaceElements):
2416       case(ReducedContactElementsZero):
2417       case(ReducedContactElementsOne):
2418              order=mesh->reducedIntegrationOrder;
2419              break;
2420       default:
2421          stringstream temp;
2422          temp << "Error - Finley does not know anything about function space type " << functionSpaceCode;
2423          throw FinleyAdapterException(temp.str());
2424      }
2425      return order;
2426    }
2427    
2428    bool MeshAdapter::supportsContactElements() const
2429    {
2430      return true;
2431    }
2432    
2433    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2434    {
2435      m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2436    }
2437    
2438    ReferenceElementSetWrapper::~ReferenceElementSetWrapper()
2439    {
2440      Finley_ReferenceElementSet_dealloc(m_refSet);
2441    }
2442    
2443    // points will be flattened
2444    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2445    {
2446          const int dim = getDim();
2447          int numPoints=points.size()/dim;
2448          int numTags=tags.size();
2449          Finley_Mesh* mesh=m_finleyMesh.get();
2450          
2451          if ( points.size() % dim != 0 )
2452          {
2453        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2454          }
2455          
2456          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2457         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2458          
2459          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2460          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2461          
2462          for (int i=0;i<numPoints;++i) {
2463           points_ptr[ i * dim     ] = points[i * dim ];
2464           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2465           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2466               tags_ptr[i]=tags[i];
2467          }
2468          
2469          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2470          checkFinleyError();
2471          
2472          TMPMEMFREE(points_ptr);
2473          TMPMEMFREE(tags_ptr);
2474    }
2475    
2476    
2477    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2478    // {
2479    //       const int dim = getDim();
2480    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2481    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2482    //       Finley_Mesh* mesh=m_finleyMesh.get();
2483    //      
2484    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2485    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2486    //      
2487    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2488    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2489    //      
2490    //       for (int i=0;i<numPoints;++i) {
2491    //     int tag_id=-1;
2492    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2493    //     if  ( numComps !=   dim ) {
2494    //                stringstream temp;            
2495    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2496    //                throw FinleyAdapterException(temp.str());
2497    //     }
2498    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2499    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2500    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2501    //    
2502    //     if ( numTags > 0) {
2503    //            boost::python::extract<string> ex_str(tags[i]);
2504    //        if  ( ex_str.check() ) {
2505    //            tag_id=getTag( ex_str());
2506    //        } else {
2507    //             boost::python::extract<int> ex_int(tags[i]);
2508    //             if ( ex_int.check() ) {
2509    //                 tag_id=ex_int();
2510    //             } else {
2511    //              stringstream temp;          
2512    //                  temp << "Error - unable to extract tag for point " << i;
2513    //              throw FinleyAdapterException(temp.str());
2514    //            }
2515    //        }
2516    //     }      
2517    //            tags_ptr[i]=tag_id;
2518    //       }
2519    //      
2520    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2521    //       checkPasoError();
2522    //      
2523    //       TMPMEMFREE(points_ptr);
2524    //       TMPMEMFREE(tags_ptr);
2525    // }
2526    
2527    /*
2528    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2529    {  
2530        boost::python::list points =  boost::python::list();
2531        boost::python::list tags =  boost::python::list();
2532        points.append(point);
2533        tags.append(tag);
2534        addDiracPoints(points, tags);
2535    }
2536    */
2537    
2538    /*
2539    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2540    {
2541            boost::python::list points =   boost::python::list();
2542            boost::python::list tags =   boost::python::list();
2543            points.append(point);
2544            tags.append(tag);
2545            addDiracPoints(points, tags);
2546    }
2547    */
2548  }  // end of namespace  }  // end of namespace

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

  ViewVC Help
Powered by ViewVC 1.1.26