/[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 4492 by caltinay, Tue Jul 2 01:44:11 2013 UTC revision 4800 by caltinay, Wed Mar 26 01:50:04 2014 UTC
# Line 1  Line 1 
1    
2  /*****************************************************************************  /*****************************************************************************
3  *  *
4  * Copyright (c) 2003-2013 by University of Queensland  * Copyright (c) 2003-2014 by University of Queensland
5  * http://www.uq.edu.au  * http://www.uq.edu.au
6  *  *
7  * Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
# Line 9  Line 9 
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)  * Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12  * Development since 2012 by School of Earth Sciences  * Development 2012-2013 by School of Earth Sciences
13    * Development from 2014 by Centre for Geoscience Computing (GeoComp)
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
# Line 19  Line 20 
20  #include "escript/Data.h"  #include "escript/Data.h"
21  #include "escript/DataFactory.h"  #include "escript/DataFactory.h"
22  #include "esysUtils/blocktimer.h"  #include "esysUtils/blocktimer.h"
23    #include "esysUtils/EsysRandom.h"
24    
25  #include <boost/python/import.hpp>  #include <boost/python/import.hpp>
26  #ifdef USE_NETCDF  #ifdef USE_NETCDF
# Line 46  const int MeshAdapter::ReducedContactEle Line 48  const int MeshAdapter::ReducedContactEle
48  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;  const int MeshAdapter::ContactElementsOne=FINLEY_CONTACT_ELEMENTS_2;
49  const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;  const int MeshAdapter::ReducedContactElementsOne=FINLEY_REDUCED_CONTACT_ELEMENTS_2;
50    
51  MeshAdapter::MeshAdapter(Finley_Mesh* finleyMesh)  MeshAdapter::MeshAdapter(Mesh* finleyMesh)
52  {  {
53      setFunctionSpaceTypeNames();      setFunctionSpaceTypeNames();
54      // need to use a null_deleter as Finley_Mesh_free deletes the pointer      // need to use a null_deleter as Finley_Mesh_free deletes the pointer
# Line 66  MeshAdapter::~MeshAdapter() Line 68  MeshAdapter::~MeshAdapter()
68  {  {
69      // 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
70      if (m_finleyMesh.unique()) {      if (m_finleyMesh.unique()) {
71          Finley_Mesh_free(m_finleyMesh.get());          delete m_finleyMesh.get();
72      }      }
73  }  }
74    
# Line 106  MeshAdapter::getMPIComm() const Line 108  MeshAdapter::getMPIComm() const
108  #endif  #endif
109  }  }
110    
111  Finley_Mesh* MeshAdapter::getFinley_Mesh() const  Mesh* MeshAdapter::getFinley_Mesh() const
112  {  {
113      return m_finleyMesh.get();      return m_finleyMesh.get();
114  }  }
115    
116  void MeshAdapter::write(const string& fileName) const  void MeshAdapter::write(const string& fileName) const
117  {  {
118      Finley_Mesh_write(m_finleyMesh.get(), fileName.c_str());      m_finleyMesh.get()->write(fileName);
119      checkFinleyError();      checkFinleyError();
120  }  }
121    
122  void MeshAdapter::Print_Mesh_Info(const bool full) const  void MeshAdapter::Print_Mesh_Info(const bool full) const
123  {  {
124      Finley_PrintMesh_Info(m_finleyMesh.get(), full);      m_finleyMesh.get()->printInfo(full);
125  }  }
126    
127  void MeshAdapter::dump(const string& fileName) const  void MeshAdapter::dump(const string& fileName) const
# Line 128  void MeshAdapter::dump(const string& fil Line 130  void MeshAdapter::dump(const string& fil
130      const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};      const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
131      NcVar *ids;      NcVar *ids;
132      int *int_ptr;      int *int_ptr;
133      Finley_Mesh *mesh = m_finleyMesh.get();      Mesh *mesh = m_finleyMesh.get();
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;
# Line 201  void MeshAdapter::dump(const string& fil Line 203  void MeshAdapter::dump(const string& fil
203              throw FinleyAdapterException(msgPrefix+"add_dim(dim_Tags)");              throw FinleyAdapterException(msgPrefix+"add_dim(dim_Tags)");
204    
205      // Attributes: MPI size, MPI rank, Name, order, reduced_order      // Attributes: MPI size, MPI rank, Name, order, reduced_order
206      if (!dataFile.add_att("mpi_size", mpi_size) )      if (!dataFile.add_att("mpi_size", mpi_size))
207          throw FinleyAdapterException(msgPrefix+"add_att(mpi_size)");          throw FinleyAdapterException(msgPrefix+"add_att(mpi_size)");
208      if (!dataFile.add_att("mpi_rank", mpi_rank) )      if (!dataFile.add_att("mpi_rank", mpi_rank))
209          throw FinleyAdapterException(msgPrefix+"add_att(mpi_rank)");          throw FinleyAdapterException(msgPrefix+"add_att(mpi_rank)");
210      if (!dataFile.add_att("Name",mesh->Name) )      if (!dataFile.add_att("Name",mesh->m_name.c_str()))
211          throw FinleyAdapterException(msgPrefix+"add_att(Name)");          throw FinleyAdapterException(msgPrefix+"add_att(Name)");
212      if (!dataFile.add_att("numDim",numDim) )      if (!dataFile.add_att("numDim",numDim))
213          throw FinleyAdapterException(msgPrefix+"add_att(order)");          throw FinleyAdapterException(msgPrefix+"add_att(order)");
214      if (!dataFile.add_att("order",mesh->integrationOrder) )      if (!dataFile.add_att("order",mesh->integrationOrder))
215          throw FinleyAdapterException(msgPrefix+"add_att(order)");          throw FinleyAdapterException(msgPrefix+"add_att(order)");
216      if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder) )      if (!dataFile.add_att("reduced_order",mesh->reducedIntegrationOrder))
217          throw FinleyAdapterException(msgPrefix+"add_att(reduced_order)");          throw FinleyAdapterException(msgPrefix+"add_att(reduced_order)");
218      if (!dataFile.add_att("numNodes",numNodes) )      if (!dataFile.add_att("numNodes",numNodes))
219          throw FinleyAdapterException(msgPrefix+"add_att(numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(numNodes)");
220      if (!dataFile.add_att("num_Elements",num_Elements) )      if (!dataFile.add_att("num_Elements",num_Elements))
221          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements)");
222      if (!dataFile.add_att("num_FaceElements",num_FaceElements) )      if (!dataFile.add_att("num_FaceElements",num_FaceElements))
223          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements)");
224      if (!dataFile.add_att("num_ContactElements",num_ContactElements) )      if (!dataFile.add_att("num_ContactElements",num_ContactElements))
225          throw FinleyAdapterException(msgPrefix+"add_att(num_ContactElements)");          throw FinleyAdapterException(msgPrefix+"add_att(num_ContactElements)");
226      if (!dataFile.add_att("num_Points",num_Points) )      if (!dataFile.add_att("num_Points",num_Points))
227          throw FinleyAdapterException(msgPrefix+"add_att(num_Points)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Points)");
228      if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes) )      if (!dataFile.add_att("num_Elements_numNodes",num_Elements_numNodes))
229          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements_numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(num_Elements_numNodes)");
230      if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )      if (!dataFile.add_att("num_FaceElements_numNodes",num_FaceElements_numNodes) )
231          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements_numNodes)");          throw FinleyAdapterException(msgPrefix+"add_att(num_FaceElements_numNodes)");
# Line 500  void MeshAdapter::dump(const string& fil Line 502  void MeshAdapter::dump(const string& fil
502     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
503    
504  #else // USE_NETCDF  #else // USE_NETCDF
505      Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with netCDF. Please contact your installation manager.");      setError(IO_ERROR, "MeshAdapter::dump: not configured with netCDF. Please contact your installation manager.");
506  #endif // USE_NETCDF  #endif // USE_NETCDF
507      checkFinleyError();      checkFinleyError();
508  }  }
# Line 628  int MeshAdapter::getDiracDeltaFunctionsC Line 630  int MeshAdapter::getDiracDeltaFunctionsC
630  //  //
631  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
632  {  {
633      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
634      int numDim=Finley_Mesh_getDim(mesh);      return mesh->getDim();
     checkFinleyError();  
     return numDim;  
635  }  }
636    
637  //  //
# Line 650  pair<int,int> MeshAdapter::getDataShape( Line 650  pair<int,int> MeshAdapter::getDataShape(
650  {  {
651      int numDataPointsPerSample=0;      int numDataPointsPerSample=0;
652      int numSamples=0;      int numSamples=0;
653      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
654      switch (functionSpaceCode) {      switch (functionSpaceCode) {
655          case Nodes:          case Nodes:
656              numDataPointsPerSample=1;              numDataPointsPerSample=1;
# Line 751  void MeshAdapter::addPDEToSystem( Line 751  void MeshAdapter::addPDEToSystem(
751      if (!smat)      if (!smat)
752          throw FinleyAdapterException("finley only supports Paso system matrices.");          throw FinleyAdapterException("finley only supports Paso system matrices.");
753    
754      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
755      Paso_SystemMatrix* S = smat->getPaso_SystemMatrix();      Paso_SystemMatrix* S = smat->getPaso_SystemMatrix();
756      Assemble_PDE(mesh->Nodes, mesh->Elements, S, rhs, A, B, C, D, X, Y);      Assemble_PDE(mesh->Nodes, mesh->Elements, S, rhs, A, B, C, D, X, Y);
757      checkFinleyError();      checkFinleyError();
758    
759    
760      Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs,      Assemble_PDE(mesh->Nodes, mesh->FaceElements, S, rhs,
761              escript::Data(), escript::Data(), escript::Data(), d,              escript::Data(), escript::Data(), escript::Data(), d,
762              escript::Data(), y);              escript::Data(), y);
# Line 775  void MeshAdapter::addPDEToLumpedSystem(e Line 776  void MeshAdapter::addPDEToLumpedSystem(e
776          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
777          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
778  {  {
779      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
780      Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);      Assemble_LumpedSystem(mesh->Nodes, mesh->Elements, mat, D, useHRZ);
781      checkFinleyError();      checkFinleyError();
782    
# Line 793  void MeshAdapter::addPDEToRHS(escript::D Line 794  void MeshAdapter::addPDEToRHS(escript::D
794          const escript::Data& Y, const escript::Data& y,          const escript::Data& Y, const escript::Data& y,
795          const escript::Data& y_contact, const escript::Data& y_dirac) const          const escript::Data& y_contact, const escript::Data& y_dirac) const
796  {  {
797      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
798      Assemble_PDE(mesh->Nodes, mesh->Elements, 0, rhs,      Assemble_PDE(mesh->Nodes, mesh->Elements, 0, rhs,
799              escript::Data(), escript::Data(), escript::Data(), escript::Data(),              escript::Data(), escript::Data(), escript::Data(), escript::Data(),
800              X, Y);              X, Y);
# Line 832  void MeshAdapter::addPDEToTransportProbl Line 833  void MeshAdapter::addPDEToTransportProbl
833    
834      source.expand();      source.expand();
835    
836      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
837      Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();      Paso_TransportProblem* _tp = tpa->getPaso_TransportProblem();
838    
839      Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix, source,      Assemble_PDE(mesh->Nodes, mesh->Elements, _tp->mass_matrix, source,
# Line 873  void MeshAdapter::interpolateOnDomain(es Line 874  void MeshAdapter::interpolateOnDomain(es
874      if (targetDomain!=*this)      if (targetDomain!=*this)
875          throw FinleyAdapterException("Error - Illegal domain of interpolation target.");          throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
876    
877      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
878      switch(in.getFunctionSpace().getTypeCode()) {      switch(in.getFunctionSpace().getTypeCode()) {
879          case Nodes:          case Nodes:
880              switch(target.getFunctionSpace().getTypeCode()) {              switch(target.getFunctionSpace().getTypeCode()) {
# Line 1135  void MeshAdapter::setToX(escript::Data& Line 1136  void MeshAdapter::setToX(escript::Data&
1136      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));      const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1137      if (argDomain!=*this)      if (argDomain!=*this)
1138          throw FinleyAdapterException("Error - Illegal domain of data point locations");          throw FinleyAdapterException("Error - Illegal domain of data point locations");
1139      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1140      // in case of appropriate function space we can do the job directly:      // in case of appropriate function space we can do the job directly:
1141      if (arg.getFunctionSpace().getTypeCode()==Nodes) {      if (arg.getFunctionSpace().getTypeCode()==Nodes) {
1142          Assemble_NodeCoordinates(mesh->Nodes, arg);          Assemble_NodeCoordinates(mesh->Nodes, arg);
# Line 1156  void MeshAdapter::setToNormal(escript::D Line 1157  void MeshAdapter::setToNormal(escript::D
1157      const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));      const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1158      if (normalDomain!=*this)      if (normalDomain!=*this)
1159          throw FinleyAdapterException("Error - Illegal domain of normal locations");          throw FinleyAdapterException("Error - Illegal domain of normal locations");
1160      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1161      switch(normal.getFunctionSpace().getTypeCode()) {      switch(normal.getFunctionSpace().getTypeCode()) {
1162          case Nodes:          case Nodes:
1163              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for nodes");
# Line 1172  void MeshAdapter::setToNormal(escript::D Line 1173  void MeshAdapter::setToNormal(escript::D
1173              break;              break;
1174          case FaceElements:          case FaceElements:
1175          case ReducedFaceElements:          case ReducedFaceElements:
1176              Assemble_setNormal(mesh->Nodes, mesh->FaceElements, normal);              Assemble_getNormal(mesh->Nodes, mesh->FaceElements, normal);
1177              break;              break;
1178          case Points:          case Points:
1179              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for point elements");
# Line 1181  void MeshAdapter::setToNormal(escript::D Line 1182  void MeshAdapter::setToNormal(escript::D
1182          case ContactElementsZero:          case ContactElementsZero:
1183          case ReducedContactElementsOne:          case ReducedContactElementsOne:
1184          case ReducedContactElementsZero:          case ReducedContactElementsZero:
1185              Assemble_setNormal(mesh->Nodes, mesh->ContactElements, normal);              Assemble_getNormal(mesh->Nodes, mesh->ContactElements, normal);
1186              break;              break;
1187          case DegreesOfFreedom:          case DegreesOfFreedom:
1188              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");              throw FinleyAdapterException("Error - Finley does not support surface normal vectors for degrees of freedom.");
# Line 1216  void MeshAdapter::setToIntegrals(vector< Line 1217  void MeshAdapter::setToIntegrals(vector<
1217          throw FinleyAdapterException("Error - Illegal domain of integration kernel");          throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1218    
1219      double blocktimer_start = blocktimer_time();      double blocktimer_start = blocktimer_time();
1220      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1221      switch(arg.getFunctionSpace().getTypeCode()) {      switch(arg.getFunctionSpace().getTypeCode()) {
1222          case Nodes:          case Nodes:
1223              {              {
# Line 1276  void MeshAdapter::setToGradient(escript: Line 1277  void MeshAdapter::setToGradient(escript:
1277      if (gradDomain!=*this)      if (gradDomain!=*this)
1278          throw FinleyAdapterException("Error - Illegal domain of gradient");          throw FinleyAdapterException("Error - Illegal domain of gradient");
1279    
1280      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1281      escript::Data nodeData;      escript::Data nodeData;
1282      if (getMPISize()>1) {      if (getMPISize()>1) {
1283          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {          if (arg.getFunctionSpace().getTypeCode() == DegreesOfFreedom) {
# Line 1333  void MeshAdapter::setToGradient(escript: Line 1334  void MeshAdapter::setToGradient(escript:
1334  //  //
1335  void MeshAdapter::setToSize(escript::Data& size) const  void MeshAdapter::setToSize(escript::Data& size) const
1336  {  {
1337      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1338      switch(size.getFunctionSpace().getTypeCode()) {      switch(size.getFunctionSpace().getTypeCode()) {
1339          case Nodes:          case Nodes:
1340              throw FinleyAdapterException("Error - Size of nodes is not supported.");              throw FinleyAdapterException("Error - Size of nodes is not supported.");
# Line 1378  void MeshAdapter::setToSize(escript::Dat Line 1379  void MeshAdapter::setToSize(escript::Dat
1379  //  //
1380  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
1381  {  {
1382      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1383      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));      const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1384      if (newDomain!=*this)      if (newDomain!=*this)
1385          throw FinleyAdapterException("Error - Illegal domain of new point locations");          throw FinleyAdapterException("Error - Illegal domain of new point locations");
1386      if (new_x.getFunctionSpace() == continuousFunction(*this)) {      if (new_x.getFunctionSpace() == continuousFunction(*this)) {
1387          Finley_Mesh_setCoordinates(mesh, new_x);          mesh->setCoordinates(new_x);
1388      } else {      } else {
1389          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");          throw FinleyAdapterException("As of escript version 3.3 SetX() only accepts ContinuousFunction arguments. Please interpolate.");
1390      }      }
# Line 1396  bool MeshAdapter::ownSample(int fs_code, Line 1397  bool MeshAdapter::ownSample(int fs_code,
1397  #ifdef ESYS_MPI  #ifdef ESYS_MPI
1398          int myFirstNode=0, myLastNode=0, k=0;          int myFirstNode=0, myLastNode=0, k=0;
1399          int* globalNodeIndex=0;          int* globalNodeIndex=0;
1400          Finley_Mesh* mesh_p=m_finleyMesh.get();          Mesh* mesh_p=m_finleyMesh.get();
1401          /*          /*
1402           * this method is only used by saveDataCSV which would use the returned           * this method is only used by saveDataCSV which would use the returned
1403           * values for reduced nodes wrongly so this case is disabled for now           * values for reduced nodes wrongly so this case is disabled for now
1404          if (fs_code == FINLEY_REDUCED_NODES) {          if (fs_code == FINLEY_REDUCED_NODES) {
1405              myFirstNode = Finley_NodeFile_getFirstReducedNode(mesh_p->Nodes);              myFirstNode = NodeFile_getFirstReducedNode(mesh_p->Nodes);
1406              myLastNode = Finley_NodeFile_getLastReducedNode(mesh_p->Nodes);              myLastNode = NodeFile_getLastReducedNode(mesh_p->Nodes);
1407              globalNodeIndex = Finley_NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);              globalNodeIndex = NodeFile_borrowGlobalReducedNodesIndex(mesh_p->Nodes);
1408          } else          } else
1409          */          */
1410          if (fs_code == FINLEY_NODES) {          if (fs_code == FINLEY_NODES) {
# Line 1454  escript::ASM_ptr MeshAdapter::newSystemM Line 1455  escript::ASM_ptr MeshAdapter::newSystemM
1455      }      }
1456    
1457      // generate matrix:      // generate matrix:
1458      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(      paso::SystemMatrixPattern* fsystemMatrixPattern=
1459              getFinley_Mesh(), reduceRowOrder, reduceColOrder);          getFinley_Mesh()->getPattern(reduceRowOrder, reduceColOrder);
1460      checkFinleyError();      checkFinleyError();
1461      Paso_SystemMatrix* fsystemMatrix;      Paso_SystemMatrix* fsystemMatrix;
1462      const int trilinos = 0;      const int trilinos = 0;
# Line 1468  escript::ASM_ptr MeshAdapter::newSystemM Line 1469  escript::ASM_ptr MeshAdapter::newSystemM
1469                  row_blocksize, column_blocksize, FALSE);                  row_blocksize, column_blocksize, FALSE);
1470      }      }
1471      checkPasoError();      checkPasoError();
1472      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      paso::SystemMatrixPattern_free(fsystemMatrixPattern);
1473      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);      SystemMatrixAdapter* sma=new SystemMatrixAdapter(fsystemMatrix, row_blocksize, row_functionspace, column_blocksize, column_functionspace);
1474      return escript::ASM_ptr(sma);      return escript::ASM_ptr(sma);
1475  }  }
# Line 1493  escript::ATP_ptr MeshAdapter::newTranspo Line 1494  escript::ATP_ptr MeshAdapter::newTranspo
1494      }      }
1495    
1496      // generate transport problem:      // generate transport problem:
1497      Paso_SystemMatrixPattern* fsystemMatrixPattern=Finley_getPattern(      paso::SystemMatrixPattern* fsystemMatrixPattern=
1498              getFinley_Mesh(), reduceOrder, reduceOrder);          getFinley_Mesh()->getPattern(reduceOrder, reduceOrder);
1499      checkFinleyError();      checkFinleyError();
1500      Paso_TransportProblem* transportProblem;      Paso_TransportProblem* transportProblem;
1501      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);      transportProblem=Paso_TransportProblem_alloc(fsystemMatrixPattern, blocksize);
1502      checkPasoError();      checkPasoError();
1503      Paso_SystemMatrixPattern_free(fsystemMatrixPattern);      paso::SystemMatrixPattern_free(fsystemMatrixPattern);
1504      TransportProblemAdapter* tpa=new TransportProblemAdapter(      TransportProblemAdapter* tpa=new TransportProblemAdapter(
1505              transportProblem, blocksize, functionspace);              transportProblem, blocksize, functionspace);
1506      return escript::ATP_ptr(tpa);      return escript::ATP_ptr(tpa);
# Line 1845  bool MeshAdapter::operator!=(const escri Line 1846  bool MeshAdapter::operator!=(const escri
1846    
1847  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
1848  {  {
1849      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1850      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,      return SystemMatrixAdapter::getSystemMatrixTypeId(solver, preconditioner,
1851                 package, symmetry, mesh->MPIInfo);                 package, symmetry, mesh->MPIInfo);
1852  }  }
1853    
1854  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
1855  {  {
1856      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1857      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,      return TransportProblemAdapter::getTransportTypeId(solver, preconditioner,
1858                  package, symmetry, mesh->MPIInfo);                  package, symmetry, mesh->MPIInfo);
1859  }  }
# Line 1875  escript::Data MeshAdapter::getSize() con Line 1876  escript::Data MeshAdapter::getSize() con
1876  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  const int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
1877  {  {
1878      int *out = NULL;      int *out = NULL;
1879      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1880      switch (functionSpaceType) {      switch (functionSpaceType) {
1881          case Nodes:          case Nodes:
1882              out=mesh->Nodes->Id;              out=mesh->Nodes->Id;
# Line 1917  const int* MeshAdapter::borrowSampleRefe Line 1918  const int* MeshAdapter::borrowSampleRefe
1918  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
1919  {  {
1920      int out=0;      int out=0;
1921      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1922      switch (functionSpaceType) {      switch (functionSpaceType) {
1923          case Nodes:          case Nodes:
1924              out=mesh->Nodes->Tag[sampleNo];              out=mesh->Nodes->Tag[sampleNo];
# Line 1960  int MeshAdapter::getTagFromSampleNo(int Line 1961  int MeshAdapter::getTagFromSampleNo(int
1961    
1962  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const  void MeshAdapter::setTags(const int functionSpaceType, const int newTag, const escript::Data& mask) const
1963  {  {
1964      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
1965      switch(functionSpaceType) {      switch(functionSpaceType) {
1966          case Nodes:          case Nodes:
1967              mesh->Nodes->setTags(newTag, mask);              mesh->Nodes->setTags(newTag, mask);
# Line 1998  void MeshAdapter::setTags(const int func Line 1999  void MeshAdapter::setTags(const int func
1999    
2000  void MeshAdapter::setTagMap(const string& name, int tag)  void MeshAdapter::setTagMap(const string& name, int tag)
2001  {  {
2002      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2003      Finley_Mesh_addTagMap(mesh, name.c_str(), tag);      mesh->addTagMap(name.c_str(), tag);
2004      checkFinleyError();      checkFinleyError();
2005  }  }
2006    
2007  int MeshAdapter::getTag(const string& name) const  int MeshAdapter::getTag(const string& name) const
2008  {  {
2009      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2010      int tag = Finley_Mesh_getTag(mesh, name.c_str());      int tag = mesh->getTag(name.c_str());
2011      checkFinleyError();      checkFinleyError();
2012      return tag;      return tag;
2013  }  }
2014    
2015  bool MeshAdapter::isValidTagName(const string& name) const  bool MeshAdapter::isValidTagName(const string& name) const
2016  {  {
2017      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2018      return Finley_Mesh_isValidTagName(mesh, name.c_str());      return mesh->isValidTagName(name.c_str());
2019  }  }
2020    
2021  string MeshAdapter::showTagNames() const  string MeshAdapter::showTagNames() const
2022  {  {
2023      stringstream temp;      stringstream temp;
2024      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2025      TagMap::const_iterator it = mesh->tagMap.begin();      TagMap::const_iterator it = mesh->tagMap.begin();
2026      while (it != mesh->tagMap.end()) {      while (it != mesh->tagMap.end()) {
2027          temp << it->first;          temp << it->first;
# Line 2033  string MeshAdapter::showTagNames() const Line 2034  string MeshAdapter::showTagNames() const
2034    
2035  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const  int MeshAdapter::getNumberOfTagsInUse(int functionSpaceCode) const
2036  {  {
2037      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2038      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2039          case Nodes:          case Nodes:
2040              return mesh->Nodes->tagsInUse.size();              return mesh->Nodes->tagsInUse.size();
# Line 2067  int MeshAdapter::getNumberOfTagsInUse(in Line 2068  int MeshAdapter::getNumberOfTagsInUse(in
2068    
2069  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const  const int* MeshAdapter::borrowListOfTagsInUse(int functionSpaceCode) const
2070  {  {
2071      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2072      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2073          case Nodes:          case Nodes:
2074              return &mesh->Nodes->tagsInUse[0];              return &mesh->Nodes->tagsInUse[0];
# Line 2120  bool MeshAdapter::canTag(int functionSpa Line 2121  bool MeshAdapter::canTag(int functionSpa
2121    
2122  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const  escript::AbstractDomain::StatusType MeshAdapter::getStatus() const
2123  {  {
2124      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2125      return Finley_Mesh_getStatus(mesh);      return mesh->getStatus();
2126  }  }
2127    
2128  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const  int MeshAdapter::getApproximationOrder(const int functionSpaceCode) const
2129  {  {
2130      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2131      int order =-1;      int order =-1;
2132      switch(functionSpaceCode) {      switch(functionSpaceCode) {
2133          case Nodes:          case Nodes:
# Line 2163  bool MeshAdapter::supportsContactElement Line 2164  bool MeshAdapter::supportsContactElement
2164      return true;      return true;
2165  }  }
2166    
2167  ReferenceElementSetWrapper::ReferenceElementSetWrapper(ElementTypeId id,  escript::Data MeshAdapter::randomFill(const escript::DataTypes::ShapeType& shape,
2168          index_t order, index_t reducedOrder)         const escript::FunctionSpace& what, long seed, const boost::python::tuple& filter) const
2169  {  {
2170      m_refSet = ReferenceElementSet_alloc(id, order, reducedOrder);      escript::Data towipe(0, shape, what, true);
2171        // since we just made this object, no sharing is possible and we don't need to check for
2172        // exlusive write
2173        escript::DataTypes::ValueType& dv=towipe.getExpandedVectorReference();
2174        const size_t dvsize=dv.size();
2175        esysUtils::randomFillArray(seed, &(dv[0]), dvsize);
2176        return towipe;  
2177  }  }
2178    
 ReferenceElementSetWrapper::~ReferenceElementSetWrapper()  
 {  
     ReferenceElementSet_dealloc(m_refSet);  
 }  
2179    
2180  void MeshAdapter::addDiracPoints(const vector<double>& points,  void MeshAdapter::addDiracPoints(const vector<double>& points,
2181                                   const vector<int>& tags) const                                   const vector<int>& tags) const
# Line 2181  void MeshAdapter::addDiracPoints(const v Line 2184  void MeshAdapter::addDiracPoints(const v
2184      const int dim = getDim();      const int dim = getDim();
2185      int numPoints=points.size()/dim;      int numPoints=points.size()/dim;
2186      int numTags=tags.size();      int numTags=tags.size();
2187      Finley_Mesh* mesh=m_finleyMesh.get();      Mesh* mesh=m_finleyMesh.get();
2188    
2189      if ( points.size() % dim != 0 ) {      if ( points.size() % dim != 0 ) {
2190          throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");          throw FinleyAdapterException("Error - number of coords does not appear to be a multiple of dimension.");
# Line 2190  void MeshAdapter::addDiracPoints(const v Line 2193  void MeshAdapter::addDiracPoints(const v
2193      if ((numTags > 0) && (numPoints != numTags))      if ((numTags > 0) && (numPoints != numTags))
2194          throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");          throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");
2195    
2196      Finley_Mesh_addPoints(mesh, numPoints, &points[0], &tags[0]);      mesh->addPoints(numPoints, &points[0], &tags[0]);
2197      checkFinleyError();      checkFinleyError();
2198  }  }
2199    
# Line 2199  void MeshAdapter::addDiracPoints(const v Line 2202  void MeshAdapter::addDiracPoints(const v
2202  //       const int dim = getDim();  //       const int dim = getDim();
2203  //       int numPoints=boost::python::extract<int>(points.attr("__len__")());  //       int numPoints=boost::python::extract<int>(points.attr("__len__")());
2204  //       int numTags=boost::python::extract<int>(tags.attr("__len__")());  //       int numTags=boost::python::extract<int>(tags.attr("__len__")());
2205  //       Finley_Mesh* mesh=m_finleyMesh.get();  //       Mesh* mesh=m_finleyMesh.get();
2206  //  //
2207  //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )  //       if  ( (numTags > 0) && ( numPoints !=  numTags ) )
2208  //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");  //       throw FinleyAdapterException("Error - if tags are given number of tags and points must match.");

Legend:
Removed from v.4492  
changed lines
  Added in v.4800

  ViewVC Help
Powered by ViewVC 1.1.26