/[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 3269 by jfenwick, Wed Oct 13 03:21:50 2010 UTC revision 4346 by jfenwick, Tue Apr 2 04:46:45 2013 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2010 by University of Queensland  * Copyright (c) 2003-2013 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 ESYS_MPI  
 #include <mpi.h>  
 #include "esysUtils/Esys_MPI.h"  
 #endif  
 extern "C" {  
24  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
25  }  
26    #include <boost/python/import.hpp>
27    
28  using namespace std;  using namespace std;
29  using namespace escript;  using namespace escript;
30    using namespace paso;
31    
32  namespace finley {  namespace finley {
33    
# Line 649  int MeshAdapter::getReducedSolutionCode( Line 649  int MeshAdapter::getReducedSolutionCode(
649     return ReducedDegreesOfFreedom;     return ReducedDegreesOfFreedom;
650  }  }
651    
652  int MeshAdapter::getDiracDeltaFunctionCode() const  int MeshAdapter::getDiracDeltaFunctionsCode() const
653  {  {
654     return Points;     return Points;
655  }  }
# Line 773  void MeshAdapter::addPDEToSystem( Line 773  void MeshAdapter::addPDEToSystem(
773                                   AbstractSystemMatrix& mat, escript::Data& rhs,                                   AbstractSystemMatrix& mat, escript::Data& rhs,
774                                   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,
775                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
776                                   const escript::Data& d_contact,const escript::Data& y_contact) const                                   const escript::Data& d_contact,const escript::Data& y_contact,
777                                     const escript::Data& d_dirac,const escript::Data& y_dirac) const
778  {  {
779     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);
780     if (smat==0)     if (smat==0)
781     {     {
782      throw FinleyAdapterException("finley only supports its own system matrices.");      throw FinleyAdapterException("finley only supports Paso system matrices.");
783     }     }
784     escriptDataC _rhs=rhs.getDataC();     escriptDataC _rhs=rhs.getDataC();
785     escriptDataC _A  =A.getDataC();     escriptDataC _A  =A.getDataC();
# Line 791  void MeshAdapter::addPDEToSystem( Line 792  void MeshAdapter::addPDEToSystem(
792     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
793     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
794     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
795       escriptDataC _d_dirac=d_dirac.getDataC();
796       escriptDataC _y_dirac=y_dirac.getDataC();
797    
798     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
799    
# Line 802  void MeshAdapter::addPDEToSystem( Line 805  void MeshAdapter::addPDEToSystem(
805    
806     Finley_Assemble_PDE(mesh->Nodes,mesh->ContactElements, smat->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 );
807     checkFinleyError();     checkFinleyError();
808    
809        Finley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &_rhs , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
810       checkFinleyError();
811  }  }
812    
813  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(
814                                          escript::Data& mat,                                          escript::Data& mat,
815                                          const escript::Data& D,                                          const escript::Data& D,
816                                          const escript::Data& d) const                                          const escript::Data& d,
817                                            const escript::Data& d_dirac,
818                                            const bool useHRZ) const
819  {  {
820     escriptDataC _mat=mat.getDataC();     escriptDataC _mat=mat.getDataC();
821     escriptDataC _D=D.getDataC();     escriptDataC _D=D.getDataC();
822     escriptDataC _d=d.getDataC();     escriptDataC _d=d.getDataC();
823       escriptDataC _d_dirac=d_dirac.getDataC();
824    
825     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
826    
827     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&_mat, &_D, useHRZ);
828     checkFinleyError();     checkFinleyError();
829        
830     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d);     Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&_mat, &_d, useHRZ);
831       checkFinleyError();
832    
833       Finley_Assemble_LumpedSystem(mesh->Nodes,mesh->Points,&_mat, &_d_dirac, useHRZ);
834     checkFinleyError();     checkFinleyError();
835    
836  }  }
# Line 827  void  MeshAdapter::addPDEToLumpedSystem( Line 839  void  MeshAdapter::addPDEToLumpedSystem(
839  //  //
840  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
841  //  //
842  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
843  {  {
844     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
845    
# Line 836  void MeshAdapter::addPDEToRHS( escript:: Line 848  void MeshAdapter::addPDEToRHS( escript::
848     escriptDataC _Y=Y.getDataC();     escriptDataC _Y=Y.getDataC();
849     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
850     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
851       escriptDataC _y_dirac=y_dirac.getDataC();
852    
853     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 );
854     checkFinleyError();     checkFinleyError();
# Line 845  void MeshAdapter::addPDEToRHS( escript:: Line 858  void MeshAdapter::addPDEToRHS( escript::
858    
859     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 );
860     checkFinleyError();     checkFinleyError();
861    
862       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, 0, &_rhs , 0, 0, 0, 0, 0, &_y_dirac );
863       checkFinleyError();
864    
865  }  }
866  //  //
867  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
# Line 854  void MeshAdapter::addPDEToTransportProbl Line 871  void MeshAdapter::addPDEToTransportProbl
871                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,
872                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,
873                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
874                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact,
875                                               const escript::Data& d_dirac, const escript::Data& y_dirac) const
876  {  {
877     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);
878     if (tpa==0)     if (tpa==0)
879     {     {
880      throw FinleyAdapterException("finley only supports its own transport problems.");      throw FinleyAdapterException("finley only supports Paso transport problems.");
881     }     }
882    
883    
# Line 877  void MeshAdapter::addPDEToTransportProbl Line 895  void MeshAdapter::addPDEToTransportProbl
895     escriptDataC _y=y.getDataC();     escriptDataC _y=y.getDataC();
896     escriptDataC _d_contact=d_contact.getDataC();     escriptDataC _d_contact=d_contact.getDataC();
897     escriptDataC _y_contact=y_contact.getDataC();     escriptDataC _y_contact=y_contact.getDataC();
898       escriptDataC _d_dirac=d_dirac.getDataC();
899       escriptDataC _y_dirac=y_dirac.getDataC();
900    
901    
902     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
903     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();     Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
# Line 892  void MeshAdapter::addPDEToTransportProbl Line 913  void MeshAdapter::addPDEToTransportProbl
913    
914     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 );
915     checkFinleyError();     checkFinleyError();
916    
917       Finley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &_source , 0, 0, 0, &_d_dirac, 0, &_y_dirac );
918       checkFinleyError();
919    
920  }  }
921    
922  //  //
# Line 1076  void MeshAdapter::interpolateOnDomain(es Line 1101  void MeshAdapter::interpolateOnDomain(es
1101        break;        break;
1102        case(Points):        case(Points):
1103        if (getMPISize()>1) {        if (getMPISize()>1) {
1104           escript::Data temp=escript::Data( in,  continuousFunction(*this) );           //escript::Data temp=escript::Data( in,  continuousFunction(*this) );
1105           escriptDataC _in2 = temp.getDataC();           //escriptDataC _in2 = temp.getDataC();
1106        } else {        } else {
1107           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);           Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&_in,&_target);
1108        }        }
# Line 1487  void MeshAdapter::setNewX(const escript: Line 1512  void MeshAdapter::setNewX(const escript:
1512         tmp = new_x.getDataC();         tmp = new_x.getDataC();
1513         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);
1514     } else {     } else {
1515         escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );         throw FinleyAdapterException("As of version escript3.3 - SetNewX only accepts ContinuousFunction arguments please interpolate.");    
1516    /*       escript::Data new_x_inter=escript::Data( new_x,  continuousFunction(*this) );
1517         tmp = new_x_inter.getDataC();         tmp = new_x_inter.getDataC();
1518         Finley_Mesh_setCoordinates(mesh,&tmp);         Finley_Mesh_setCoordinates(mesh,&tmp);*/
    }  
    checkFinleyError();  
 }  
   
 //  
 // Helper for the save* methods. Extracts optional data variable names and the  
 // corresponding pointers from python dictionary. Caller must free arrays.  
 //  
 void MeshAdapter::extractArgsFromDict(const boost::python::dict& arg, int& numData, char**& names, escriptDataC*& data, escriptDataC**& dataPtr) const  
 {  
    numData = boost::python::extract<int>(arg.attr("__len__")());  
    /* win32 refactor */  
    names = (numData>0) ? TMPMEMALLOC(numData, char*) : (char**)NULL;  
    data = (numData>0) ? TMPMEMALLOC(numData,escriptDataC) : (escriptDataC*)NULL;  
    dataPtr = (numData>0) ? TMPMEMALLOC(numData,escriptDataC*) : (escriptDataC**)NULL;  
   
    boost::python::list keys=arg.keys();  
    for (int i=0; i<numData; ++i) {  
       string n=boost::python::extract<string>(keys[i]);  
       escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);  
       if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)  
          throw FinleyAdapterException("Error: Data must be defined on same Domain");  
       data[i] = d.getDataC();  
       dataPtr[i] = &(data[i]);  
       names[i] = TMPMEMALLOC(n.length()+1, char);  
       strcpy(names[i], n.c_str());  
    }  
 }  
   
 //  
 // saves mesh and optionally data arrays in openDX format  
 //  
 void MeshAdapter::saveDX(const string& filename,const boost::python::dict& arg) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    escriptDataC **ptr_data;  
   
    extractArgsFromDict(arg, num_data, names, data, ptr_data);  
    Finley_Mesh_saveDX(filename.c_str(), m_finleyMesh.get(), num_data, names, ptr_data);  
    checkFinleyError();  
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
1519     }     }
    TMPMEMFREE(names);  
   
    return;  
 }  
   
 //  
 // saves mesh and optionally data arrays in VTK format  
 //  
 void MeshAdapter::saveVTK(const string& filename,const boost::python::dict& arg,  const string& metadata, const string& metadata_schema) const  
 {  
    int num_data;  
    char **names;  
    escriptDataC *data;  
    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());  
1520     checkFinleyError();     checkFinleyError();
   
    /* win32 refactor */  
    TMPMEMFREE(data);  
    TMPMEMFREE(ptr_data);  
    for(int i=0; i<num_data; i++)  
    {  
       TMPMEMFREE(names[i]);  
    }  
    TMPMEMFREE(names);  
1521  }  }
1522    
1523  bool MeshAdapter::ownSample(int fs_code, index_t id) const  bool MeshAdapter::ownSample(int fs_code, index_t id) const
1524  {  {
1525        if (getMPISize() > 1) {
1526  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1527      index_t myFirstNode=0, myLastNode=0, k=0;          index_t myFirstNode=0, myLastNode=0, k=0;
1528      index_t* globalNodeIndex=0;          index_t* globalNodeIndex=0;
1529      Finley_Mesh* mesh_p=m_finleyMesh.get();          Finley_Mesh* mesh_p=m_finleyMesh.get();
1530      if (fs_code == FINLEY_REDUCED_NODES)          /*
1531      {           * this method is only used by saveDataCSV which would use the returned
1532      myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);           * values for reduced nodes wrongly so this case is disabled for now
1533      myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);          if (fs_code == FINLEY_REDUCED_NODES)
1534      globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);          {
1535      }              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);
1536      else              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);
1537      {              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1538      myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);          }
1539      myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);          else
1540      globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);          */
1541      }          if (fs_code == FINLEY_NODES)
1542      k=globalNodeIndex[id];          {
1543      return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );              myFirstNode = Finley_NodeFile_getFirstNode(mesh_p->Nodes);
1544                myLastNode = Finley_NodeFile_getLastNode(mesh_p->Nodes);
1545                globalNodeIndex = Finley_NodeFile_borrowGlobalNodesIndex(mesh_p->Nodes);
1546            }
1547            else
1548            {
1549                throw FinleyAdapterException("Unsupported function space type for ownSample()");
1550            }
1551    
1552            k=globalNodeIndex[id];
1553            return static_cast<bool>( (myFirstNode <= k) && (k < myLastNode) );
1554  #endif  #endif
1555        }
1556      return true;      return true;
1557  }  }
1558    
1559    
   
1560  //  //
1561  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros
1562  //  //
# Line 1613  ASM_ptr MeshAdapter::newSystemMatrix( Line 1575  ASM_ptr MeshAdapter::newSystemMatrix(
1575        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of row function space does not match the domain of matrix generator.");
1576     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1577     if (col_domain!=*this)     if (col_domain!=*this)
1578        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of column function space does not match the domain of matrix generator.");
1579     // is the function space type right     // is the function space type right
1580     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {     if (row_functionspace.getTypeCode()==DegreesOfFreedom) {
1581        reduceRowOrder=0;        reduceRowOrder=0;
# Line 1654  ASM_ptr MeshAdapter::newSystemMatrix( Line 1616  ASM_ptr MeshAdapter::newSystemMatrix(
1616  // creates a TransportProblemAdapter  // creates a TransportProblemAdapter
1617  //  //
1618  ATP_ptr MeshAdapter::newTransportProblem(  ATP_ptr MeshAdapter::newTransportProblem(
                                                          const bool useBackwardEuler,  
1619                                                           const int blocksize,                                                           const int blocksize,
1620                                                           const escript::FunctionSpace& functionspace,                                                           const escript::FunctionSpace& functionspace,
1621                                                           const int type) const                                                           const int type) const
# Line 1677  ATP_ptr MeshAdapter::newTransportProblem Line 1638  ATP_ptr MeshAdapter::newTransportProblem
1638     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);     Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(getFinley_Mesh(),reduceOrder,reduceOrder);
1639     checkFinleyError();     checkFinleyError();
1640     Paso_TransportProblem* transportProblem;     Paso_TransportProblem* transportProblem;
1641     transportProblem=Paso_TransportProblem_alloc(useBackwardEuler,fsystemMatrixPattern,blocksize);     transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern,blocksize);
1642     checkPasoError();     checkPasoError();
1643     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);     Paso_SystemMatrixPattern_free(fsystemMatrixPattern);
1644     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);     TransportProblemAdapter* tpa=new TransportProblemAdapter(transportProblem,blocksize,functionspace);
1645     return ATP_ptr(tpa);     return ATP_ptr(tpa);
1646  //   return TransportProblemAdapter(transportProblem,useBackwardEuler,blocksize,functionspace);  //   return TransportProblemAdapter(transportProblem, blocksize,functionspace);
1647  }  }
1648    
1649  //  //
# Line 2018  bool MeshAdapter::probeInterpolationOnDo Line 1979  bool MeshAdapter::probeInterpolationOnDo
1979     return false;     return false;
1980  }  }
1981    
1982    signed char MeshAdapter::preferredInterpolationOnDomain(int functionSpaceType_source,int functionSpaceType_target) const
1983    {
1984        if (probeInterpolationOnDomain(functionSpaceType_source, functionSpaceType_target))
1985        {
1986            return 1;
1987        }
1988        else if (probeInterpolationOnDomain(functionSpaceType_target, functionSpaceType_source))
1989        {
1990            return -1;
1991        }
1992        return 0;
1993    }  
1994      
1995  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const  bool MeshAdapter::probeInterpolationACross(int functionSpaceType_source,const AbstractDomain& targetDomain, int functionSpaceType_target) const
1996  {  {
1997     return false;     return false;
# Line 2041  bool MeshAdapter::operator!=(const Abstr Line 2015  bool MeshAdapter::operator!=(const Abstr
2015  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
2016  {  {
2017     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2018     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
2019     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2020  }  }
2021    
2022  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
2023  {  {
2024     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2025     int out=Paso_TransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0, mesh->MPIInfo);     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
2026     checkPasoError();             package, symmetry, mesh->MPIInfo);
    return out;  
2027  }  }
2028    
2029  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
# Line 2231  void MeshAdapter::setTagMap(const string Line 2204  void MeshAdapter::setTagMap(const string
2204  {  {
2205     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2206     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);     Finley_Mesh_addTagMap(mesh, name.c_str(),tag);
2207     checkPasoError();     checkFinleyError();
2208     // throwStandardException("MeshAdapter::set TagMap is not implemented.");     // throwStandardException("MeshAdapter::set TagMap is not implemented.");
2209  }  }
2210    
# Line 2240  int MeshAdapter::getTag(const string& na Line 2213  int MeshAdapter::getTag(const string& na
2213     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
2214     int tag=0;     int tag=0;
2215     tag=Finley_Mesh_getTag(mesh, name.c_str());     tag=Finley_Mesh_getTag(mesh, name.c_str());
2216     checkPasoError();     checkFinleyError();
2217     // throwStandardException("MeshAdapter::getTag is not implemented.");     // throwStandardException("MeshAdapter::getTag is not implemented.");
2218     return tag;     return tag;
2219  }  }
# Line 2428  ReferenceElementSetWrapper::~ReferenceEl Line 2401  ReferenceElementSetWrapper::~ReferenceEl
2401    Finley_ReferenceElementSet_dealloc(m_refSet);    Finley_ReferenceElementSet_dealloc(m_refSet);
2402  }  }
2403    
2404    // points will be flattened
2405    void MeshAdapter:: addDiracPoints(const std::vector<double>& points, const std::vector<int>& tags) const
2406    {
2407          const int dim = getDim();
2408          int numPoints=points.size()/dim;
2409          int numTags=tags.size();
2410          Finley_Mesh* mesh=m_finleyMesh.get();
2411          
2412          if ( points.size() % dim != 0 )
2413          {
2414        throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
2415          }
2416          
2417          if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2418         throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2419          
2420          double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2421          int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2422          
2423          for (int i=0;i<numPoints;++i) {
2424           points_ptr[ i * dim     ] = points[i * dim ];
2425           if ( dim > 1 ) points_ptr[ i * dim + 1 ] = points[i * dim + 1];
2426           if ( dim > 2 ) points_ptr[ i * dim + 2 ] = points[i * dim + 2];    
2427               tags_ptr[i]=tags[i];
2428          }
2429          
2430          Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2431          checkFinleyError();
2432          
2433          TMPMEMFREE(points_ptr);
2434          TMPMEMFREE(tags_ptr);
2435    }
2436    
2437    
2438    // void MeshAdapter:: addDiracPoints(const boost::python::list& points, const boost::python::list& tags) const
2439    // {
2440    //       const int dim = getDim();
2441    //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2442    //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2443    //       Finley_Mesh* mesh=m_finleyMesh.get();
2444    //      
2445    //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2446    //   throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2447    //      
2448    //       double* points_ptr=TMPMEMALLOC(numPoints * dim, double);
2449    //       int*    tags_ptr= TMPMEMALLOC(numPoints, int);
2450    //      
2451    //       for (int i=0;i<numPoints;++i) {
2452    //     int tag_id=-1;
2453    //     int numComps=boost::python::extract<int>(points[i].attr("__len__")());
2454    //     if  ( numComps !=   dim ) {
2455    //                stringstream temp;            
2456    //                temp << "Error - illegal number of components " << numComps << " for point " << i;              
2457    //                throw FinleyAdapterException(temp.str());
2458    //     }
2459    //     points_ptr[ i * dim     ] = boost::python::extract<double>(points[i][0]);
2460    //     if ( dim > 1 ) points_ptr[ i * dim + 1 ] = boost::python::extract<double>(points[i][1]);
2461    //     if ( dim > 2 ) points_ptr[ i * dim + 2 ] = boost::python::extract<double>(points[i][2]);
2462    //    
2463    //     if ( numTags > 0) {
2464    //            boost::python::extract<string> ex_str(tags[i]);
2465    //        if  ( ex_str.check() ) {
2466    //            tag_id=getTag( ex_str());
2467    //        } else {
2468    //             boost::python::extract<int> ex_int(tags[i]);
2469    //             if ( ex_int.check() ) {
2470    //                 tag_id=ex_int();
2471    //             } else {
2472    //              stringstream temp;          
2473    //                  temp << "Error - unable to extract tag for point " << i;
2474    //              throw FinleyAdapterException(temp.str());
2475    //            }
2476    //        }
2477    //     }      
2478    //            tags_ptr[i]=tag_id;
2479    //       }
2480    //      
2481    //       Finley_Mesh_addPoints(mesh, numPoints, points_ptr, tags_ptr);
2482    //       checkPasoError();
2483    //      
2484    //       TMPMEMFREE(points_ptr);
2485    //       TMPMEMFREE(tags_ptr);
2486    // }
2487    
2488    /*
2489    void MeshAdapter:: addDiracPoint( const boost::python::list& point, const int tag) const
2490    {  
2491        boost::python::list points =  boost::python::list();
2492        boost::python::list tags =  boost::python::list();
2493        points.append(point);
2494        tags.append(tag);
2495        addDiracPoints(points, tags);
2496    }
2497    */
2498    
2499    /*
2500    void MeshAdapter:: addDiracPointWithTagName( const boost::python::list& point, const std::string& tag) const
2501    {
2502            boost::python::list points =   boost::python::list();
2503            boost::python::list tags =   boost::python::list();
2504            points.append(point);
2505            tags.append(tag);
2506            addDiracPoints(points, tags);
2507    }
2508    */
2509  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.3269  
changed lines
  Added in v.4346

  ViewVC Help
Powered by ViewVC 1.1.26