/[escript]/trunk/ripley/src/RipleyDomain.cpp
ViewVC logotype

Diff of /trunk/ripley/src/RipleyDomain.cpp

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

revision 3755 by caltinay, Thu Jan 5 06:51:31 2012 UTC revision 3785 by caltinay, Wed Jan 25 04:00:33 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 83  string RipleyDomain::functionSpaceTypeAs Line 83  string RipleyDomain::functionSpaceTypeAs
83          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85          case Nodes: return "Ripley_Nodes";          case Nodes: return "Ripley_Nodes";
86          case ReducedNodes: return "Ripley_Reduced_Nodes";          case ReducedNodes: return "Ripley_ReducedNodes";
87          case Elements: return "Ripley_Elements";          case Elements: return "Ripley_Elements";
88          case ReducedElements: return "Ripley_Reduced_Elements";          case ReducedElements: return "Ripley_ReducedElements";
89          case FaceElements: return "Ripley_Face_Elements";          case FaceElements: return "Ripley_FaceElements";
90          case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";          case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91          case Points: return "Ripley_Points";          case Points: return "Ripley_Points";
92          default:          default:
93              break;              break;
# Line 234  bool RipleyDomain::commonFunctionSpace(c Line 234  bool RipleyDomain::commonFunctionSpace(c
234  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235                                                int fsType_target) const                                                int fsType_target) const
236  {  {
237      if (fsType_target != Nodes &&      if (!isValidFunctionSpaceType(fsType_target)) {
             fsType_target != ReducedNodes &&  
             fsType_target != ReducedDegreesOfFreedom &&  
             fsType_target != DegreesOfFreedom &&  
             fsType_target != Elements &&  
             fsType_target != ReducedElements &&  
             fsType_target != FaceElements &&  
             fsType_target != ReducedFaceElements &&  
             fsType_target != Points) {  
238          stringstream msg;          stringstream msg;
239          msg << "probeInterpolationOnDomain(): Invalid functionspace type "          msg << "probeInterpolationOnDomain(): Invalid functionspace type "
240              << fsType_target << " for " << getDescription();              << fsType_target << " for " << getDescription();
# Line 260  bool RipleyDomain::probeInterpolationOnD Line 252  bool RipleyDomain::probeInterpolationOnD
252          case Elements:          case Elements:
253              return (fsType_target==Elements ||              return (fsType_target==Elements ||
254                      fsType_target==ReducedElements);                      fsType_target==ReducedElements);
         case ReducedElements:  
             return (fsType_target==ReducedElements);  
255          case FaceElements:          case FaceElements:
256              return (fsType_target==FaceElements ||              return (fsType_target==FaceElements ||
257                      fsType_target==ReducedFaceElements);                      fsType_target==ReducedFaceElements);
258            case ReducedElements:
259          case ReducedFaceElements:          case ReducedFaceElements:
             return (fsType_target==ReducedFaceElements);  
260          case Points:          case Points:
261              return (fsType_target==Points);              return (fsType_target==fsType_source);
262    
263          default: {          default: {
264              stringstream msg;              stringstream msg;
# Line 307  void RipleyDomain::interpolateOnDomain(e Line 297  void RipleyDomain::interpolateOnDomain(e
297          throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");          throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
298      } else if ((inFS==Elements && outFS==ReducedElements)      } else if ((inFS==Elements && outFS==ReducedElements)
299              || (inFS==FaceElements && outFS==ReducedFaceElements)) {              || (inFS==FaceElements && outFS==ReducedFaceElements)) {
300          averageData(target, *const_cast<escript::Data*>(&in));          if (in.actsExpanded())
301                averageData(target, *const_cast<escript::Data*>(&in));
302            else
303                copyData(target, *const_cast<escript::Data*>(&in));
304      } else {      } else {
305          switch (inFS) {          switch (inFS) {
306              case Nodes:              case Nodes:
# Line 368  void RipleyDomain::interpolateOnDomain(e Line 361  void RipleyDomain::interpolateOnDomain(e
361                          if (getMPISize()==1) {                          if (getMPISize()==1) {
362                              interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);                              interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363                          } else {                          } else {
364                              escript::Data contIn=escript::Data(in, continuousFunction(*this));                              escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366                              interpolateNodesOnElements(target, contIn, outFS==ReducedElements);                              interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367                          }                          }
368                          break;                          break;
# Line 378  void RipleyDomain::interpolateOnDomain(e Line 372  void RipleyDomain::interpolateOnDomain(e
372                          if (getMPISize()==1) {                          if (getMPISize()==1) {
373                              interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);                              interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374                          } else {                          } else {
375                              escript::Data contIn=escript::Data(in, continuousFunction(*this));                              escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376                              interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);                                       escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378                          }                          }
379                          break;                          break;
380    
# Line 427  void RipleyDomain::setToX(escript::Data& Line 422  void RipleyDomain::setToX(escript::Data&
422      }      }
423  }  }
424    
425    void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426    {
427        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
428                *(arg.getFunctionSpace().getDomain()));
429        if (argDomain != *this)
430            throw RipleyException("setToGradient: Illegal domain of gradient argument");
431        const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
432                *(grad.getFunctionSpace().getDomain()));
433        if (gradDomain != *this)
434            throw RipleyException("setToGradient: Illegal domain of gradient");
435    
436        switch (grad.getFunctionSpace().getTypeCode()) {
437            case Elements:
438            case ReducedElements:
439            case FaceElements:
440            case ReducedFaceElements:
441                break;
442            default: {
443                stringstream msg;
444                msg << "setToGradient(): not supported for "
445                    << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
446                throw RipleyException(msg.str());
447            }
448        }
449    
450        if (getMPISize()>1) {
451            if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
452                escript::Data contArg(arg, escript::continuousFunction(*this));
453                assembleGradient(grad, contArg);
454            } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
455                escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
456                assembleGradient(grad, contArg);
457            } else {
458                assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459            }
460        } else {
461            assembleGradient(grad, *const_cast<escript::Data*>(&arg));
462        }
463    }
464    
465    void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
466    {
467        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
468                *(arg.getFunctionSpace().getDomain()));
469        if (argDomain != *this)
470            throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
471    
472        switch (arg.getFunctionSpace().getTypeCode()) {
473            case Nodes:
474            case ReducedNodes:
475            case DegreesOfFreedom:
476            case ReducedDegreesOfFreedom:
477                {
478                    escript::Data funcArg(arg, escript::function(*this));
479                    assembleIntegrate(integrals, funcArg);
480                }
481                break;
482            case Elements:
483            case ReducedElements:
484            case FaceElements:
485            case ReducedFaceElements:
486                assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
487                break;
488            default: {
489                stringstream msg;
490                msg << "setToIntegrals(): not supported for "
491                    << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
492                throw RipleyException(msg.str());
493            }
494        }
495    
496    }
497    
498  bool RipleyDomain::isCellOriented(int fsType) const  bool RipleyDomain::isCellOriented(int fsType) const
499  {  {
500      switch(fsType) {      switch(fsType) {
501          case Nodes:          case Nodes:
502            case ReducedNodes:
503          case DegreesOfFreedom:          case DegreesOfFreedom:
504          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
505              return false;              return false;
# Line 461  bool RipleyDomain::canTag(int fsType) co Line 530  bool RipleyDomain::canTag(int fsType) co
530          case DegreesOfFreedom:          case DegreesOfFreedom:
531          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
532          case Points:          case Points:
533            case ReducedNodes:
534              return false;              return false;
535          default:          default:
536              break;              break;
# Line 620  escript::ASM_ptr RipleyDomain::newSystem Line 690  escript::ASM_ptr RipleyDomain::newSystem
690      bool reduceColOrder=false;      bool reduceColOrder=false;
691      // is the domain right?      // is the domain right?
692      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
693      if (row_domain!=*this)      if (row_domain != *this)
694          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
695      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
696      if (col_domain!=*this)      if (col_domain != *this)
697          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
698      // is the function space type right?      // is the function space type right?
699      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
# Line 640  escript::ASM_ptr RipleyDomain::newSystem Line 710  escript::ASM_ptr RipleyDomain::newSystem
710      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
711              row_blocksize, column_blocksize, FALSE);              row_blocksize, column_blocksize, FALSE);
712      paso::checkPasoError();      paso::checkPasoError();
     Paso_SystemMatrixPattern_free(pattern);  
713      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
714                  row_functionspace, column_blocksize, column_functionspace));                  row_functionspace, column_blocksize, column_functionspace));
715      return sma;      return sma;
# Line 664  void RipleyDomain::addPDEToSystem( Line 733  void RipleyDomain::addPDEToSystem(
733      if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))      if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
734          throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");          throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
735    
736      //TODO: more input param checks (shape, function space etc)      int fsType=UNKNOWN;
737        fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
738        fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
739        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
740        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
741        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
742        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
743        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
744        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
745        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
746        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
747        if (fsType==UNKNOWN)
748            return;
749    
750        const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
751    
752        //TODO: more input param checks (shape, etc)
753    
754      Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();      Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
755    
# Line 678  void RipleyDomain::addPDEToSystem( Line 763  void RipleyDomain::addPDEToSystem(
763          throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");          throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
764      //TODO: more system matrix checks      //TODO: more system matrix checks
765    
766      if (numEq==1)      if (numEq==1) {
767          assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);          if (reducedOrder) {
768      else              assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
769          assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);              assemblePDEBoundarySingleReduced(S, rhs, d, y);
770            } else {
771                assemblePDESingle(S, rhs, A, B, C, D, X, Y);
772                assemblePDEBoundarySingle(S, rhs, d, y);
773            }
774        } else {
775            if (reducedOrder) {
776                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
777                assemblePDEBoundarySystemReduced(S, rhs, d, y);
778            } else {
779                assemblePDESystem(S, rhs, A, B, C, D, X, Y);
780                assemblePDEBoundarySystem(S, rhs, d, y);
781            }
782        }
783    }
784    
785    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
786            const escript::Data& Y, const escript::Data& y,
787            const escript::Data& y_contact, const escript::Data& y_dirac) const
788    {
789        if (!y_contact.isEmpty())
790            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
791    
792        if (rhs.isEmpty()) {
793            if (!X.isEmpty() || !Y.isEmpty())
794                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
795            else
796                return;
797        }
798    
799        if (rhs.getDataPointSize() == 1) {
800            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
801            assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
802        } else {
803            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
804            assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
805        }
806  }  }
807    
808  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
# Line 780  void RipleyDomain::updateTagsInUse(int f Line 901  void RipleyDomain::updateTagsInUse(int f
901      }      }
902  }  }
903    
904    //protected
905    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
906            const IndexVector& index, const dim_t M, const dim_t N) const
907    {
908        // paso will manage the memory
909        index_t* indexC = MEMALLOC(index.size(), index_t);
910        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
911        copy(index.begin(), index.end(), indexC);
912        copy(ptr.begin(), ptr.end(), ptrC);
913        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
914    }
915    
916    //protected
917    Paso_Pattern* RipleyDomain::createMainPattern() const
918    {
919        IndexVector ptr(1,0);
920        IndexVector index;
921    
922        for (index_t i=0; i<getNumDOF(); i++) {
923            // add the DOF itself
924            index.push_back(i);
925            const dim_t num=insertNeighbourNodes(index, i);
926            ptr.push_back(ptr.back()+num+1);
927        }
928    
929        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
930    }
931    
932    //protected
933    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
934            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
935    {
936        IndexVector ptr(1,0);
937        IndexVector index;
938        for (index_t i=0; i<getNumDOF(); i++) {
939            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
940            ptr.push_back(ptr.back()+colIndices[i].size());
941        }
942    
943        const dim_t M=ptr.size()-1;
944        *colPattern=createPasoPattern(ptr, index, M, N);
945    
946        IndexVector rowPtr(1,0);
947        IndexVector rowIndex;
948        for (dim_t id=0; id<N; id++) {
949            dim_t n=0;
950            for (dim_t i=0; i<M; i++) {
951                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
952                    if (index[j]==id) {
953                        rowIndex.push_back(i);
954                        n++;
955                        break;
956                    }
957                }
958            }
959            rowPtr.push_back(rowPtr.back()+n);
960        }
961    
962        // M and N are now reversed
963        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
964    }
965    
966    //protected
967    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
968           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
969           dim_t num_Sol, const vector<double>& array) const
970    {
971        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
972        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
973        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
974        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
975    
976        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
977        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
978        double* mainBlock_val = mat->mainBlock->val;
979        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
980        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
981        double* col_coupleBlock_val = mat->col_coupleBlock->val;
982        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
983        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
984        double* row_coupleBlock_val = mat->row_coupleBlock->val;
985    
986        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
987            // down columns of array
988            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
989                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
990                // only look at the matrix rows stored on this processor
991                if (i_row < numMyRows) {
992                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
993                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
994                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
995                            if (i_col < numMyCols) {
996                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
997                                    if (mainBlock_index[k] == i_col) {
998                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
999                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1000                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1001                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1002                                                mainBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1003                                            }
1004                                        }
1005                                        break;
1006                                    }
1007                                }
1008                            } else {
1009                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1010                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1011                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1012                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1013                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1014                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1015                                                col_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1016                                            }
1017                                        }
1018                                        break;
1019                                    }
1020                                }
1021                            }
1022                        }
1023                    }
1024                } else {
1025                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1026                        // across rows of array
1027                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1028                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1029                            if (i_col < numMyCols) {
1030                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1031                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1032                                {
1033                                    if (row_coupleBlock_index[k] == i_col) {
1034                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1035                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1036                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1037                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1038                                            row_coupleBlock_val[k*mat->block_size+ir+mat->row_block_size*ic] += array[INDEX4(i_Eq, i_Sol, k_Eq, k_Sol, num_Eq, num_Sol, nodes_Eq.size())];
1039                                            }
1040                                        }
1041                                        break;
1042                                    }
1043                                }
1044                            }
1045                        }
1046                    }
1047                }
1048            }
1049        }
1050    }
1051    
1052  //  //
1053  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
1054  //  //
# Line 825  void RipleyDomain::setToSize(escript::Da Line 1094  void RipleyDomain::setToSize(escript::Da
1094      throw RipleyException("setToSize() not implemented");      throw RipleyException("setToSize() not implemented");
1095  }  }
1096    
 void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const  
 {  
     throw RipleyException("setToGradient() not implemented");  
 }  
   
 void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  
 {  
     throw RipleyException("setToIntegrals() not implemented");  
 }  
   
1097  bool RipleyDomain::ownSample(int fsType, index_t id) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1098  {  {
1099      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1100  }  }
1101    
 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  
         const escript::Data& D, const escript::Data& d,  
         const escript::Data& d_dirac, const bool useHRZ) const  
 {  
     throw RipleyException("addPDEToLumpedSystem() not implemented");  
 }  
   
 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,  
         const escript::Data& Y, const escript::Data& y,  
         const escript::Data& y_contact, const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToRHS() not implemented");  
 }  
   
1102  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1103          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1104          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 904  IndexVector RipleyDomain::getNodeDistrib Line 1149  IndexVector RipleyDomain::getNodeDistrib
1149      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1150  }  }
1151    
1152    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1153    {
1154        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1155    }
1156    
1157  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1158  {  {
1159      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 929  dim_t RipleyDomain::getNumDOF() const Line 1179  dim_t RipleyDomain::getNumDOF() const
1179      throw RipleyException("getNumDOF() not implemented");      throw RipleyException("getNumDOF() not implemented");
1180  }  }
1181    
1182    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1183    {
1184        throw RipleyException("insertNeighbourNodes() not implemented");
1185    }
1186    
1187  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1188  {  {
1189      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1190  }  }
1191    
1192    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1193    {
1194        throw RipleyException("assembleGradient() not implemented");
1195    }
1196    
1197    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1198    {
1199        throw RipleyException("assembleIntegrate() not implemented");
1200    }
1201    
1202  void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,  void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1203          const escript::Data& A, const escript::Data& B, const escript::Data& C,          const escript::Data& A, const escript::Data& B, const escript::Data& C,
1204          const escript::Data& D, const escript::Data& X, const escript::Data& Y,          const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1205  {  {
1206      throw RipleyException("assemblePDESingle() not implemented");      throw RipleyException("assemblePDESingle() not implemented");
1207  }  }
1208    
1209    void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1210          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1211    {
1212        throw RipleyException("assemblePDEBoundarySingle() not implemented");
1213    }
1214    
1215    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1216            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1217            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1218            const escript::Data& Y) const
1219    {
1220        throw RipleyException("assemblePDESingleReduced() not implemented");
1221    }
1222    
1223    void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1224          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1225    {
1226        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1227    }
1228    
1229  void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,  void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1230          const escript::Data& A, const escript::Data& B, const escript::Data& C,          const escript::Data& A, const escript::Data& B, const escript::Data& C,
1231          const escript::Data& D, const escript::Data& X, const escript::Data& Y,          const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
         const escript::Data& d, const escript::Data& y) const  
1232  {  {
1233      throw RipleyException("assemblePDESystem() not implemented");      throw RipleyException("assemblePDESystem() not implemented");
1234  }  }
1235    
1236    void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1237          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1238    {
1239        throw RipleyException("assemblePDEBoundarySystem() not implemented");
1240    }
1241    
1242    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1243            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1244            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1245            const escript::Data& Y) const
1246    {
1247        throw RipleyException("assemblePDESystemReduced() not implemented");
1248    }
1249    
1250    void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1251          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1252    {
1253        throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1254    }
1255    
1256  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1257  {  {
1258      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");

Legend:
Removed from v.3755  
changed lines
  Added in v.3785

  ViewVC Help
Powered by ViewVC 1.1.26