/[escript]/trunk-mpi-branch/finley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /trunk-mpi-branch/finley/src/CPPAdapter/MeshAdapter.cpp

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

revision 1222 by ksteube, Fri Jun 15 03:45:48 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 13  Line 13 
13   ******************************************************************************   ******************************************************************************
14  */  */
15    
 #ifdef PASO_MPI  
 #include <mpi.h>  
 #endif  
16  #include "MeshAdapter.h"  #include "MeshAdapter.h"
   
17  #include "escript/Data.h"  #include "escript/Data.h"
18  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
19  extern "C" {  extern "C" {
20  #include "escript/blocktimer.h"  #include "escript/blocktimer.h"
21  }  }
22    #include <vector>
23    
24  using namespace std;  using namespace std;
25  using namespace escript;  using namespace escript;
# Line 50  MeshAdapter::MeshAdapter(Finley_Mesh* fi Line 47  MeshAdapter::MeshAdapter(Finley_Mesh* fi
47  {  {
48    setFunctionSpaceTypeNames();    setFunctionSpaceTypeNames();
49    //    //
50    // need to use a null_deleter as Finley_Mesh_dealloc deletes the pointer    // need to use a null_deleter as Finley_Mesh_free deletes the pointer
51    // for us.    // for us.
52    m_finleyMesh.reset(finleyMesh,null_deleter());    m_finleyMesh.reset(finleyMesh,null_deleter());
53  }  }
# Line 69  MeshAdapter::~MeshAdapter() Line 66  MeshAdapter::~MeshAdapter()
66    // I hope the case for the pointer being zero has been taken care of.    // I hope the case for the pointer being zero has been taken care of.
67    //  cout << "In MeshAdapter destructor." << endl;    //  cout << "In MeshAdapter destructor." << endl;
68    if (m_finleyMesh.unique()) {    if (m_finleyMesh.unique()) {
69      //   cout << "Calling dealloc." << endl;      Finley_Mesh_free(m_finleyMesh.get());
     Finley_Mesh_dealloc(m_finleyMesh.get());  
     //   cout << "Finished dealloc." << endl;  
70    }    }
71  }  }
72    
73    int MeshAdapter::getMPISize() const
74    {
75       return m_finleyMesh.get()->MPIInfo->size;
76    }
77    int MeshAdapter::getMPIRank() const
78    {
79       return m_finleyMesh.get()->MPIInfo->rank;
80    }
81    
82    
83  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
84     return m_finleyMesh.get();     return m_finleyMesh.get();
85  }  }
# Line 223  pair<int,int> MeshAdapter::getDataShape( Line 228  pair<int,int> MeshAdapter::getDataShape(
228     switch (functionSpaceCode) {     switch (functionSpaceCode) {
229        case(Nodes):        case(Nodes):
230             numDataPointsPerSample=1;             numDataPointsPerSample=1;
231             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;             numSamples=Finley_NodeFile_getNumNodes(mesh->Nodes);
232             break;             break;
233        case(ReducedNodes):        case(ReducedNodes):
            /* TODO: add ReducedNodes */  
234             numDataPointsPerSample=1;             numDataPointsPerSample=1;
235             if (mesh->Nodes!=NULL) numSamples=mesh->Nodes->numNodes;             numSamples=Finley_NodeFile_getNumReducedNodes(mesh->Nodes);
            throw FinleyAdapterException("Error - ReducedNodes is not supported yet.");  
236             break;             break;
237        case(Elements):        case(Elements):
238             if (mesh->Elements!=NULL) {             if (mesh->Elements!=NULL) {
# Line 288  pair<int,int> MeshAdapter::getDataShape( Line 291  pair<int,int> MeshAdapter::getDataShape(
291        case(DegreesOfFreedom):        case(DegreesOfFreedom):
292             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
293               numDataPointsPerSample=1;               numDataPointsPerSample=1;
294  #ifndef PASO_MPI               numSamples=Finley_NodeFile_getNumDegreesOfFreedom(mesh->Nodes);
              numSamples=mesh->Nodes->numDegreesOfFreedom;  
 #else  
              numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;  
 #endif  
295             }             }
296             break;             break;
297        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
298             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
299               numDataPointsPerSample=1;               numDataPointsPerSample=1;
300  #ifndef PASO_MPI               numSamples=Finley_NodeFile_getNumReducedDegreesOfFreedom(mesh->Nodes);
              numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;  
 #else  
              numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;  
 #endif  
301             }             }
302             break;             break;
303        default:        default:
# Line 347  void MeshAdapter::addPDEToSystem( Line 342  void MeshAdapter::addPDEToSystem(
342     checkFinleyError();     checkFinleyError();
343  }  }
344    
345    void  MeshAdapter::addPDEToLumpedSystem(
346                         escript::Data& mat,
347                         const escript::Data& D,
348                         const escript::Data& d) const
349    {
350       escriptDataC _mat=mat.getDataC();
351       escriptDataC _D=D.getDataC();
352       escriptDataC _d=d.getDataC();
353    
354       Finley_Mesh* mesh=m_finleyMesh.get();
355    
356       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);
357       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);
358    
359       checkFinleyError();
360    }
361    
362    
363  //  //
364  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
365  //  //
# Line 521  void MeshAdapter::interpolateOnDomain(es Line 534  void MeshAdapter::interpolateOnDomain(es
534             case(ReducedNodes):             case(ReducedNodes):
535                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);                Finley_Assemble_CopyNodalData(mesh->Nodes,&_target,&_in);
536                break;                break;
 #ifndef PASO_MPI  
537             case(Elements):             case(Elements):
538             case(ReducedElements):             case(ReducedElements):
539                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&_in,&_target);
# Line 539  void MeshAdapter::interpolateOnDomain(es Line 551  void MeshAdapter::interpolateOnDomain(es
551             case(ReducedContactElementsOne):             case(ReducedContactElementsOne):
552                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&_in,&_target);
553               break;               break;
 #else  
            /* need to copy Degrees of freedom data to nodal data so that the external values are available */  
            case(Elements):  
            case(ReducedElements):  
           { /* limit scope for variable nodeTemp */  
               escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );  
               Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&_target);  
           }  
               break;  
            case(FaceElements):  
            case(ReducedFaceElements):  
           {  
               escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );  
               Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&_target);  
           }  
               break;  
            case(Points):  
           {  
               escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );  
               Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&_target);  
           }  
               break;  
            case(ContactElementsZero):  
            case(ContactElementsOne):  
            case(ReducedContactElementsZero):  
            case(ReducedContactElementsOne):  
           {  
               escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );  
               Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&_target);  
           }  
              break;  
 #endif  
554             default:             default:
555               stringstream temp;               stringstream temp;
556               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
# Line 726  void MeshAdapter::setToIntegrals(std::ve Line 706  void MeshAdapter::setToIntegrals(std::ve
706    
707    double blocktimer_start = blocktimer_time();    double blocktimer_start = blocktimer_time();
708    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
709      escriptDataC _temp;
710      escript::Data temp;
711    escriptDataC _arg=arg.getDataC();    escriptDataC _arg=arg.getDataC();
712    switch(arg.getFunctionSpace().getTypeCode()) {    switch(arg.getFunctionSpace().getTypeCode()) {
713       case(Nodes):       case(Nodes):
714          /* TODO */          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
715          throw FinleyAdapterException("Error - Integral of data on nodes is not supported.");          _temp=temp.getDataC();
716            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
717          break;          break;
718       case(ReducedNodes):       case(ReducedNodes):
719          /* TODO */          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
720          throw FinleyAdapterException("Error - Integral of data on reduced nodes is not supported.");          _temp=temp.getDataC();
721            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
722          break;          break;
723       case(Elements):       case(Elements):
724          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_arg,&integrals[0]);
# Line 764  void MeshAdapter::setToIntegrals(std::ve Line 748  void MeshAdapter::setToIntegrals(std::ve
748          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);          Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
749          break;          break;
750       case(DegreesOfFreedom):       case(DegreesOfFreedom):
751          throw FinleyAdapterException("Error - Integral of data on degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
752            _temp=temp.getDataC();
753            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
754          break;          break;
755       case(ReducedDegreesOfFreedom):       case(ReducedDegreesOfFreedom):
756          throw FinleyAdapterException("Error - Integral of data on reduced degrees of freedom is not supported.");          temp=escript::Data( arg, function(asAbstractContinuousDomain()) );
757            _temp=temp.getDataC();
758            Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
759          break;          break;
760       default:       default:
761          stringstream temp;          stringstream temp;
# Line 794  void MeshAdapter::setToGradient(escript: Line 782  void MeshAdapter::setToGradient(escript:
782    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
783    escriptDataC _grad=grad.getDataC();    escriptDataC _grad=grad.getDataC();
784    escriptDataC nodeDataC;    escriptDataC nodeDataC;
785  #ifdef PASO_MPI    escript::Data temp;
786    escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );    if (getMPISize()>1) {
787    if( arg.getFunctionSpace().getTypeCode() != Nodes )        if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
788    {          temp=escript::Data( arg,  continuousFunction(asAbstractContinuousDomain()) );
789      Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));          nodeDataC = temp.getDataC();
790      nodeDataC = nodeTemp.getDataC();        } else if( arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom ) {
791            temp=escript::Data( arg,  reducedContinuousFunction(asAbstractContinuousDomain()) );
792            nodeDataC = temp.getDataC();
793          } else {
794            nodeDataC = arg.getDataC();
795          }
796      } else {
797         nodeDataC = arg.getDataC();
798    }    }
   else  
     nodeDataC = arg.getDataC();  
 #else  
   nodeDataC = arg.getDataC();  
 #endif  
799    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
800         case(Nodes):         case(Nodes):
801            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
# Line 935  void MeshAdapter::saveDX(const std::stri Line 925  void MeshAdapter::saveDX(const std::stri
925    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;    char* *c_names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;
926    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;    escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
927    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;    escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
928      std::vector<escript::Data> data_in;
929    
930      boost::python::list keys=arg.keys();      boost::python::list keys=arg.keys();
931      for (int i=0;i<num_data;++i) {      for (int i=0;i<num_data;++i) {
932           escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);           data_in[i]=boost::python::extract<escript::Data&>(arg[keys[i]]);
933           if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)           if (dynamic_cast<const MeshAdapter&>(data_in[i].getFunctionSpace().getDomain()) !=*this)
934               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");               throw FinleyAdapterException("Error  in saveDX: Data must be defined on same Domain");
935           data[i]=d.getDataC();           if ( (data_in[i].getFunctionSpace().getTypeCode() == DegreesOfFreedom) ||
936                  (data_in[i].getFunctionSpace().getTypeCode() == ReducedDegreesOfFreedom) ||
937                  (data_in[i].getFunctionSpace().getTypeCode() == Nodes)  ) {
938                data_in[i]=escript::Data(data_in[i],reducedContinuousFunction(asAbstractContinuousDomain()));
939                data[i]=data_in[i].getDataC();
940             } else {
941                data[i]=data_in[i].getDataC();
942             }
943           ptr_data[i]=&(data[i]);           ptr_data[i]=&(data[i]);
944           std::string n=boost::python::extract<std::string>(keys[i]);           std::string n=boost::python::extract<std::string>(keys[i]);
945           c_names[i]=&(names[i][0]);           c_names[i]=&(names[i][0]);
# Line 1000  void MeshAdapter::saveVTK(const std::str Line 998  void MeshAdapter::saveVTK(const std::str
998              strcpy(c_names[i],n.c_str());              strcpy(c_names[i],n.c_str());
999           }           }
1000      }      }
 #ifndef PASO_MPI      
1001      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);      Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);
 #else  
     Finley_Mesh_saveVTK_MPIO(filename.c_str(),m_finleyMesh.get(),num_data,c_names,ptr_data);  
 #endif  
1002    
1003  checkFinleyError();    checkFinleyError();
1004    /* win32 refactor */    /* win32 refactor */
1005    TMPMEMFREE(c_names);    TMPMEMFREE(c_names);
1006    TMPMEMFREE(data);    TMPMEMFREE(data);
# Line 1068  SystemMatrixAdapter MeshAdapter::newSyst Line 1062  SystemMatrixAdapter MeshAdapter::newSyst
1062        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);        fsystemMatrix=Paso_SystemMatrix_alloc(type,fsystemMatrixPattern,row_blocksize,column_blocksize);
1063      }      }
1064      checkPasoError();      checkPasoError();
1065      Paso_SystemMatrixPattern_dealloc(fsystemMatrixPattern);      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1066      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);      return SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace,column_blocksize,column_functionspace);
1067  }  }
1068    
# Line 1309  int* MeshAdapter::borrowSampleReferenceI Line 1303  int* MeshAdapter::borrowSampleReferenceI
1303    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
1304    switch (functionSpaceType) {    switch (functionSpaceType) {
1305      case(Nodes):      case(Nodes):
1306        if (mesh->Nodes!=NULL) {        out=mesh->Nodes->Id;
1307          out=mesh->Nodes->Id;        break;
         break;  
       }  
1308      case(ReducedNodes):      case(ReducedNodes):
1309        throw FinleyAdapterException("Error -  ReducedNodes not supported yet.");        out=mesh->Nodes->reducedNodesId;
1310        break;        break;
1311      case(Elements):      case(Elements):
1312        out=mesh->Elements->Id;        out=mesh->Elements->Id;
# Line 1344  int* MeshAdapter::borrowSampleReferenceI Line 1336  int* MeshAdapter::borrowSampleReferenceI
1336        out=mesh->ContactElements->Id;        out=mesh->ContactElements->Id;
1337        break;        break;
1338      case(DegreesOfFreedom):      case(DegreesOfFreedom):
1339        out=mesh->Nodes->degreeOfFreedomId;        out=mesh->Nodes->degreesOfFreedomId;
1340        break;        break;
1341      case(ReducedDegreesOfFreedom):      case(ReducedDegreesOfFreedom):
1342        out=mesh->Nodes->reducedDegreeOfFreedomId;        out=mesh->Nodes->reducedDegreesOfFreedomId;
1343        break;        break;
1344      default:      default:
1345        stringstream temp;        stringstream temp;

Legend:
Removed from v.1222  
changed lines
  Added in v.1223

  ViewVC Help
Powered by ViewVC 1.1.26