/[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 532 by gross, Wed Feb 15 09:45:53 2006 UTC revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC
# Line 15  Line 15 
15    
16  #include "MeshAdapter.h"  #include "MeshAdapter.h"
17    
18  #include "Data.h"  #include "escript/Data.h"
19  #include "DataFactory.h"  #include "escript/DataFactory.h"
20    
21  using namespace std;  using namespace std;
22  using namespace escript;  using namespace escript;
# Line 76  void MeshAdapter::write(const std::strin Line 76  void MeshAdapter::write(const std::strin
76    checkFinleyError();    checkFinleyError();
77  }  }
78    
 // void MeshAdapter::getTagList(int functionSpaceType,  
 //                  int* numTags) const  
 // {  
 //   Finley_Mesh_tagList(m_finleyMesh.get(),functionSpaceType,tagList,numTags);  
 //   return;  
 // }  
   
79  string MeshAdapter::getDescription() const  string MeshAdapter::getDescription() const
80  {  {
81    return "FinleyMesh";    return "FinleyMesh";
# Line 167  int MeshAdapter::getDiracDeltaFunctionCo Line 160  int MeshAdapter::getDiracDeltaFunctionCo
160  }  }
161    
162  //  //
 // returns a pointer to the tag list of samples of functionSpaceType  
 //  
 void MeshAdapter::getTagList(int functionSpaceType, int** tagList,  
                  int* numTags) const  
 {  
   *tagList=NULL;  
   *numTags=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *tagList=mesh->Nodes->Tag;  
       *numTags=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *tagList=mesh->Elements->Tag;  
       *numTags=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *tagList=mesh->FaceElements->Tag;  
       *numTags=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *tagList=mesh->Points->Tag;  
       *numTags=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *tagList=mesh->ContactElements->Tag;  
       *numTags=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *tagList=NULL;  
       *numTags=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*tagList==NULL) {  
     stringstream temp;  
     temp << "Error - no tags available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
 // returns a pointer to the reference no list of samples of functionSpaceType  
 //  
 void MeshAdapter::getReferenceNoList(int functionSpaceType, int** referenceNoList,  
                  int* numReferenceNo) const  
 {  
   *referenceNoList=NULL;  
   *numReferenceNo=0;  
   Finley_Mesh* mesh=m_finleyMesh.get();  
   switch (functionSpaceType) {  
   case(Nodes):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=mesh->Nodes->Id;  
       *numReferenceNo=mesh->Nodes->numNodes;  
     }  
     break;  
   case(Elements):  
     if (mesh->Elements!=NULL) {  
       *referenceNoList=mesh->Elements->Id;  
       *numReferenceNo=mesh->Elements->numElements;  
     }  
     break;  
   case(FaceElements):  
     if (mesh->FaceElements!=NULL) {  
       *referenceNoList=mesh->FaceElements->Id;  
       *numReferenceNo=mesh->FaceElements->numElements;  
     }  
     break;  
   case(Points):  
     if (mesh->Points!=NULL) {  
       *referenceNoList=mesh->Points->Id;  
       *numReferenceNo=mesh->Points->numElements;  
     }  
     break;  
   case(ContactElementsZero):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(ContactElementsOne):  
     if (mesh->ContactElements!=NULL) {  
       *referenceNoList=mesh->ContactElements->Id;  
       *numReferenceNo=mesh->ContactElements->numElements;  
     }  
     break;  
   case(DegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   case(ReducedDegreesOfFreedom):  
     if (mesh->Nodes!=NULL) {  
       *referenceNoList=NULL;  
       *numReferenceNo=0;  
     }  
     break;  
   default:  
     stringstream temp;  
     temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
     break;  
   }  
   if (*referenceNoList==NULL) {  
     stringstream temp;  
     temp << "Error - reference number list available for " << functionSpaceType << " for domain: " << getDescription();  
     throw FinleyAdapterException(temp.str());  
   }  
   return;  
 }  
   
 //  
163  // return the spatial dimension of the Mesh:  // return the spatial dimension of the Mesh:
164  //  //
165  int MeshAdapter::getDim() const  int MeshAdapter::getDim() const
# Line 367  pair<int,int> MeshAdapter::getDataShape( Line 216  pair<int,int> MeshAdapter::getDataShape(
216        case(DegreesOfFreedom):        case(DegreesOfFreedom):
217             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
218               numDataPointsPerSample=1;               numDataPointsPerSample=1;
219    #ifndef PASO_MPI
220               numSamples=mesh->Nodes->numDegreesOfFreedom;               numSamples=mesh->Nodes->numDegreesOfFreedom;
221    #else
222                 numSamples=mesh->Nodes->degreeOfFreedomDistribution->numLocal;
223    #endif
224             }             }
225             break;             break;
226        case(ReducedDegreesOfFreedom):        case(ReducedDegreesOfFreedom):
227             if (mesh->Nodes!=NULL) {             if (mesh->Nodes!=NULL) {
228               numDataPointsPerSample=1;               numDataPointsPerSample=1;
229    #ifndef PASO_MPI
230               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;               numSamples=mesh->Nodes->reducedNumDegreesOfFreedom;
231    #else
232                 numSamples=mesh->Nodes->reducedDegreeOfFreedomDistribution->numLocal;
233    #endif
234             }             }
235             break;             break;
236        default:        default:
# Line 394  void MeshAdapter::addPDEToSystem( Line 251  void MeshAdapter::addPDEToSystem(
251                       const escript::Data& d, const escript::Data& y,                       const escript::Data& d, const escript::Data& y,
252                       const escript::Data& d_contact,const escript::Data& y_contact) const                       const escript::Data& d_contact,const escript::Data& y_contact) const
253  {  {
254      escriptDataC _rhs=rhs.getDataC();
255      escriptDataC _A  =A.getDataC();
256      escriptDataC _B=B.getDataC();
257      escriptDataC _C=C.getDataC();
258      escriptDataC _D=D.getDataC();
259      escriptDataC _X=X.getDataC();
260      escriptDataC _Y=Y.getDataC();
261      escriptDataC _d=d.getDataC();
262      escriptDataC _y=y.getDataC();
263      escriptDataC _d_contact=d_contact.getDataC();
264      escriptDataC _y_contact=y_contact.getDataC();
265    
266     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
267     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(),&(rhs.getDataC()),     Finley_Assemble_PDE(mesh->Nodes,mesh->Elements,mat.getPaso_SystemMatrix(), &_rhs, &_A, &_B, &_C, &_D, &_X, &_Y );
268                         &(A.getDataC()),&(B.getDataC()),&(C.getDataC()),&(D.getDataC()),&(X.getDataC()),&(Y.getDataC()));  
269     checkFinleyError();     checkFinleyError();
270     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements,  
271                    mat.getPaso_SystemMatrix(),     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->FaceElements, mat.getPaso_SystemMatrix(), &_rhs, &_d, &_y, Finley_Assemble_handelShapeMissMatch_Mean_out);
                   &(rhs.getDataC()),  
                                   &(d.getDataC()),&(y.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Mean_out);  
272     checkFinleyError();     checkFinleyError();
273     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements,  
274                    mat.getPaso_SystemMatrix(),     Finley_Assemble_RobinCondition(mesh->Nodes,mesh->ContactElements, mat.getPaso_SystemMatrix(), &_rhs , &_d_contact, &_y_contact ,             Finley_Assemble_handelShapeMissMatch_Step_out);
                   &(rhs.getDataC()),  
                                   &(d_contact.getDataC()),  
                   &(y_contact.getDataC()),  
                                   Finley_Assemble_handelShapeMissMatch_Step_out);  
275     checkFinleyError();     checkFinleyError();
276  }  }
277    
278  //  //
279  // adds linear PDE of second order into the right hand side only  // adds linear PDE of second order into the right hand side only
280  //  //
281  void MeshAdapter::addPDEToRHS( escript::Data& rhs,  void MeshAdapter::addPDEToRHS( escript::Data& rhs, const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const
                      const  escript::Data& X,const  escript::Data& Y, const escript::Data& y, const escript::Data& y_contact) const  
282  {  {
283     Finley_Mesh* mesh=m_finleyMesh.get();     Finley_Mesh* mesh=m_finleyMesh.get();
284    
# Line 441  void MeshAdapter::interpolateOnDomain(es Line 302  void MeshAdapter::interpolateOnDomain(es
302  {  {
303    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());    const MeshAdapter& inDomain=dynamic_cast<const MeshAdapter&>(in.getFunctionSpace().getDomain());
304    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());    const MeshAdapter& targetDomain=dynamic_cast<const MeshAdapter&>(target.getFunctionSpace().getDomain());
305    if (inDomain!=*this)    if (inDomain!=*this)  
306       throw FinleyAdapterException("Error - Illegal domain of interpolant.");      throw FinleyAdapterException("Error - Illegal domain of interpolant.");
307    if (targetDomain!=*this)    if (targetDomain!=*this)
308       throw FinleyAdapterException("Error - Illegal domain of interpolation target.");      throw FinleyAdapterException("Error - Illegal domain of interpolation target.");
309    
310    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
311    switch(in.getFunctionSpace().getTypeCode()) {    switch(in.getFunctionSpace().getTypeCode()) {
# Line 508  void MeshAdapter::interpolateOnDomain(es Line 369  void MeshAdapter::interpolateOnDomain(es
369             break;             break;
370          }          }
371          break;          break;
372       case(DegreesOfFreedom):       case(DegreesOfFreedom):      
373          switch(target.getFunctionSpace().getTypeCode()) {          switch(target.getFunctionSpace().getTypeCode()) {
374             case(ReducedDegreesOfFreedom):             case(ReducedDegreesOfFreedom):
375             case(DegreesOfFreedom):             case(DegreesOfFreedom):
376             case(Nodes):             case(Nodes):
377                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));                Finley_Assemble_CopyNodalData(mesh->Nodes,&(target.getDataC()),&(in.getDataC()));
378                break;                break;
379    #ifndef PASO_MPI
380             case(Elements):             case(Elements):
381                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(in.getDataC()),&(target.getDataC()));
382                break;                break;
# Line 528  void MeshAdapter::interpolateOnDomain(es Line 390  void MeshAdapter::interpolateOnDomain(es
390             case(ContactElementsOne):             case(ContactElementsOne):
391                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));                Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(in.getDataC()),&(target.getDataC()));
392               break;               break;
393    #else
394               /* need to copy Degrees of freedom data to nodal data so that the external values are available */
395               case(Elements):
396               {
397                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
398                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Elements,&(nodeTemp.getDataC()),&(target.getDataC()));
399                  break;
400               }
401               case(FaceElements):
402               {
403                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
404                  Finley_Assemble_interpolate(mesh->Nodes,mesh->FaceElements,&(nodeTemp.getDataC()),&(target.getDataC()));
405                  break;
406               }
407               case(Points):
408               {
409                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
410                  Finley_Assemble_interpolate(mesh->Nodes,mesh->Points,&(nodeTemp.getDataC()),&(target.getDataC()));
411                  break;
412               }
413               case(ContactElementsZero):
414               case(ContactElementsOne):
415               {
416                  escript::Data nodeTemp( in, continuousFunction(asAbstractContinuousDomain()) );
417                  Finley_Assemble_interpolate(mesh->Nodes,mesh->ContactElements,&(nodeTemp.getDataC()),&(target.getDataC()));
418                 break;
419               }
420    #endif
421             default:             default:
422               stringstream temp;               stringstream temp;
423               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();               temp << "Error - Interpolation On Domain: Finley does not know anything about function space type " << target.getFunctionSpace().getTypeCode();
# Line 706  void MeshAdapter::setToGradient(escript: Line 596  void MeshAdapter::setToGradient(escript:
596       throw FinleyAdapterException("Error - Illegal domain of gradient");       throw FinleyAdapterException("Error - Illegal domain of gradient");
597    
598    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
599      escriptDataC nodeDataC;
600    #ifdef PASO_MPI
601      escript::Data nodeTemp( arg, continuousFunction(asAbstractContinuousDomain()) );
602      if( arg.getFunctionSpace().getTypeCode() != Nodes )
603      {
604        Finley_Assemble_CopyNodalData(mesh->Nodes,&(nodeTemp.getDataC()),&(arg.getDataC()));
605        nodeDataC = nodeTemp.getDataC();
606      }
607      else
608        nodeDataC = arg.getDataC();
609    #else
610      nodeDataC = arg.getDataC();
611    #endif
612    switch(grad.getFunctionSpace().getTypeCode()) {    switch(grad.getFunctionSpace().getTypeCode()) {
613         case(Nodes):         case(Nodes):
614            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");            throw FinleyAdapterException("Error - Gradient at nodes is not supported.");
615            break;            break;
616         case(Elements):         case(Elements):
617            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->Elements,&(grad.getDataC()),&nodeDataC);
618            break;            break;
619         case(FaceElements):         case(FaceElements):
620            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->FaceElements,&(grad.getDataC()),&nodeDataC);
621            break;            break;
622         case(Points):         case(Points):
623            throw FinleyAdapterException("Error - Gradient at points is not supported.");            throw FinleyAdapterException("Error - Gradient at points is not supported.");
624            break;            break;
625         case(ContactElementsZero):         case(ContactElementsZero):
626            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
627            break;            break;
628         case(ContactElementsOne):         case(ContactElementsOne):
629            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&(arg.getDataC()));            Finley_Assemble_gradient(mesh->Nodes,mesh->ContactElements,&(grad.getDataC()),&nodeDataC);
630            break;            break;
631         case(DegreesOfFreedom):         case(DegreesOfFreedom):
632            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");            throw FinleyAdapterException("Error - Gradient at degrees of freedom is not supported.");
# Line 783  void MeshAdapter::setToSize(escript::Dat Line 686  void MeshAdapter::setToSize(escript::Dat
686  void MeshAdapter::setNewX(const escript::Data& new_x)  void MeshAdapter::setNewX(const escript::Data& new_x)
687  {  {
688    Finley_Mesh* mesh=m_finleyMesh.get();    Finley_Mesh* mesh=m_finleyMesh.get();
689      escriptDataC tmp;
690    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());    const MeshAdapter& newDomain=dynamic_cast<const MeshAdapter&>(new_x.getFunctionSpace().getDomain());
691    if (newDomain!=*this)    if (newDomain!=*this)
692       throw FinleyAdapterException("Error - Illegal domain of new point locations");       throw FinleyAdapterException("Error - Illegal domain of new point locations");
693    Finley_Mesh_setCoordinates(mesh,&(new_x.getDataC()));    tmp = new_x.getDataC();
694      Finley_Mesh_setCoordinates(mesh,&tmp);
695    checkFinleyError();    checkFinleyError();
696  }  }
697    
# Line 1061  escript::Data MeshAdapter::getSize() con Line 966  escript::Data MeshAdapter::getSize() con
966    
967  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getTagFromSampleNo(int functionSpaceType, int sampleNo) const
968  {  {
969    int* tagList;    int out=0;
970    int numTags;    Finley_Mesh* mesh=m_finleyMesh.get();
971    getTagList(functionSpaceType, &tagList, &numTags);    switch (functionSpaceType) {
972    return tagList[sampleNo];    case(Nodes):
973        out=mesh->Nodes->Tag[sampleNo];
974        break;
975      case(Elements):
976        out=mesh->Elements->Tag[sampleNo];
977        break;
978      case(FaceElements):
979        out=mesh->FaceElements->Tag[sampleNo];
980        break;
981      case(Points):
982        out=mesh->Points->Tag[sampleNo];
983        break;
984      case(ContactElementsZero):
985        out=mesh->ContactElements->Tag[sampleNo];
986        break;
987      case(ContactElementsOne):
988        out=mesh->ContactElements->Tag[sampleNo];
989        break;
990      case(DegreesOfFreedom):
991        break;
992      case(ReducedDegreesOfFreedom):
993        break;
994      default:
995        stringstream temp;
996        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
997        throw FinleyAdapterException(temp.str());
998        break;
999      }
1000      return out;
1001  }  }
   
1002  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const  int MeshAdapter::getReferenceNoFromSampleNo(int functionSpaceType, int sampleNo) const
1003  {  {
1004    int* referenceNoList;    int out=0,i;
1005    int numReferenceNo;    Finley_Mesh* mesh=m_finleyMesh.get();
1006    getReferenceNoList(functionSpaceType, &referenceNoList, &numReferenceNo);    switch (functionSpaceType) {
1007    return referenceNoList[sampleNo];    case(Nodes):
1008        if (mesh->Nodes!=NULL) {
1009          out=mesh->Nodes->Id[sampleNo];
1010          break;
1011        }
1012      case(Elements):
1013        out=mesh->Elements->Id[sampleNo];
1014        break;
1015      case(FaceElements):
1016        out=mesh->FaceElements->Id[sampleNo];
1017        break;
1018      case(Points):
1019        out=mesh->Points->Id[sampleNo];
1020        break;
1021      case(ContactElementsZero):
1022        out=mesh->ContactElements->Id[sampleNo];
1023        break;
1024      case(ContactElementsOne):
1025        out=mesh->ContactElements->Id[sampleNo];
1026        break;
1027      case(DegreesOfFreedom):
1028        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1029           if (mesh->Nodes->degreeOfFreedom[i]==sampleNo) {
1030              out=mesh->Nodes->Id[i];
1031              break;
1032           }
1033        }
1034        break;
1035      case(ReducedDegreesOfFreedom):
1036        for (i=0;i<mesh->Nodes->numNodes; ++i) {
1037           if (mesh->Nodes->reducedDegreeOfFreedom[i]==sampleNo) {
1038              out=mesh->Nodes->Id[i];
1039              break;
1040           }
1041        }
1042        break;
1043      default:
1044        stringstream temp;
1045        temp << "Error - Invalid function space type: " << functionSpaceType << " for domain: " << getDescription();
1046        throw FinleyAdapterException(temp.str());
1047        break;
1048      }
1049      return out;
1050  }  }
1051    
1052  }  // end of namespace  }  // end of namespace

Legend:
Removed from v.532  
changed lines
  Added in v.751

  ViewVC Help
Powered by ViewVC 1.1.26