/[escript]/branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/diaplayground/finley/src/CPPAdapter/MeshAdapter.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3258 by gross, Thu Mar 18 01:45:55 2010 UTC revision 3259 by jfenwick, Mon Oct 11 01:48:14 2010 UTC
# Line 18  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21  #ifdef PASO_MPI  #ifdef ESYS_MPI
22  #include <mpi.h>  #include <mpi.h>
23  #include "paso/Paso_MPI.h"  #include "esysUtils/Esys_MPI.h"
24  #endif  #endif
25  extern "C" {  extern "C" {
26  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
# Line 85  int MeshAdapter::getMPIRank() const Line 85  int MeshAdapter::getMPIRank() const
85  }  }
86  void MeshAdapter::MPIBarrier() const  void MeshAdapter::MPIBarrier() const
87  {  {
88  #ifdef PASO_MPI  #ifdef ESYS_MPI
89     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);     MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
90  #endif  #endif
91     return;     return;
# Line 96  bool MeshAdapter::onMasterProcessor() co Line 96  bool MeshAdapter::onMasterProcessor() co
96  }  }
97    
98    
99  #ifdef PASO_MPI  #ifdef ESYS_MPI
100    MPI_Comm    MPI_Comm
101  #else  #else
102    unsigned int    unsigned int
103  #endif  #endif
104  MeshAdapter::getMPIComm() const  MeshAdapter::getMPIComm() const
105  {  {
106  #ifdef PASO_MPI  #ifdef ESYS_MPI
107      return m_finleyMesh->MPIInfo->comm;      return m_finleyMesh->MPIInfo->comm;
108  #else  #else
109      return 0;      return 0;
# Line 149  void MeshAdapter::dump(const string& fil Line 149  void MeshAdapter::dump(const string& fil
149     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
150     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
151     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
152  #ifdef PASO_MPI  #ifdef ESYS_MPI
153     MPI_Status status;     MPI_Status status;
154  #endif  #endif
155    
156  /* Incoming token indicates it's my turn to write */  /* Incoming token indicates it's my turn to write */
157  #ifdef PASO_MPI  #ifdef ESYS_MPI
158     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);
159  #endif  #endif
160    
161     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Esys_MPI_appendRankToFileName(fileName.c_str(),
162                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
163    
164     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
# Line 529  void MeshAdapter::dump(const string& fil Line 529  void MeshAdapter::dump(const string& fil
529     }     }
530    
531  /* Send token to next MPI process so he can take his turn */  /* Send token to next MPI process so he can take his turn */
532  #ifdef PASO_MPI  #ifdef ESYS_MPI
533     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);
534  #endif  #endif
535    
# Line 1044  void MeshAdapter::interpolateOnDomain(es Line 1044  void MeshAdapter::interpolateOnDomain(es
1044        case(Elements):        case(Elements):
1045        case(ReducedElements):        case(ReducedElements):
1046        if (getMPISize()>1) {        if (getMPISize()>1) {
1047           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1048           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1049           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1050        } else {        } else {
# Line 1054  void MeshAdapter::interpolateOnDomain(es Line 1054  void MeshAdapter::interpolateOnDomain(es
1054        case(FaceElements):        case(FaceElements):
1055        case(ReducedFaceElements):        case(ReducedFaceElements):
1056        if (getMPISize()>1) {        if (getMPISize()>1) {
1057           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1058           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1059           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1060        
# Line 1064  void MeshAdapter::interpolateOnDomain(es Line 1064  void MeshAdapter::interpolateOnDomain(es
1064        break;        break;
1065        case(Points):        case(Points):
1066        if (getMPISize()>1) {        if (getMPISize()>1) {
1067           escript::Data temp=escript::Data( in,  continuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1068           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1069        } else {        } else {
1070           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
# Line 1075  void MeshAdapter::interpolateOnDomain(es Line 1075  void MeshAdapter::interpolateOnDomain(es
1075        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1076        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
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           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1081        } else {        } else {
# Line 1113  void MeshAdapter::interpolateOnDomain(es Line 1113  void MeshAdapter::interpolateOnDomain(es
1113        case(Elements):        case(Elements):
1114        case(ReducedElements):        case(ReducedElements):
1115        if (getMPISize()>1) {        if (getMPISize()>1) {
1116           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1117           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1118           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in2,&_target);
1119        } else {        } else {
# Line 1123  void MeshAdapter::interpolateOnDomain(es Line 1123  void MeshAdapter::interpolateOnDomain(es
1123        case(FaceElements):        case(FaceElements):
1124        case(ReducedFaceElements):        case(ReducedFaceElements):
1125        if (getMPISize()>1) {        if (getMPISize()>1) {
1126           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1127           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1128           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&_in2,&_target);
1129        } else {        } else {
# Line 1132  void MeshAdapter::interpolateOnDomain(es Line 1132  void MeshAdapter::interpolateOnDomain(es
1132        break;        break;
1133        case(Points):        case(Points):
1134        if (getMPISize()>1) {        if (getMPISize()>1) {
1135           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1136           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1137           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in2,&_target);
1138        } else {        } else {
# Line 1144  void MeshAdapter::interpolateOnDomain(es Line 1144  void MeshAdapter::interpolateOnDomain(es
1144        case(ReducedContactElementsZero):        case(ReducedContactElementsZero):
1145        case(ReducedContactElementsOne):        case(ReducedContactElementsOne):
1146        if (getMPISize()>1) {        if (getMPISize()>1) {
1147           escript::Data temp=escript::Data( in,  reducedContinuousFunction(asAbstractContinuousDomain()) );           escript::Data temp=escript::Data( in,  reducedContinuousFunction(*this) );
1148           escriptDataC _in2 = temp.getDataC();           escriptDataC _in2 = temp.getDataC();
1149           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in2,&_target);
1150        } else {        } else {
# Line 1181  void MeshAdapter::setToX(escript::Data& Line 1181  void MeshAdapter::setToX(escript::Data&
1181        escriptDataC _arg=arg.getDataC();        escriptDataC _arg=arg.getDataC();
1182        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_arg);
1183     } else {     } else {
1184        escript::Data tmp_data=Vector(0.0,continuousFunction(asAbstractContinuousDomain()),true);        escript::Data tmp_data=Vector(0.0,continuousFunction(*this),true);
1185        escriptDataC _tmp_data=tmp_data.getDataC();        escriptDataC _tmp_data=tmp_data.getDataC();
1186        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);        Finley_Assemble_NodeCoordinates(mesh->Nodes,&_tmp_data);
1187        // this is then interpolated onto arg:        // this is then interpolated onto arg:
# Line 1275  void MeshAdapter::setToIntegrals(vector< Line 1275  void MeshAdapter::setToIntegrals(vector<
1275     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1276     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1277     case(Nodes):     case(Nodes):
1278     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1279     _temp=temp.getDataC();     _temp=temp.getDataC();
1280     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1281     break;     break;
1282     case(ReducedNodes):     case(ReducedNodes):
1283     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1284     _temp=temp.getDataC();     _temp=temp.getDataC();
1285     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1286     break;     break;
# Line 1312  void MeshAdapter::setToIntegrals(vector< Line 1312  void MeshAdapter::setToIntegrals(vector<
1312     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1313     break;     break;
1314     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1315     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1316     _temp=temp.getDataC();     _temp=temp.getDataC();
1317     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1318     break;     break;
1319     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1320     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(*this) );
1321     _temp=temp.getDataC();     _temp=temp.getDataC();
1322     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1323     break;     break;
# Line 1349  void MeshAdapter::setToGradient(escript: Line 1349  void MeshAdapter::setToGradient(escript:
1349     escript::Data temp;     escript::Data temp;
1350     if (getMPISize()>1) {     if (getMPISize()>1) {
1351        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
1352           temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  continuousFunction(*this) );
1353           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1354        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {        } else if( arg.getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom ) {
1355           temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );           temp=escript::Data( arg,  reducedContinuousFunction(*this) );
1356           nodeDataC = temp.getDataC();           nodeDataC = temp.getDataC();
1357        } else {        } else {
1358           nodeDataC = arg.getDataC();           nodeDataC = arg.getDataC();
# Line 1471  void MeshAdapter::setNewX(const escript: Line 1471  void MeshAdapter::setNewX(const escript:
1471     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1472     if (newDomain!=*this)     if (newDomain!=*this)
1473        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1474     if ( new_x.getFunctionSpace() == continuousFunction(asAbstractContinuousDomain()) ) {     if ( new_x.getFunctionSpace() == continuousFunction(*this) ) {
1475         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1476         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1477     } else {     } else {
1478         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(asAbstractContinuousDomain()) );         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1479         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1480         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1481     }     }
# Line 1559  void MeshAdapter::saveVTK(const string& Line 1559  void MeshAdapter::saveVTK(const string&
1559    
1560  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1561  {  {
1562  #ifdef PASO_MPI  #ifdef ESYS_MPI
1563      index_t myFirstNode=0, myLastNode=0, k=0;      index_t myFirstNode=0, myLastNode=0, k=0;
1564      index_t* globalNodeIndex=0;      index_t* globalNodeIndex=0;
1565      Finley_Mesh* mesh_p=m_finleyMesh.get();      Finley_Mesh* mesh_p=m_finleyMesh.get();
# Line 1586  bool MeshAdapter::ownSample(int fs_code, Line 1586  bool MeshAdapter::ownSample(int fs_code,
1586  //  //
1587  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1588  //  //
1589  SystemMatrixAdapter MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(
1590                                                   const int row_blocksize,                                                   const int row_blocksize,
1591                                                   const escript::FunctionSpace& row_functionspace,                                                   const escript::FunctionSpace& row_functionspace,
1592                                                   const int column_blocksize,                                                   const int column_blocksize,
# Line 1633  SystemMatrixAdapter MeshAdapter::newSyst Line 1633  SystemMatrixAdapter MeshAdapter::newSyst
1633     }     }
1634     checkPasoError();     checkPasoError();
1635     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1636     return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1637       return ASM_ptr(sma);
1638    //   return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1639  }  }
1640    
1641  //  //
1642  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1643  //  //
1644  TransportProblemAdapter MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
1645                                                           const bool useBackwardEuler,                                                           const bool useBackwardEuler,
1646                                                           const int blocksize,                                                           const int blocksize,
1647                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
# Line 1666  TransportProblemAdapter MeshAdapter::new Line 1668  TransportProblemAdapter MeshAdapter::new
1668     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);
1669     checkPasoError();     checkPasoError();
1670     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1671     return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1672       return ATP_ptr(tpa);
1673    //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);
1674  }  }
1675    
1676  //  //
# Line 2039  int MeshAdapter::getTransportTypeId(cons Line 2043  int MeshAdapter::getTransportTypeId(cons
2043    
2044  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
2045  {  {
2046     return continuousFunction(asAbstractContinuousDomain()).getX();     return continuousFunction(*this).getX();
2047  }  }
2048    
2049  escript::Data MeshAdapter::getNormal() const  escript::Data MeshAdapter::getNormal() const
2050  {  {
2051     return functionOnBoundary(asAbstractContinuousDomain()).getNormal();     return functionOnBoundary(*this).getNormal();
2052  }  }
2053    
2054  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
2055  {  {
2056     return escript::function(asAbstractContinuousDomain()).getSize();     return escript::function(*this).getSize();
2057  }  }
2058    
2059  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2397  int MeshAdapter::getApproximationOrder(c Line 2401  int MeshAdapter::getApproximationOrder(c
2401    return order;    return order;
2402  }  }
2403    
2404  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id, index_t order, index_t reducedOrder)  bool MeshAdapter::supportsContactElements() const
2405    {
2406      return true;
2407    }
2408    
2409    ReferenceElementSetWrapper::ReferenceElementSetWrapper(Finley_ElementTypeId id, index_t order, index_t reducedOrder)
2410  {  {
2411    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);    m_refSet = Finley_ReferenceElementSet_alloc(id, order, reducedOrder);
2412  }  }

Legend:
Removed from v.3258  
changed lines
  Added in v.3259

  ViewVC Help
Powered by ViewVC 1.1.26