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

Diff of /trunk/finley/src/CPPAdapter/MeshAdapter.cpp

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

revision 1755 by ksteube, Mon Sep 8 01:34:40 2008 UTC revision 2150 by caltinay, Wed Dec 10 05:57:12 2008 UTC
# Line 1  Line 1 
1    
 /* $Id$ */  
   
2  /*******************************************************  /*******************************************************
3   *  *
4   *           Copyright 2003-2007 by ACceSS MNRF  * Copyright (c) 2003-2008 by University of Queensland
5   *       Copyright 2007 by University of Queensland  * Earth Systems Science Computational Center (ESSCC)
6   *  * http://www.uq.edu.au/esscc
7   *                http://esscc.uq.edu.au  *
8   *        Primary Business: Queensland, Australia  * Primary Business: Queensland, Australia
9   *  Licensed under the Open Software License version 3.0  * Licensed under the Open Software License version 3.0
10   *     http://www.opensource.org/licenses/osl-3.0.php  * http://www.opensource.org/licenses/osl-3.0.php
11   *  *
12   *******************************************************/  *******************************************************/
13    
14    
15  #include "MeshAdapter.h"  #include "MeshAdapter.h"
16  #include "escript/Data.h"  #include "escript/Data.h"
# Line 19  Line 18 
18  #ifdef USE_NETCDF  #ifdef USE_NETCDF
19  #include <netcdfcpp.h>  #include <netcdfcpp.h>
20  #endif  #endif
21    #ifdef PASO_MPI
22    #include <mpi.h>
23    #include "paso/Paso_MPI.h"
24    #endif
25  extern "C" {  extern "C" {
26  #include "escript/blocktimer.h"  #include "esysUtils/blocktimer.h"
27  }  }
28  #include <vector>  #include <vector>
29    
 #define IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH  256  
   
30  using namespace std;  using namespace std;
31  using namespace escript;  using namespace escript;
32    
# Line 83  int MeshAdapter::getMPIRank() const Line 84  int MeshAdapter::getMPIRank() const
84  {  {
85     return m_finleyMesh.get()->MPIInfo->rank;     return m_finleyMesh.get()->MPIInfo->rank;
86  }  }
87    void MeshAdapter::MPIBarrier() const
88    {
89    #ifdef PASO_MPI
90       MPI_Barrier(m_finleyMesh.get()->MPIInfo->comm);
91    #endif
92       return;
93    }
94    bool MeshAdapter::onMasterProcessor() const
95    {
96       return m_finleyMesh.get()->MPIInfo->rank == 0;
97    }
98    
99    
100  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {  Finley_Mesh* MeshAdapter::getFinley_Mesh() const {
# Line 106  void MeshAdapter::Print_Mesh_Info(const Line 118  void MeshAdapter::Print_Mesh_Info(const
118  void MeshAdapter::dump(const std::string& fileName) const  void MeshAdapter::dump(const std::string& fileName) const
119  {  {
120  #ifdef USE_NETCDF  #ifdef USE_NETCDF
121     const NcDim* ncdims[12];     const NcDim* ncdims[12] = {NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL};
122     NcVar *ids;     NcVar *ids;
123     int *int_ptr;     int *int_ptr;
124     Finley_Mesh *mesh = m_finleyMesh.get();     Finley_Mesh *mesh = m_finleyMesh.get();
# Line 123  void MeshAdapter::dump(const std::string Line 135  void MeshAdapter::dump(const std::string
135     int num_Elements_numNodes        = mesh->Elements->numNodes;     int num_Elements_numNodes        = mesh->Elements->numNodes;
136     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;     int num_FaceElements_numNodes    = mesh->FaceElements->numNodes;
137     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;     int num_ContactElements_numNodes = mesh->ContactElements->numNodes;
138    #ifdef PASO_MPI
139       MPI_Status status;
140    #endif
141    
142    /* Incoming token indicates it's my turn to write */
143    #ifdef PASO_MPI
144       if (mpi_rank>0) MPI_Recv(&num_Tags, 0, MPI_INT, mpi_rank-1, 81800, mesh->MPIInfo->comm, &status);
145    #endif
146    
147     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),     char *newFileName = Paso_MPI_appendRankToFileName(fileName.c_str(),
148                                                       mpi_size, mpi_rank);                                                       mpi_size, mpi_rank);
149    
150     /* Figure out how much storage is required for tags */     /* Figure out how much storage is required for tags */
151     tag_map = mesh->TagMap;     tag_map = mesh->TagMap;
152       num_Tags = 0;
153     if (tag_map) {     if (tag_map) {
154        while (tag_map) {        while (tag_map) {
155           num_Tags++;           num_Tags++;
# Line 275  void MeshAdapter::dump(const std::string Line 296  void MeshAdapter::dump(const std::string
296        int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];        int_ptr = &mesh->Nodes->degreesOfFreedomDistribution->first_component[0];
297        if (! (ids->put(int_ptr, mpi_size+1)) )        if (! (ids->put(int_ptr, mpi_size+1)) )
298           throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_DofDistribution to netCDF buffer failed: " + *newFileName);
 { int i; printf("ksteube wrote Nodes_DofDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", mesh->Nodes->degreesOfFreedomDistribution->first_component[i]); printf("\n"); }  
299    
300        // Nodes nodeDistribution        // Nodes nodeDistribution
301        if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )        if (! ( ids = dataFile.add_var("Nodes_NodeDistribution", ncInt, ncdims[2])) )
# Line 283  void MeshAdapter::dump(const std::string Line 303  void MeshAdapter::dump(const std::string
303        int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];        int_ptr = &mesh->Nodes->nodesDistribution->first_component[0];
304        if (! (ids->put(int_ptr, mpi_size+1)) )        if (! (ids->put(int_ptr, mpi_size+1)) )
305           throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Nodes_NodeDistribution to netCDF buffer failed: " + *newFileName);
 { int i; printf("ksteube wrote Nodes_NodeDistribution:"); for(i=0; i<mpi_size+1; i++) printf(" %d", mesh->Nodes->nodesDistribution->first_component[i]); printf("\n"); }  
306    
307     }     }
308    
# Line 291  void MeshAdapter::dump(const std::string Line 310  void MeshAdapter::dump(const std::string
310    
311     if (num_Elements>0) {     if (num_Elements>0) {
312    
       // Temp storage to gather node IDs  
       int *Elements_Nodes = TMPMEMALLOC(num_Elements*num_Elements_numNodes,int);  
   
313        // Elements_Id        // Elements_Id
314        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )        if (! ( ids = dataFile.add_var("Elements_Id", ncInt, ncdims[3])) )
315           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Elements_Id to netCDF file failed: " + *newFileName);
# Line 323  void MeshAdapter::dump(const std::string Line 339  void MeshAdapter::dump(const std::string
339           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Color to netCDF buffer failed: " + *newFileName);
340    
341        // Elements_Nodes        // Elements_Nodes
       for (int i=0; i<num_Elements; i++)  
          for (int j=0; j<num_Elements_numNodes; j++)  
             Elements_Nodes[INDEX2(j,i,num_Elements_numNodes)] = mesh->Nodes->Id[mesh->Elements->Nodes[INDEX2(j,i,num_Elements_numNodes)]];  
342        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )        if (! ( ids = dataFile.add_var("Elements_Nodes", ncInt, ncdims[3], ncdims[7]) ) )
343           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Elements_Nodes to netCDF file failed: " + *newFileName);
344        if (! (ids->put(&(Elements_Nodes[0]), num_Elements, num_Elements_numNodes)) )        if (! (ids->put(&(mesh->Elements->Nodes[0]), num_Elements, num_Elements_numNodes)) )
345           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Elements_Nodes to netCDF buffer failed: " + *newFileName);
346    
       TMPMEMFREE(Elements_Nodes);  
   
347     }     }
348    
349     // // // // // Face_Elements // // // // //     // // // // // Face_Elements // // // // //
350    
351     if (num_FaceElements>0) {     if (num_FaceElements>0) {
352    
       // Temp storage to gather node IDs  
       int *FaceElements_Nodes = TMPMEMALLOC(num_FaceElements*num_FaceElements_numNodes,int);  
   
353        // FaceElements_Id        // FaceElements_Id
354        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )        if (! ( ids = dataFile.add_var("FaceElements_Id", ncInt, ncdims[4])) )
355           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Id to netCDF file failed: " + *newFileName);
# Line 371  void MeshAdapter::dump(const std::string Line 379  void MeshAdapter::dump(const std::string
379           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Color to netCDF buffer failed: " + *newFileName);
380    
381        // FaceElements_Nodes        // FaceElements_Nodes
       for (int i=0; i<num_FaceElements; i++)  
          for (int j=0; j<num_FaceElements_numNodes; j++)  
             FaceElements_Nodes[INDEX2(j,i,num_FaceElements_numNodes)] = mesh->Nodes->Id[mesh->FaceElements->Nodes[INDEX2(j,i,num_FaceElements_numNodes)]];  
382        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )        if (! ( ids = dataFile.add_var("FaceElements_Nodes", ncInt, ncdims[4], ncdims[8]) ) )
383           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending FaceElements_Nodes to netCDF file failed: " + *newFileName);
384        if (! (ids->put(&(FaceElements_Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )        if (! (ids->put(&(mesh->FaceElements->Nodes[0]), num_FaceElements, num_FaceElements_numNodes)) )
385           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy FaceElements_Nodes to netCDF buffer failed: " + *newFileName);
386    
       TMPMEMFREE(FaceElements_Nodes);  
   
387     }     }
388    
389     // // // // // Contact_Elements // // // // //     // // // // // Contact_Elements // // // // //
390    
391     if (num_ContactElements>0) {     if (num_ContactElements>0) {
392    
       // Temp storage to gather node IDs  
       int *ContactElements_Nodes = TMPMEMALLOC(num_ContactElements*num_ContactElements_numNodes,int);  
   
393        // ContactElements_Id        // ContactElements_Id
394        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )        if (! ( ids = dataFile.add_var("ContactElements_Id", ncInt, ncdims[5])) )
395           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Id to netCDF file failed: " + *newFileName);
# Line 419  void MeshAdapter::dump(const std::string Line 419  void MeshAdapter::dump(const std::string
419           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Color to netCDF buffer failed: " + *newFileName);
420    
421        // ContactElements_Nodes        // ContactElements_Nodes
       for (int i=0; i<num_ContactElements; i++)  
          for (int j=0; j<num_ContactElements_numNodes; j++)  
             ContactElements_Nodes[INDEX2(j,i,num_ContactElements_numNodes)] = mesh->Nodes->Id[mesh->ContactElements->Nodes[INDEX2(j,i,num_ContactElements_numNodes)]];  
422        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )        if (! ( ids = dataFile.add_var("ContactElements_Nodes", ncInt, ncdims[5], ncdims[9]) ) )
423           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending ContactElements_Nodes to netCDF file failed: " + *newFileName);
424        if (! (ids->put(&(ContactElements_Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )        if (! (ids->put(&(mesh->ContactElements->Nodes[0]), num_ContactElements, num_ContactElements_numNodes)) )
425           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy ContactElements_Nodes to netCDF buffer failed: " + *newFileName);
426    
       TMPMEMFREE(ContactElements_Nodes);  
   
427     }     }
428    
429     // // // // // Points // // // // //     // // // // // Points // // // // //
# Line 437  void MeshAdapter::dump(const std::string Line 432  void MeshAdapter::dump(const std::string
432    
433        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");        fprintf(stderr, "\n\n\nWARNING: MeshAdapter::dump has not been tested with Point elements\n\n\n");
434    
       // Temp storage to gather node IDs  
       int *Points_Nodes = TMPMEMALLOC(num_Points,int);  
   
435        // Points_Id        // Points_Id
436        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )        if (! ( ids = dataFile.add_var("Points_Id", ncInt, ncdims[6])) )
437           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Id to netCDF file failed: " + *newFileName);
# Line 470  void MeshAdapter::dump(const std::string Line 462  void MeshAdapter::dump(const std::string
462    
463        // Points_Nodes        // Points_Nodes
464        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]        // mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]]
       for (int i=0; i<num_Points; i++)  
          Points_Nodes[i] = mesh->Nodes->Id[mesh->Points->Nodes[INDEX2(0,i,1)]];  
465        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )        if (! ( ids = dataFile.add_var("Points_Nodes", ncInt, ncdims[6]) ) )
466           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: appending Points_Nodes to netCDF file failed: " + *newFileName);
467        if (! (ids->put(&(Points_Nodes[0]), num_Points)) )        if (! (ids->put(&(mesh->Points->Nodes[0]), num_Points)) )
468           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);           throw DataException("Error - MeshAdapter::dump: copy Points_Nodes to netCDF buffer failed: " + *newFileName);
469    
       TMPMEMFREE(Points_Nodes);  
   
470     }     }
471    
472     // // // // // TagMap // // // // //     // // // // // TagMap // // // // //
# Line 525  void MeshAdapter::dump(const std::string Line 513  void MeshAdapter::dump(const std::string
513    
514     }     }
515    
516    /* Send token to next MPI process so he can take his turn */
517    #ifdef PASO_MPI
518       if (mpi_rank<mpi_size-1) MPI_Send(&num_Tags, 0, MPI_INT, mpi_rank+1, 81800, mesh->MPIInfo->comm);
519    #endif
520    
521     // NetCDF file is closed by destructor of NcFile object     // NetCDF file is closed by destructor of NcFile object
522    
523  #else  #else
524     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");     Finley_setError(IO_ERROR, "MeshAdapter::dump: not configured with NetCDF. Please contact your installation manager.");
525  #endif  /* USE_NETCDF */  #endif  /* USE_NETCDF */
# Line 651  int MeshAdapter::getDiracDeltaFunctionCo Line 644  int MeshAdapter::getDiracDeltaFunctionCo
644  //  //
645  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
646  {  {
647     int numDim=Finley_Mesh_getDim(m_finleyMesh.get());     Finley_Mesh* mesh=m_finleyMesh.get();
648       int numDim=Finley_Mesh_getDim(mesh);
649     checkFinleyError();     checkFinleyError();
650     return numDim;     return numDim;
651  }  }
# Line 840  void MeshAdapter::addPDEToTransportProbl Line 834  void MeshAdapter::addPDEToTransportProbl
834                                             const escript::Data& d, const escript::Data& y,                                             const escript::Data& d, const escript::Data& y,
835                                             const escript::Data& d_contact,const escript::Data& y_contact) const                                             const escript::Data& d_contact,const escript::Data& y_contact) const
836  {  {
837     DataArrayView::ShapeType shape;     DataTypes::ShapeType shape;
838     source.expand();     source.expand();
839     escriptDataC _source=source.getDataC();     escriptDataC _source=source.getDataC();
840     escriptDataC _M=M.getDataC();     escriptDataC _M=M.getDataC();
# Line 876  void MeshAdapter::addPDEToTransportProbl Line 870  void MeshAdapter::addPDEToTransportProbl
870  //  //
871  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const  void MeshAdapter::interpolateOnDomain(escript::Data& target,const escript::Data& in) const
872  {  {
873     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());     const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(*(in.getFunctionSpace().getDomain()));
874     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(*(target.getFunctionSpace().getDomain()));
875     if (inDomain!=*this)       if (inDomain!=*this)  
876        throw FinleyAdapterException("Error - Illegal domain of interpolant.");        throw FinleyAdapterException("Error - Illegal domain of interpolant.");
877     if (targetDomain!=*this)     if (targetDomain!=*this)
# Line 1161  void MeshAdapter::interpolateOnDomain(es Line 1155  void MeshAdapter::interpolateOnDomain(es
1155  //  //
1156  void MeshAdapter::setToX(escript::Data& arg) const  void MeshAdapter::setToX(escript::Data& arg) const
1157  {  {
1158     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1159     if (argDomain!=*this)     if (argDomain!=*this)
1160        throw FinleyAdapterException("Error - Illegal domain of data point locations");        throw FinleyAdapterException("Error - Illegal domain of data point locations");
1161     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1184  void MeshAdapter::setToX(escript::Data& Line 1178  void MeshAdapter::setToX(escript::Data&
1178  //  //
1179  void MeshAdapter::setToNormal(escript::Data& normal) const  void MeshAdapter::setToNormal(escript::Data& normal) const
1180  {  {
1181     const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());  /*   const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(normal.getFunctionSpace().getDomain());*/
1182       const MeshAdapter& normalDomain=dynamic_cast<const MeshAdapter&>(*(normal.getFunctionSpace().getDomain()));
1183     if (normalDomain!=*this)     if (normalDomain!=*this)
1184        throw FinleyAdapterException("Error - Illegal domain of normal locations");        throw FinleyAdapterException("Error - Illegal domain of normal locations");
1185     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
# Line 1239  void MeshAdapter::setToNormal(escript::D Line 1234  void MeshAdapter::setToNormal(escript::D
1234  //  //
1235  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const  void MeshAdapter::interpolateACross(escript::Data& target,const escript::Data& source) const
1236  {  {
1237     const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());     const_Domain_ptr targetDomain_p=target.getFunctionSpace().getDomain();
1238     if (targetDomain!=*this)     const MeshAdapter* targetDomain=dynamic_cast<const MeshAdapter*>(targetDomain_p.get());
1239       if (targetDomain!=this)
1240        throw FinleyAdapterException("Error - Illegal domain of interpolation target");        throw FinleyAdapterException("Error - Illegal domain of interpolation target");
1241    
1242     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");     throw FinleyAdapterException("Error - Finley does not allow interpolation across domains yet.");
# Line 1251  void MeshAdapter::interpolateACross(escr Line 1247  void MeshAdapter::interpolateACross(escr
1247  //  //
1248  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const  void MeshAdapter::setToIntegrals(std::vector<double>& integrals,const escript::Data& arg) const
1249  {  {
1250     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1251     if (argDomain!=*this)     if (argDomain!=*this)
1252        throw FinleyAdapterException("Error - Illegal domain of integration kernel");        throw FinleyAdapterException("Error - Illegal domain of integration kernel");
1253    
# Line 1262  void MeshAdapter::setToIntegrals(std::ve Line 1258  void MeshAdapter::setToIntegrals(std::ve
1258     escriptDataC _arg=arg.getDataC();     escriptDataC _arg=arg.getDataC();
1259     switch(arg.getFunctionSpace().getTypeCode()) {     switch(arg.getFunctionSpace().getTypeCode()) {
1260     case(Nodes):     case(Nodes):
1261     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1262     _temp=temp.getDataC();     _temp=temp.getDataC();
1263     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1264     break;     break;
1265     case(ReducedNodes):     case(ReducedNodes):
1266     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1267     _temp=temp.getDataC();     _temp=temp.getDataC();
1268     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1269     break;     break;
# Line 1299  void MeshAdapter::setToIntegrals(std::ve Line 1295  void MeshAdapter::setToIntegrals(std::ve
1295     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->ContactElements,&_arg,&integrals[0]);
1296     break;     break;
1297     case(DegreesOfFreedom):     case(DegreesOfFreedom):
1298     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1299     _temp=temp.getDataC();     _temp=temp.getDataC();
1300     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1301     break;     break;
1302     case(ReducedDegreesOfFreedom):     case(ReducedDegreesOfFreedom):
1303     temp=escript::Data( arg, function(asAbstractContinuousDomain()) );     temp=escript::Data( arg, escript::function(asAbstractContinuousDomain()) );
1304     _temp=temp.getDataC();     _temp=temp.getDataC();
1305     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);     Finley_Assemble_integrate(mesh->Nodes,mesh->Elements,&_temp,&integrals[0]);
1306     break;     break;
# Line 1323  void MeshAdapter::setToIntegrals(std::ve Line 1319  void MeshAdapter::setToIntegrals(std::ve
1319  //  //
1320  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const  void MeshAdapter::setToGradient(escript::Data& grad,const escript::Data& arg) const
1321  {  {
1322     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(arg.getFunctionSpace().getDomain());     const MeshAdapter& argDomain=dynamic_cast<const MeshAdapter&>(*(arg.getFunctionSpace().getDomain()));
1323     if (argDomain!=*this)     if (argDomain!=*this)
1324        throw FinleyAdapterException("Error - Illegal domain of gradient argument");        throw FinleyAdapterException("Error - Illegal domain of gradient argument");
1325     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(grad.getFunctionSpace().getDomain());     const MeshAdapter& gradDomain=dynamic_cast<const MeshAdapter&>(*(grad.getFunctionSpace().getDomain()));
1326     if (gradDomain!=*this)     if (gradDomain!=*this)
1327        throw FinleyAdapterException("Error - Illegal domain of gradient");        throw FinleyAdapterException("Error - Illegal domain of gradient");
1328    
# Line 1453  void MeshAdapter::setNewX(const escript: Line 1449  void MeshAdapter::setNewX(const escript:
1449  {  {
1450     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
1451     escriptDataC tmp;     escriptDataC tmp;
1452     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());     const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(*(new_x.getFunctionSpace().getDomain()));
1453     if (newDomain!=*this)     if (newDomain!=*this)
1454        throw FinleyAdapterException("Error - Illegal domain of new point locations");        throw FinleyAdapterException("Error - Illegal domain of new point locations");
1455     tmp = new_x.getDataC();     tmp = new_x.getDataC();
# Line 1464  void MeshAdapter::setNewX(const escript: Line 1460  void MeshAdapter::setNewX(const escript:
1460  // saves a data array in openDX format:  // saves a data array in openDX format:
1461  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveDX(const std::string& filename,const boost::python::dict& arg) const
1462  {  {
    unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;  
1463     const int num_data=boost::python::extract<int>(arg.attr("__len__")());     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1464     /* win32 refactor */     /* win32 refactor */
1465     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     char **names = (num_data>0) ? TMPMEMALLOC(num_data, char*) : (char**)NULL;
    for(int i=0;i<num_data;i++)  
    {  
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
1466    
1467     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1468     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
1469    
1470     boost::python::list keys=arg.keys();     boost::python::list keys=arg.keys();
1471     for (int i=0;i<num_data;++i) {     for (int i=0; i<num_data; ++i) {
1472        std::string n=boost::python::extract<std::string>(keys[i]);        std::string n=boost::python::extract<std::string>(keys[i]);
1473        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1474        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1475           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error in saveDX: Data must be defined on same Domain");
1476        data[i]=d.getDataC();        data[i]=d.getDataC();
1477        ptr_data[i]= &(data[i]);        ptr_data[i]= &(data[i]);
1478        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1479           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i], n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1480     }     }
1481     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);     Finley_Mesh_saveDX(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);
1482     checkFinleyError();     checkFinleyError();
# Line 1497  void MeshAdapter::saveDX(const std::stri Line 1484  void MeshAdapter::saveDX(const std::stri
1484     /* win32 refactor */     /* win32 refactor */
1485     TMPMEMFREE(data);     TMPMEMFREE(data);
1486     TMPMEMFREE(ptr_data);     TMPMEMFREE(ptr_data);
1487     for(int i=0;i<num_data;i++)     for(int i=0; i<num_data; i++)
1488     {     {
1489        TMPMEMFREE(names[i]);        TMPMEMFREE(names[i]);
1490     }     }
# Line 1509  void MeshAdapter::saveDX(const std::stri Line 1496  void MeshAdapter::saveDX(const std::stri
1496  // saves a data array in openVTK format:  // saves a data array in openVTK format:
1497  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const  void MeshAdapter::saveVTK(const std::string& filename,const boost::python::dict& arg) const
1498  {  {
    unsigned int MAX_namelength=IS_THERE_A_REASON_FOR_THIS_MAGIC_NAME_LENGTH;  
1499     const int num_data=boost::python::extract<int>(arg.attr("__len__")());     const int num_data=boost::python::extract<int>(arg.attr("__len__")());
1500     /* win32 refactor */     /* win32 refactor */
1501     char* *names = (num_data>0) ? TMPMEMALLOC(num_data,char*) : (char**)NULL;     char **names = (num_data>0) ? TMPMEMALLOC(num_data, char*) : (char**)NULL;
    for(int i=0;i<num_data;i++)  
    {  
       names[i] = (MAX_namelength>0) ? TMPMEMALLOC(MAX_namelength,char) : (char*)NULL;  
    }  
1502    
1503     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;     escriptDataC *data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC) : (escriptDataC*)NULL;
1504     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;     escriptDataC* *ptr_data = (num_data>0) ? TMPMEMALLOC(num_data,escriptDataC*) : (escriptDataC**)NULL;
# Line 1525  void MeshAdapter::saveVTK(const std::str Line 1507  void MeshAdapter::saveVTK(const std::str
1507     for (int i=0;i<num_data;++i) {     for (int i=0;i<num_data;++i) {
1508        std::string n=boost::python::extract<std::string>(keys[i]);        std::string n=boost::python::extract<std::string>(keys[i]);
1509        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);        escript::Data& d=boost::python::extract<escript::Data&>(arg[keys[i]]);
1510        if (dynamic_cast<const MeshAdapter&>(d.getFunctionSpace().getDomain()) !=*this)        if (dynamic_cast<const MeshAdapter&>(*(d.getFunctionSpace().getDomain())) !=*this)
1511           throw FinleyAdapterException("Error  in saveVTK: Data must be defined on same Domain");           throw FinleyAdapterException("Error in saveVTK: Data must be defined on same Domain");
1512        data[i]=d.getDataC();        data[i]=d.getDataC();
1513        ptr_data[i]=&(data[i]);        ptr_data[i]=&(data[i]);
1514        if (n.length()>MAX_namelength-1) {        names[i] = TMPMEMALLOC(n.length()+1, char);
1515           strncpy(names[i],n.c_str(),MAX_namelength-1);        strcpy(names[i],n.c_str());
          names[i][MAX_namelength-1]='\0';  
       } else {  
          strcpy(names[i],n.c_str());  
       }  
1516     }     }
1517     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);     Finley_Mesh_saveVTK(filename.c_str(),m_finleyMesh.get(),num_data,names,ptr_data);
1518    
# Line 1563  SystemMatrixAdapter MeshAdapter::newSyst Line 1541  SystemMatrixAdapter MeshAdapter::newSyst
1541     int reduceRowOrder=0;     int reduceRowOrder=0;
1542     int reduceColOrder=0;     int reduceColOrder=0;
1543     // is the domain right?     // is the domain right?
1544     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(row_functionspace.getDomain());     const MeshAdapter& row_domain=dynamic_cast<const MeshAdapter&>(*(row_functionspace.getDomain()));
1545     if (row_domain!=*this)     if (row_domain!=*this)
1546        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.");
1547     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(column_functionspace.getDomain());     const MeshAdapter& col_domain=dynamic_cast<const MeshAdapter&>(*(column_functionspace.getDomain()));
1548     if (col_domain!=*this)     if (col_domain!=*this)
1549        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");        throw FinleyAdapterException("Error - domain of columnn function space does not match the domain of matrix generator.");
1550     // is the function space type right     // is the function space type right
# Line 1611  TransportProblemAdapter MeshAdapter::new Line 1589  TransportProblemAdapter MeshAdapter::new
1589  {  {
1590     int reduceOrder=0;     int reduceOrder=0;
1591     // is the domain right?     // is the domain right?
1592     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(functionspace.getDomain());     const MeshAdapter& domain=dynamic_cast<const MeshAdapter&>(*(functionspace.getDomain()));
1593     if (domain!=*this)     if (domain!=*this)
1594        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");        throw FinleyAdapterException("Error - domain of function space does not match the domain of transport problem generator.");
1595     // is the function space type right     // is the function space type right
# Line 1842  bool MeshAdapter::operator!=(const Abstr Line 1820  bool MeshAdapter::operator!=(const Abstr
1820     return !(operator==(other));     return !(operator==(other));
1821  }  }
1822    
1823  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int package, const bool symmetry) const  int MeshAdapter::getSystemMatrixTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1824  {  {
1825     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);     int out=Paso_SystemMatrix_getSystemMatrixTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1826       checkPasoError();
1827       return out;
1828    }
1829    int MeshAdapter::getTransportTypeId(const int solver, const int preconditioner, const int package, const bool symmetry) const
1830    {
1831       int out=Paso_FCTransportProblem_getTypeId(SystemMatrixAdapter::mapOptionToPaso(solver),SystemMatrixAdapter::mapOptionToPaso(preconditioner), SystemMatrixAdapter::mapOptionToPaso(package),symmetry?1:0);
1832     checkPasoError();     checkPasoError();
1833     return out;     return out;
1834  }  }
# Line 1861  escript::Data MeshAdapter::getNormal() c Line 1845  escript::Data MeshAdapter::getNormal() c
1845    
1846  escript::Data MeshAdapter::getSize() const  escript::Data MeshAdapter::getSize() const
1847  {  {
1848     return function(asAbstractContinuousDomain()).getSize();     return escript::function(asAbstractContinuousDomain()).getSize();
1849  }  }
1850    
1851  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const  int* MeshAdapter::borrowSampleReferenceIDs(int functionSpaceType) const
# Line 2144  int* MeshAdapter::borrowListOfTagsInUse( Line 2128  int* MeshAdapter::borrowListOfTagsInUse(
2128  }  }
2129    
2130    
2131    bool MeshAdapter::canTag(int functionSpaceCode) const
2132    {
2133      switch(functionSpaceCode) {
2134       case(Nodes):
2135       case(Elements):
2136       case(ReducedElements):
2137       case(FaceElements):
2138       case(ReducedFaceElements):
2139       case(Points):
2140       case(ContactElementsZero):
2141       case(ReducedContactElementsZero):
2142       case(ContactElementsOne):
2143       case(ReducedContactElementsOne):
2144              return true;
2145       case(ReducedNodes):
2146       case(DegreesOfFreedom):
2147       case(ReducedDegreesOfFreedom):
2148          return false;
2149       default:
2150        return false;
2151      }
2152    }
2153    
2154    
2155  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.1755  
changed lines
  Added in v.2150

  ViewVC Help
Powered by ViewVC 1.1.26