/[escript]/branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapter.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/CPPAdapter/MeshAdapter.cpp

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

revision 5932 by caltinay, Fri Feb 5 03:37:49 2016 UTC revision 5933 by caltinay, Wed Feb 17 23:53:30 2016 UTC
# Line 18  Line 18 
18  #include "esysUtils/first.h"  #include "esysUtils/first.h"
19    
20  #include "MeshAdapter.h"  #include "MeshAdapter.h"
 #include "escript/Data.h"  
 #include "escript/DataFactory.h"  
21  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
22  #include "esysUtils/EsysRandom.h"  #include "esysUtils/EsysRandom.h"
23    #include <escript/Data.h>
24    #include <escript/DataFactory.h>
25    #include <escript/SolverOptions.h>
26    
27    #include <paso/SystemMatrix.h>
28    #include <paso/Transport.h>
29    
30  #ifdef USE_NETCDF  #ifdef USE_NETCDF
31  #include <netcdfcpp.h>  #include <netcdfcpp.h>
# Line 30  Line 34 
34    
35  using namespace std;  using namespace std;
36  using namespace escript;  using namespace escript;
 using namespace paso;  
37  namespace bp = boost::python;  namespace bp = boost::python;
38    
39  namespace dudley {  namespace dudley {
# Line 129  void MeshAdapter::dump(const string& fil Line 132  void MeshAdapter::dump(const string& fil
132     Dudley_Mesh *mesh = m_dudleyMesh.get();     Dudley_Mesh *mesh = m_dudleyMesh.get();
133     Dudley_TagMap* tag_map;     Dudley_TagMap* tag_map;
134     int num_Tags = 0;     int num_Tags = 0;
135     int mpi_size             = mesh->MPIInfo->size;     int mpi_size                         = mesh->MPIInfo->size;
136     int mpi_rank             = mesh->MPIInfo->rank;     int mpi_rank                         = mesh->MPIInfo->rank;
137     int numDim               = mesh->Nodes->numDim;     int numDim                           = mesh->Nodes->numDim;
138     int numNodes             = mesh->Nodes->numNodes;     int numNodes                         = mesh->Nodes->numNodes;
139     int num_Elements         = mesh->Elements->numElements;     int num_Elements                     = mesh->Elements->numElements;
140     int num_FaceElements         = mesh->FaceElements->numElements;     int num_FaceElements                 = mesh->FaceElements->numElements;
141     int num_Points           = mesh->Points->numElements;     int num_Points                       = mesh->Points->numElements;
142     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes            = mesh->Elements->numNodes;
143     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes        = mesh->FaceElements->numNodes;
144  #ifdef ESYS_MPI  #ifdef ESYS_MPI
145     MPI_Status status;     MPI_Status status;
146  #endif  #endif
# Line 474  void MeshAdapter::dump(const string& fil Line 477  void MeshAdapter::dump(const string& fil
477    
478  #else  #else
479     Dudley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Dudley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
480  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
481     checkDudleyError();     checkDudleyError();
482  }  }
483    
# Line 676  pair<int,int> MeshAdapter::getDataShape( Line 679  pair<int,int> MeshAdapter::getDataShape(
679  //  //
680  // 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:
681  //  //
682  void MeshAdapter::addPDEToSystem(  void MeshAdapter::addPDEToSystem(AbstractSystemMatrix& mat, escript::Data& rhs,
683                                   AbstractSystemMatrix& mat, escript::Data& rhs,                                   const escript::Data& A, const escript::Data& B,
684                                   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& C, const escript::Data& D,
685                                     const escript::Data& X, const escript::Data& Y,
686                                   const escript::Data& d, const escript::Data& y,                                   const escript::Data& d, const escript::Data& y,
687                   const escript::Data& d_contact, const escript::Data& y_contact,                                   const escript::Data& d_contact,
688                                   const escript::Data& d_dirac,const escript::Data& y_dirac) const                                   const escript::Data& y_contact,
689  {                                   const escript::Data& d_dirac,
690      if (!d_contact.isEmpty())                                   const escript::Data& y_dirac) const
691      {  {
692      throw DudleyAdapterException("Dudley does not support d_contact");      if (!d_contact.isEmpty() || !y_contact.isEmpty())
693      }          throw DudleyAdapterException("Dudley does not support contact elements");
694      if (!y_contact.isEmpty())  
695      {      paso::SystemMatrix* smat = dynamic_cast<paso::SystemMatrix*>(&mat);
696      throw DudleyAdapterException("Dudley does not support y_contact");      if (smat) {
697            paso::SystemMatrix_ptr S(smat->shared_from_this());
698            Dudley_Mesh* mesh = m_dudleyMesh.get();
699    
700            Dudley_Assemble_PDE(mesh->Nodes, mesh->Elements, S, &rhs, &A, &B, &C, &D, &X, &Y);
701            checkDudleyError();
702    
703            Dudley_Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, &rhs, 0, 0, 0, &d, 0, &y);
704            checkDudleyError();
705    
706            Dudley_Assemble_PDE(mesh->Nodes,mesh->Points, S, &rhs, 0, 0, 0, &d_dirac, 0, &y_dirac);
707            checkDudleyError();
708            return;
709      }      }
710     SystemMatrixAdapter* smat=dynamic_cast<SystemMatrixAdapter*>(&mat);      throw DudleyAdapterException("Dudley only accepts Paso system matrices");
    if (smat==0)  
    {  
     throw DudleyAdapterException("Dudley only accepts Paso system matrices");  
    }  
   
    Dudley_Mesh* mesh=m_dudleyMesh.get();  
   
    Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,smat->getPaso_SystemMatrix(), &rhs, &A, &B, &C, &D, &X, &Y );  
    checkDudleyError();  
   
    Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, smat->getPaso_SystemMatrix(), &rhs, 0, 0, 0, &d, 0, &y );  
    checkDudleyError();  
   
    Dudley_Assemble_PDE(mesh->Nodes,mesh->Points, smat->getPaso_SystemMatrix(), &rhs, 0, 0, 0, &d_dirac, 0, &y_dirac );  
    checkDudleyError();  
711  }  }
712    
713  void  MeshAdapter::addPDEToLumpedSystem(  void  MeshAdapter::addPDEToLumpedSystem(escript::Data& mat,
                                         escript::Data& mat,  
714                                          const escript::Data& D,                                          const escript::Data& D,
715                                          const escript::Data& d,                                          const escript::Data& d,
716                                          const escript::Data& d_dirac,                                          const escript::Data& d_dirac,
717                      const bool useHRZ) const                                          bool useHRZ) const
718  {  {
719     Dudley_Mesh* mesh=m_dudleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
720    
721     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements,&mat, &D, useHRZ);     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->Elements, &mat, &D, useHRZ);
722     checkDudleyError();     checkDudleyError();
723        
724     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&mat, &d, useHRZ);     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements, &mat, &d, useHRZ);
725     checkDudleyError();     checkDudleyError();
726    
727     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements,&mat, &d_dirac, useHRZ);     Dudley_Assemble_LumpedSystem(mesh->Nodes,mesh->FaceElements, &mat, &d_dirac, useHRZ);
728     checkDudleyError();     checkDudleyError();
729    
730  }  }
# Line 733  void  MeshAdapter::addPDEToLumpedSystem( Line 733  void  MeshAdapter::addPDEToLumpedSystem(
733  //  //
734  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
735  //  //
736  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  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
737  {  {
738     if (!y_contact.isEmpty())     if (!y_contact.isEmpty())
739     {     {
740      throw DudleyAdapterException("Dudley does not support y_contact");          throw DudleyAdapterException("Dudley does not support y_contact");
741     }     }
742     Dudley_Mesh* mesh=m_dudleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
743    
744     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, paso::SystemMatrix_ptr(), &rhs, 0, 0, 0, 0, &X, &Y);     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements, escript::ASM_ptr(), &rhs,
745                           NULL, NULL, NULL, NULL, &X, &Y);
746     checkDudleyError();     checkDudleyError();
747    
748     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, paso::SystemMatrix_ptr(), &rhs, 0, 0, 0, 0, 0, &y );     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, escript::ASM_ptr(),
749                           &rhs, NULL, NULL, NULL, NULL, NULL, &y);
750     checkDudleyError();     checkDudleyError();
751    
752     Dudley_Assemble_PDE(mesh->Nodes,mesh->Points, paso::SystemMatrix_ptr(), &rhs, 0, 0, 0, 0, 0, &y_dirac );     Dudley_Assemble_PDE(mesh->Nodes,mesh->Points, escript::ASM_ptr(), &rhs,
753                           NULL, NULL, NULL, NULL, NULL, &y_dirac);
754     checkDudleyError();     checkDudleyError();
755  }  }
756  //  //
757  // adds PDE of second order into a transport problem  // adds PDE of second order into a transport problem
758  //  //
759  void MeshAdapter::addPDEToTransportProblem(  void MeshAdapter::addPDEToTransportProblem(
760                                             AbstractTransportProblem& tp, escript::Data& source, const escript::Data& M,          AbstractTransportProblem& tp, escript::Data& source,
761                                             const escript::Data& A, const escript::Data& B, const escript::Data& C,          const escript::Data& M, const escript::Data& A, const escript::Data& B,
762                                             const  escript::Data& D,const  escript::Data& X,const  escript::Data& Y,          const escript::Data& C, const  escript::Data& D,const escript::Data& X,
763                                             const escript::Data& d, const escript::Data& y,          const escript::Data& Y, const escript::Data& d, const escript::Data& y,
764                         const escript::Data& d_contact,const escript::Data& y_contact,          const escript::Data& d_contact,const escript::Data& y_contact,
765                         const escript::Data& d_dirac, const escript::Data& y_dirac) const          const escript::Data& d_dirac, const escript::Data& y_dirac) const
766  {  {
767      if (!d_contact.isEmpty())      if (!d_contact.isEmpty())
768      {      {
769      throw DudleyAdapterException("Dudley does not support d_contact");          throw DudleyAdapterException("Dudley does not support d_contact");
770      }      }
771      if (!y_contact.isEmpty())      if (!y_contact.isEmpty())
772      {      {
773      throw DudleyAdapterException("Dudley does not support y_contact");          throw DudleyAdapterException("Dudley does not support y_contact");
774      }        }  
775     TransportProblemAdapter* tpa=dynamic_cast<TransportProblemAdapter*>(&tp);      paso::TransportProblem* ptp = dynamic_cast<paso::TransportProblem*>(&tp);
776     if (tpa==0)      if (!ptp)
777     {      {
778      throw DudleyAdapterException("Dudley only accepts Paso transport problems");          throw DudleyAdapterException("Dudley only accepts Paso transport problems");
779     }      }
780     DataTypes::ShapeType shape;      DataTypes::ShapeType shape;
781     source.expand();      source.expand();
782    
783     Dudley_Mesh* mesh=m_dudleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
    paso::TransportProblem_ptr _tp(tpa->getPaso_TransportProblem());  
784    
785     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->mass_matrix, &source, 0, 0, 0, &M, 0, 0 );     Dudley_Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowMassMatrix(),
786                           &source, 0, 0, 0, &M, 0, 0);
787     checkDudleyError();     checkDudleyError();
788    
789     Dudley_Assemble_PDE(mesh->Nodes,mesh->Elements,_tp->transport_matrix, &source, &A, &B, &C, &D, &X, &Y );     Dudley_Assemble_PDE(mesh->Nodes, mesh->Elements, ptp->borrowTransportMatrix(),
790                           &source, &A, &B, &C, &D, &X, &Y);
791     checkDudleyError();     checkDudleyError();
792    
793     Dudley_Assemble_PDE(mesh->Nodes,mesh->FaceElements, _tp->transport_matrix, &source, 0, 0, 0, &d, 0, &y );     Dudley_Assemble_PDE(mesh->Nodes, mesh->FaceElements,
794                           ptp->borrowTransportMatrix(), &source, NULL, NULL, NULL,
795                           &d, NULL, &y);
796     checkDudleyError();     checkDudleyError();
797    
798     Dudley_Assemble_PDE(mesh->Nodes,mesh->Points, _tp->transport_matrix, &source, 0, 0, 0, &d_dirac, 0, &y_dirac );     Dudley_Assemble_PDE(mesh->Nodes, mesh->Points, ptp->borrowTransportMatrix(),
799                           &source, NULL, NULL, NULL, &d_dirac, NULL, &y_dirac);
800     checkDudleyError();     checkDudleyError();
801  }  }
802    
# Line 1305  bool MeshAdapter::ownSample(int fs_code, Line 1312  bool MeshAdapter::ownSample(int fs_code,
1312    
1313    
1314  //  //
1315  // creates a SystemMatrixAdapter stiffness matrix an initializes it with zeros  // creates a stiffness matrix an initializes it with zeros
1316  //  //
1317  ASM_ptr MeshAdapter::newSystemMatrix(  ASM_ptr MeshAdapter::newSystemMatrix(const int row_blocksize,
1318                                                   const int row_blocksize,                                       const escript::FunctionSpace& row_functionspace,
1319                                                   const escript::FunctionSpace& row_functionspace,                                       const int column_blocksize,
1320                                                   const int column_blocksize,                                       const escript::FunctionSpace& column_functionspace,
1321                                                   const escript::FunctionSpace& column_functionspace,                                       const int type) const
                                                  const int type) const  
1322  {  {
1323     int reduceRowOrder=0;     int reduceRowOrder=0;
1324     int reduceColOrder=0;     int reduceColOrder=0;
# Line 1342  ASM_ptr MeshAdapter::newSystemMatrix( Line 1348  ASM_ptr MeshAdapter::newSystemMatrix(
1348    
1349     paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder));     paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(getDudley_Mesh(),reduceRowOrder,reduceColOrder));
1350     checkDudleyError();     checkDudleyError();
1351     paso::SystemMatrix_ptr fsystemMatrix;     paso::SystemMatrix_ptr sm;
1352     int trilinos = 0;     int trilinos = 0;
1353     if (trilinos) {     if (trilinos) {
1354  #ifdef TRILINOS  #ifdef TRILINOS
# Line 1350  ASM_ptr MeshAdapter::newSystemMatrix( Line 1356  ASM_ptr MeshAdapter::newSystemMatrix(
1356  #endif  #endif
1357     }     }
1358     else {     else {
1359        fsystemMatrix.reset(new paso::SystemMatrix(type, fsystemMatrixPattern, row_blocksize, column_blocksize,false));        sm.reset(new paso::SystemMatrix(type, fsystemMatrixPattern,
1360                      row_blocksize, column_blocksize, false, row_functionspace,
1361                      column_functionspace));
1362     }     }
1363     checkPasoError();     checkPasoError();
1364     SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix,row_blocksize,row_functionspace, column_blocksize,column_functionspace);     return sm;
    return ASM_ptr(sma);  
1365  }  }
1366    
1367  //  //
1368  // creates a TransportProblemAdapter  // creates a TransportProblem
1369  //  //
1370  ATP_ptr MeshAdapter::newTransportProblem(const int blocksize,  ATP_ptr MeshAdapter::newTransportProblem(int blocksize,
1371                                           const escript::FunctionSpace& fs,                                           const escript::FunctionSpace& fs,
1372                                           const int type) const                                           int type) const
1373  {  {
1374     int reduceOrder=0;     int reduceOrder=0;
1375     // is the domain right?     // is the domain right?
# Line 1379  ATP_ptr MeshAdapter::newTransportProblem Line 1386  ATP_ptr MeshAdapter::newTransportProblem
1386     }     }
1387     // generate matrix:     // generate matrix:
1388    
1389     paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(getDudley_Mesh(),reduceOrder,reduceOrder));     paso::SystemMatrixPattern_ptr fsystemMatrixPattern(Dudley_getPattern(
1390                   getDudley_Mesh(),reduceOrder,reduceOrder));
1391     checkDudleyError();     checkDudleyError();
1392     paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(     paso::TransportProblem_ptr transportProblem(new paso::TransportProblem(
1393                                              fsystemMatrixPattern, blocksize));                                              fsystemMatrixPattern, blocksize,
1394                                                fs));
1395     checkPasoError();     checkPasoError();
1396     AbstractTransportProblem* atp=new TransportProblemAdapter(transportProblem, blocksize, fs);     return transportProblem;
    return ATP_ptr(atp);  
1397  }  }
1398    
1399  //  //
# Line 1424  bool Line 1432  bool
1432  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const  MeshAdapter::commonFunctionSpace(const vector<int>& fs, int& resultcode) const
1433  {  {
1434     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]     /* The idea is to use equivalence classes. [Types which can be interpolated back and forth]
1435      class 1: DOF <-> Nodes          class 1: DOF <-> Nodes
1436      class 2: ReducedDOF <-> ReducedNodes          class 2: ReducedDOF <-> ReducedNodes
1437      class 3: Points          class 3: Points
1438      class 4: Elements          class 4: Elements
1439      class 5: ReducedElements          class 5: ReducedElements
1440      class 6: FaceElements          class 6: FaceElements
1441      class 7: ReducedFaceElements          class 7: ReducedFaceElements
1442      class 8: ContactElementZero <-> ContactElementOne          class 8: ContactElementZero <-> ContactElementOne
1443      class 9: ReducedContactElementZero <-> ReducedContactElementOne          class 9: ReducedContactElementZero <-> ReducedContactElementOne
1444    
1445     There is also a set of lines. Interpolation is possible down a line but not between lines.     There is also a set of lines. Interpolation is possible down a line but not between lines.
1446     class 1 and 2 belong to all lines so aren't considered.     class 1 and 2 belong to all lines so aren't considered.
1447      line 0: class 3          line 0: class 3
1448      line 1: class 4,5          line 1: class 4,5
1449      line 2: class 6,7          line 2: class 6,7
1450      line 3: class 8,9          line 3: class 8,9
1451    
1452     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.
1453     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.
# Line 1449  MeshAdapter::commonFunctionSpace(const v Line 1457  MeshAdapter::commonFunctionSpace(const v
1457          return false;          return false;
1458      }      }
1459      vector<int> hasclass(10);      vector<int> hasclass(10);
1460      vector<int> hasline(4);      vector<int> hasline(4);    
1461      bool hasnodes=false;      bool hasnodes=false;
1462      bool hasrednodes=false;      bool hasrednodes=false;
1463      for (int i=0;i<fs.size();++i)      for (int i=0;i<fs.size();++i)
1464      {      {
1465      switch(fs[i])          switch(fs[i])
1466      {          {
1467      case(Nodes):   hasnodes=true;   // no break is deliberate          case(Nodes):   hasnodes=true;   // no break is deliberate
1468      case(DegreesOfFreedom):          case(DegreesOfFreedom):
1469          hasclass[1]=1;                  hasclass[1]=1;
1470          break;                  break;
1471      case(ReducedNodes):    hasrednodes=true;    // no break is deliberate          case(ReducedNodes):    hasrednodes=true;        // no break is deliberate
1472      case(ReducedDegreesOfFreedom):          case(ReducedDegreesOfFreedom):
1473          hasclass[2]=1;                  hasclass[2]=1;
1474          break;                  break;
1475      case(Points):          case(Points):
1476          hasline[0]=1;                  hasline[0]=1;
1477          hasclass[3]=1;                  hasclass[3]=1;
1478          break;                  break;
1479      case(Elements):          case(Elements):
1480          hasclass[4]=1;                  hasclass[4]=1;
1481          hasline[1]=1;                  hasline[1]=1;
1482          break;                  break;
1483      case(ReducedElements):          case(ReducedElements):
1484          hasclass[5]=1;                  hasclass[5]=1;
1485          hasline[1]=1;                  hasline[1]=1;
1486          break;                  break;
1487      case(FaceElements):          case(FaceElements):
1488          hasclass[6]=1;                  hasclass[6]=1;
1489          hasline[2]=1;                  hasline[2]=1;
1490          break;                  break;
1491      case(ReducedFaceElements):          case(ReducedFaceElements):
1492          hasclass[7]=1;                  hasclass[7]=1;
1493          hasline[2]=1;                  hasline[2]=1;
1494          break;                  break;
1495      default:          default:
1496          return false;                  return false;
1497      }          }
1498      }      }
1499      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];
1500      // fail if we have more than one leaf group      // fail if we have more than one leaf group
1501    
1502      if (totlines>1)      if (totlines>1)
1503      {      {
1504      return false;   // there are at least two branches we can't interpolate between          return false;   // there are at least two branches we can't interpolate between
1505      }      }
1506      else if (totlines==1)      else if (totlines==1)
1507      {      {
1508      if (hasline[0]==1)      // we have points          if (hasline[0]==1)              // we have points
1509      {          {
1510          resultcode=Points;              resultcode=Points;
1511      }          }
1512      else if (hasline[1]==1)          else if (hasline[1]==1)
1513      {          {
1514          if (hasclass[5]==1)              if (hasclass[5]==1)
1515          {              {
1516          resultcode=ReducedElements;                  resultcode=ReducedElements;
1517          }              }
1518          else              else
1519          {              {
1520          resultcode=Elements;                  resultcode=Elements;
1521          }              }
1522      }          }
1523      else if (hasline[2]==1)          else if (hasline[2]==1)
1524      {          {
1525          if (hasclass[7]==1)              if (hasclass[7]==1)
1526          {              {
1527          resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
1528          }              }
1529          else              else
1530          {              {
1531          resultcode=FaceElements;                  resultcode=FaceElements;
1532          }              }
1533      }          }
1534      else    // so we must be in line3          else    // so we must be in line3
1535      {          {
1536    
1537          throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");              throw DudleyAdapterException("Programmer Error - choosing between contact elements - we should never get here.");
1538    
1539      }          }
1540      }      }
1541      else    // totlines==0      else        // totlines==0
1542      {      {
1543      if (hasclass[2]==1)          if (hasclass[2]==1)
1544      {          {
1545          // something from class 2                  // something from class 2
1546          resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);                  resultcode=(hasrednodes?ReducedNodes:ReducedDegreesOfFreedom);
1547      }          }
1548      else          else
1549      {   // something from class 1          {       // something from class 1
1550          resultcode=(hasnodes?Nodes:DegreesOfFreedom);                  resultcode=(hasnodes?Nodes:DegreesOfFreedom);
1551      }          }
1552      }      }
1553      return true;      return true;
1554  }  }
# Line 1564  bool MeshAdapter::probeInterpolationOnDo Line 1572  bool MeshAdapter::probeInterpolationOnDo
1572  {  {
1573     switch(functionSpaceType_source) {     switch(functionSpaceType_source) {
1574     case(Nodes):     case(Nodes):
1575      switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1576      case(Nodes):          case(Nodes):
1577      case(ReducedNodes):          case(ReducedNodes):
1578      case(ReducedDegreesOfFreedom):          case(ReducedDegreesOfFreedom):
1579      case(DegreesOfFreedom):          case(DegreesOfFreedom):
1580      case(Elements):          case(Elements):
1581      case(ReducedElements):          case(ReducedElements):
1582      case(FaceElements):          case(FaceElements):
1583      case(ReducedFaceElements):          case(ReducedFaceElements):
1584      case(Points):          case(Points):
1585      return true;          return true;
1586      default:          default:
1587            stringstream temp;                stringstream temp;
1588            temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;                temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1589            throw DudleyAdapterException(temp.str());                throw DudleyAdapterException(temp.str());
1590     }     }
1591     break;     break;
1592     case(ReducedNodes):     case(ReducedNodes):
1593      switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1594      case(ReducedNodes):          case(ReducedNodes):
1595      case(ReducedDegreesOfFreedom):          case(ReducedDegreesOfFreedom):
1596      case(Elements):          case(Elements):
1597      case(ReducedElements):          case(ReducedElements):
1598      case(FaceElements):          case(FaceElements):
1599      case(ReducedFaceElements):          case(ReducedFaceElements):
1600      case(Points):          case(Points):
1601      return true;          return true;
1602      case(Nodes):          case(Nodes):
1603      case(DegreesOfFreedom):          case(DegreesOfFreedom):
1604      return false;          return false;
1605      default:          default:
1606          stringstream temp;                  stringstream temp;
1607          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;                  temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1608          throw DudleyAdapterException(temp.str());                  throw DudleyAdapterException(temp.str());
1609     }     }
1610     break;     break;
1611     case(Elements):     case(Elements):
1612      if (functionSpaceType_target==Elements) {          if (functionSpaceType_target==Elements) {
1613        return true;            return true;
1614      } else if (functionSpaceType_target==ReducedElements) {          } else if (functionSpaceType_target==ReducedElements) {
1615        return true;            return true;
1616          } else {          } else {
1617            return false;            return false;
1618          }          }
1619     case(ReducedElements):     case(ReducedElements):
1620      if (functionSpaceType_target==ReducedElements) {          if (functionSpaceType_target==ReducedElements) {
1621        return true;            return true;
1622      } else {          } else {
1623            return false;            return false;
1624      }          }
1625     case(FaceElements):     case(FaceElements):
1626      if (functionSpaceType_target==FaceElements) {          if (functionSpaceType_target==FaceElements) {
1627              return true;                  return true;
1628      } else if (functionSpaceType_target==ReducedFaceElements) {          } else if (functionSpaceType_target==ReducedFaceElements) {
1629              return true;                  return true;
1630      } else {          } else {
1631              return false;                  return false;
1632      }          }
1633     case(ReducedFaceElements):     case(ReducedFaceElements):
1634      if (functionSpaceType_target==ReducedFaceElements) {          if (functionSpaceType_target==ReducedFaceElements) {
1635              return true;                  return true;
1636      } else {          } else {
1637          return false;                  return false;
1638      }          }
1639     case(Points):     case(Points):
1640      if (functionSpaceType_target==Points) {          if (functionSpaceType_target==Points) {
1641              return true;                  return true;
1642      } else {          } else {
1643              return false;                  return false;
1644      }          }
1645     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1646      switch(functionSpaceType_target) {          switch(functionSpaceType_target) {
1647      case(ReducedDegreesOfFreedom):          case(ReducedDegreesOfFreedom):
1648      case(DegreesOfFreedom):          case(DegreesOfFreedom):
1649      case(Nodes):          case(Nodes):
1650      case(ReducedNodes):          case(ReducedNodes):
1651      case(Elements):          case(Elements):
1652      case(ReducedElements):          case(ReducedElements):
1653      case(Points):          case(Points):
1654      case(FaceElements):          case(FaceElements):
1655      case(ReducedFaceElements):          case(ReducedFaceElements):
1656      return true;          return true;
1657      default:          default:
1658          stringstream temp;                  stringstream temp;
1659          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;                  temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1660          throw DudleyAdapterException(temp.str());                  throw DudleyAdapterException(temp.str());
1661      }          }
1662      break;          break;
1663     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1664     switch(functionSpaceType_target) {     switch(functionSpaceType_target) {
1665      case(ReducedDegreesOfFreedom):          case(ReducedDegreesOfFreedom):
1666      case(ReducedNodes):          case(ReducedNodes):
1667      case(Elements):          case(Elements):
1668      case(ReducedElements):          case(ReducedElements):
1669      case(FaceElements):          case(FaceElements):
1670      case(ReducedFaceElements):          case(ReducedFaceElements):
1671      case(Points):          case(Points):
1672      return true;          return true;
1673      case(Nodes):          case(Nodes):
1674      case(DegreesOfFreedom):          case(DegreesOfFreedom):
1675      return false;          return false;
1676      default:          default:
1677          stringstream temp;                  stringstream temp;
1678          temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;                  temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_target;
1679          throw DudleyAdapterException(temp.str());                  throw DudleyAdapterException(temp.str());
1680      }          }
1681      break;          break;
1682     default:     default:
1683        stringstream temp;        stringstream temp;
1684        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;        temp << "Error - Interpolation On Domain: Dudley does not know anything about function space type " << functionSpaceType_source;
# Line 1705  int MeshAdapter::getSystemMatrixTypeId(c Line 1713  int MeshAdapter::getSystemMatrixTypeId(c
1713  {  {
1714      const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);      const escript::SolverBuddy& sb = bp::extract<escript::SolverBuddy>(options);
1715    
1716      return SystemMatrixAdapter::getSystemMatrixTypeId(sb.getSolverMethod(),      return paso::SystemMatrix::getSystemMatrixTypeId(sb.getSolverMethod(),
1717                  sb.getPreconditioner(), sb.getPackage(), sb.isSymmetric(),                  sb.getPreconditioner(), sb.getPackage(), sb.isSymmetric(),
1718                  m_dudleyMesh->MPIInfo);                  m_dudleyMesh->MPIInfo);
1719  }  }
1720    
1721  int MeshAdapter::getTransportTypeId(int solver, int preconditioner, int package, bool symmetry) const  int MeshAdapter::getTransportTypeId(int solver, int preconditioner,
1722                                        int package, bool symmetry) const
1723  {  {
1724     Dudley_Mesh* mesh=m_dudleyMesh.get();     Dudley_Mesh* mesh=m_dudleyMesh.get();
1725     return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,     return paso::TransportProblem::getTypeId(solver, preconditioner, package,
1726             package, symmetry, mesh->MPIInfo);                                              symmetry, mesh->MPIInfo);
1727  }  }
1728    
1729  escript::Data MeshAdapter::getX() const  escript::Data MeshAdapter::getX() const
# Line 1977  bool MeshAdapter::canTag(int functionSpa Line 1986  bool MeshAdapter::canTag(int functionSpa
1986     case(ReducedNodes):     case(ReducedNodes):
1987     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1988     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1989        return false;            return false;
1990     default:     default:
1991      return false;          return false;
1992    }    }
1993  }  }
1994    
# Line 2036  escript::Data MeshAdapter::randomFill(co Line 2045  escript::Data MeshAdapter::randomFill(co
2045      escript::DataTypes::ValueType& dv=towipe.getExpandedVectorReference();      escript::DataTypes::ValueType& dv=towipe.getExpandedVectorReference();
2046      const size_t dvsize=dv.size();      const size_t dvsize=dv.size();
2047      esysUtils::randomFillArray(seed, &(dv[0]), dvsize);      esysUtils::randomFillArray(seed, &(dv[0]), dvsize);
2048      return towipe;            return towipe;      
2049  }  }
2050    
2051    

Legend:
Removed from v.5932  
changed lines
  Added in v.5933

  ViewVC Help
Powered by ViewVC 1.1.26