/[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 3746 by caltinay, Thu Dec 15 00:02:22 2011 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 61  bool RipleyDomain::operator==(const Abst Line 61  bool RipleyDomain::operator==(const Abst
61  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62  {  {
63      switch (fsType) {      switch (fsType) {
64            case DegreesOfFreedom:
65            case ReducedDegreesOfFreedom:
66          case Nodes:          case Nodes:
67          case ReducedNodes:          case ReducedNodes:
68          case Elements:          case Elements:
# Line 78  bool RipleyDomain::isValidFunctionSpaceT Line 80  bool RipleyDomain::isValidFunctionSpaceT
80  string RipleyDomain::functionSpaceTypeAsString(int fsType) const  string RipleyDomain::functionSpaceTypeAsString(int fsType) const
81  {  {
82      switch (fsType) {      switch (fsType) {
83            case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84            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 98  pair<int,int> RipleyDomain::getDataShape Line 102  pair<int,int> RipleyDomain::getDataShape
102          case Nodes:          case Nodes:
103          case ReducedNodes: //FIXME: reduced          case ReducedNodes: //FIXME: reduced
104              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
105            case DegreesOfFreedom:
106            case ReducedDegreesOfFreedom: //FIXME: reduced
107                return pair<int,int>(1, getNumDOF());
108          case Elements:          case Elements:
109              return pair<int,int>(ptsPerSample, getNumElements());              return pair<int,int>(ptsPerSample, getNumElements());
110          case FaceElements:          case FaceElements:
# Line 134  bool RipleyDomain::commonFunctionSpace(c Line 141  bool RipleyDomain::commonFunctionSpace(c
141     /*     /*
142      The idea is to use equivalence classes (i.e. types which can be      The idea is to use equivalence classes (i.e. types which can be
143      interpolated back and forth):      interpolated back and forth):
144      class 0: Nodes      class 0: DOF <-> Nodes
145      class 1: ReducedNodes      class 1: ReducedDOF <-> ReducedNodes
146      class 2: Points      class 2: Points
147      class 3: Elements      class 3: Elements
148      class 4: ReducedElements      class 4: ReducedElements
# Line 148  bool RipleyDomain::commonFunctionSpace(c Line 155  bool RipleyDomain::commonFunctionSpace(c
155      line 0: class 2      line 0: class 2
156      line 1: class 3,4      line 1: class 3,4
157      line 2: class 5,6      line 2: class 5,6
158    
159        For classes with multiple members (eg class 1) we have vars to record if
160        there is at least one instance. e.g. hasnodes is true if we have at least
161        one instance of Nodes.
162      */      */
163      if (fs.empty())      if (fs.empty())
164          return false;          return false;
165      vector<bool> hasclass(7, false);      vector<bool> hasclass(7, false);
166      vector<int> hasline(3, 0);      vector<int> hasline(3, 0);
167        bool hasnodes=false;
168        bool hasrednodes=false;
169      for (size_t i=0; i<fs.size(); ++i) {      for (size_t i=0; i<fs.size(); ++i) {
170          switch (fs[i]) {          switch (fs[i]) {
171              case Nodes:              case Nodes: hasnodes=true; // fall through
172                case DegreesOfFreedom:
173                  hasclass[0]=true;                  hasclass[0]=true;
174                  break;                  break;
175              case ReducedNodes:              case ReducedNodes: hasrednodes=true; // fall through
176                case ReducedDegreesOfFreedom:
177                  hasclass[1]=true;                  hasclass[1]=true;
178                  break;                  break;
179              case Points:              case Points:
# Line 166  bool RipleyDomain::commonFunctionSpace(c Line 181  bool RipleyDomain::commonFunctionSpace(c
181                  hasclass[2]=true;                  hasclass[2]=true;
182                  break;                  break;
183              case Elements:              case Elements:
184                    hasclass[3]=true;
185                  hasline[1]=1;                  hasline[1]=1;
186                  break;                  break;
187              case ReducedElements:              case ReducedElements:
# Line 173  bool RipleyDomain::commonFunctionSpace(c Line 189  bool RipleyDomain::commonFunctionSpace(c
189                  hasline[1]=1;                  hasline[1]=1;
190                  break;                  break;
191              case FaceElements:              case FaceElements:
192                    hasclass[5]=true;
193                  hasline[2]=1;                  hasline[2]=1;
194                  break;                  break;
195              case ReducedFaceElements:              case ReducedFaceElements:
# Line 190  bool RipleyDomain::commonFunctionSpace(c Line 207  bool RipleyDomain::commonFunctionSpace(c
207      if (numLines > 1)      if (numLines > 1)
208          return false;          return false;
209      else if (numLines==1) {      else if (numLines==1) {
210            // we have points
211          if (hasline[0]==1)          if (hasline[0]==1)
212              resultcode=Points;              resultcode=Points;
213          else if (hasline[1]==1) {          else if (hasline[1]==1) {
# Line 205  bool RipleyDomain::commonFunctionSpace(c Line 223  bool RipleyDomain::commonFunctionSpace(c
223          }          }
224      } else { // numLines==0      } else { // numLines==0
225          if (hasclass[1])          if (hasclass[1])
226              resultcode=ReducedNodes;              // something from class 1
227          else              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228              resultcode=Nodes;          else // something from class 0
229                resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230      }      }
231      return true;      return true;
232  }  }
# Line 215  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 != 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 230  bool RipleyDomain::probeInterpolationOnD Line 243  bool RipleyDomain::probeInterpolationOnD
243    
244      switch (fsType_source) {      switch (fsType_source) {
245          case Nodes:          case Nodes:
246            case DegreesOfFreedom:
247              return true;              return true;
248          case ReducedNodes:          case ReducedNodes:
249              return (fsType_target != Nodes);          case ReducedDegreesOfFreedom:
250                return (fsType_target != Nodes &&
251                        fsType_target != DegreesOfFreedom);
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 278  void RipleyDomain::interpolateOnDomain(e Line 292  void RipleyDomain::interpolateOnDomain(e
292      if (inFS==outFS) {      if (inFS==outFS) {
293          copyData(target, *const_cast<escript::Data*>(&in));          copyData(target, *const_cast<escript::Data*>(&in));
294      // not allowed: reduced->non-reduced      // not allowed: reduced->non-reduced
295      } else if (inFS==ReducedNodes && outFS==Nodes) {      } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
296                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
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:
307              case ReducedNodes:              case ReducedNodes: //FIXME: reduced
308                  switch (outFS) {                  switch (outFS) {
309                        case DegreesOfFreedom:
310                        case ReducedDegreesOfFreedom: //FIXME: reduced
311                            if (getMPISize()==1)
312                                copyData(target, *const_cast<escript::Data*>(&in));
313                            else
314                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
315                            break;
316    
317                      case Nodes:                      case Nodes:
318                      case ReducedNodes: //FIXME: reduced                      case ReducedNodes: //FIXME: reduced
319                          copyData(target, *const_cast<escript::Data*>(&in));                          copyData(target, *const_cast<escript::Data*>(&in));
# Line 313  void RipleyDomain::interpolateOnDomain(e Line 339  void RipleyDomain::interpolateOnDomain(e
339                          throw RipleyException(msg.str());                          throw RipleyException(msg.str());
340                  }                  }
341                  break;                  break;
342    
343                case DegreesOfFreedom:
344                case ReducedDegreesOfFreedom: //FIXME: reduced
345                    switch (outFS) {
346                        case Nodes:
347                        case ReducedNodes: //FIXME: reduced
348                            if (getMPISize()==1)
349                                copyData(target, *const_cast<escript::Data*>(&in));
350                            else
351                                dofToNodes(target, *const_cast<escript::Data*>(&in));
352                            break;
353    
354                        case DegreesOfFreedom:
355                        case ReducedDegreesOfFreedom: //FIXME: reduced
356                            copyData(target, *const_cast<escript::Data*>(&in));
357                            break;
358    
359                        case Elements:
360                        case ReducedElements:
361                            if (getMPISize()==1) {
362                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363                            } else {
364                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367                            }
368                            break;
369    
370                        case FaceElements:
371                        case ReducedFaceElements:
372                            if (getMPISize()==1) {
373                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374                            } else {
375                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376                                         escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378                            }
379                            break;
380    
381                        default:
382                            throw RipleyException(msg.str());
383                    }
384                    break;
385              default:              default:
386                  throw RipleyException(msg.str());                  throw RipleyException(msg.str());
387          }          }
# Line 353  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:
504            case ReducedDegreesOfFreedom:
505              return false;              return false;
506          case Elements:          case Elements:
507          case FaceElements:          case FaceElements:
# Line 382  bool RipleyDomain::canTag(int fsType) co Line 527  bool RipleyDomain::canTag(int fsType) co
527          case FaceElements:          case FaceElements:
528          case ReducedFaceElements:          case ReducedFaceElements:
529              return true;              return true;
530            case DegreesOfFreedom:
531            case ReducedDegreesOfFreedom:
532          case Points:          case Points:
533            case ReducedNodes:
534              return false;              return false;
535          default:          default:
536              break;              break;
# Line 542  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()==ReducedNodes)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
700          reduceRowOrder=true;          reduceRowOrder=true;
701      else if (row_functionspace.getTypeCode()!=Nodes)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
702          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
703      if (column_functionspace.getTypeCode()==ReducedNodes)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
704          reduceColOrder=true;          reduceColOrder=true;
705      else if (column_functionspace.getTypeCode()!=Nodes)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
706          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
707    
708      // generate matrix      // generate matrix
# Line 562  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;
716  }  }
717    
718    void RipleyDomain::addPDEToSystem(
719            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
720            const escript::Data& A, const escript::Data& B, const escript::Data& C,
721            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
722            const escript::Data& d, const escript::Data& y,
723            const escript::Data& d_contact, const escript::Data& y_contact,
724            const escript::Data& d_dirac,const escript::Data& y_dirac) const
725    {
726        if (!d_contact.isEmpty() || !y_contact.isEmpty())
727            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
728    
729        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
730        if (!sma)
731            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
732    
733        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
734            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
735    
736        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();
755    
756        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
757            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
758    
759        const int numEq=S->logical_row_block_size;
760        const int numComp=S->logical_col_block_size;
761    
762        if (numEq != numComp)
763            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
764        //TODO: more system matrix checks
765    
766        if (numEq==1) {
767            if (reducedOrder) {
768                assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
769                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)
809  {  {
810      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
# Line 664  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 709  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::addPDEToSystem(  
         escript::AbstractSystemMatrix& mat, escript::Data& rhs,  
         const escript::Data& A, const escript::Data& B, const escript::Data& C,  
         const escript::Data& D, const escript::Data& X, const escript::Data& Y,  
         const escript::Data& d, const escript::Data& y,  
         const escript::Data& d_contact, const escript::Data& y_contact,  
         const escript::Data& d_dirac,const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToSystem() not implemented");  
 }  
   
 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 799  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 819  dim_t RipleyDomain::getNumNodes() const Line 1174  dim_t RipleyDomain::getNumNodes() const
1174      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1175  }  }
1176    
1177    dim_t RipleyDomain::getNumDOF() const
1178    {
1179        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,
1203            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
1205    {
1206        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,
1230            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
1232    {
1233        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");
# Line 834  void RipleyDomain::interpolateNodesOnFac Line 1263  void RipleyDomain::interpolateNodesOnFac
1263      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1264  }  }
1265    
1266    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1267    {
1268        throw RipleyException("nodesToDOF() not implemented");
1269    }
1270    
1271    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1272    {
1273        throw RipleyException("dofToNodes() not implemented");
1274    }
1275    
1276  } // end of namespace ripley  } // end of namespace ripley
1277    

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

  ViewVC Help
Powered by ViewVC 1.1.26