/[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 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC revision 3761 by caltinay, Mon Jan 9 06:50:33 2012 UTC
# Line 47  RipleyDomain::~RipleyDomain() Line 47  RipleyDomain::~RipleyDomain()
47      Esys_MPIInfo_free(m_mpiInfo);      Esys_MPIInfo_free(m_mpiInfo);
48  }  }
49    
50    bool RipleyDomain::operator==(const AbstractDomain& other) const
51    {
52        const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53        if (o) {
54            return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55                    && m_elementTags==o->m_elementTags
56                    && m_faceTags==o->m_faceTags);
57        }
58        return false;
59    }
60    
61  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62  {  {
63      switch (fsType) {      switch (fsType) {
# Line 89  pair<int,int> RipleyDomain::getDataShape Line 100  pair<int,int> RipleyDomain::getDataShape
100      const int ptsPerSample = (m_numDim==2 ? 4 : 8);      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101      switch (fsType) {      switch (fsType) {
102          case Nodes:          case Nodes:
103          case ReducedNodes:          case ReducedNodes: //FIXME: reduced
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
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 101  pair<int,int> RipleyDomain::getDataShape Line 113  pair<int,int> RipleyDomain::getDataShape
113              return pair<int,int>(1, getNumElements());              return pair<int,int>(1, getNumElements());
114          case ReducedFaceElements:          case ReducedFaceElements:
115              return pair<int,int>(1, getNumFaceElements());              return pair<int,int>(1, getNumFaceElements());
             /*  
116          case Points:          case Points:
117              return pair<int,int>(1, getNumPoints());              return pair<int,int>(1, 0); //FIXME: dirac
             */  
118          default:          default:
119              break;              break;
120      }      }
121    
122      stringstream msg;      stringstream msg;
123      msg << "getDataShape(): Unsupported function space type "      msg << "getDataShape(): Unsupported function space type " << fsType
124          << functionSpaceTypeAsString(fsType) << " for " << getDescription();          << " for " << getDescription();
125      throw RipleyException(msg.str());      throw RipleyException(msg.str());
126  }  }
127    
# Line 131  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 1: DOF <-> Nodes      class 0: DOF <-> Nodes
145      class 2: ReducedDOF <-> ReducedNodes      class 1: ReducedDOF <-> ReducedNodes
146      class 3: Points      class 2: Points
147      class 4: Elements      class 3: Elements
148      class 5: ReducedElements      class 4: ReducedElements
149      class 6: FaceElements      class 5: FaceElements
150      class 7: ReducedFaceElements      class 6: ReducedFaceElements
     class 8: ContactElementZero <-> ContactElementOne  
     class 9: ReducedContactElementZero <-> ReducedContactElementOne  
151    
152      There is also a set of lines. Interpolation is possible down a line but not      There is also a set of lines. Interpolation is possible down a line but not
153      between lines.      between lines.
154      class 1 and 2 belong to all lines so aren't considered.      class 0 and 1 belong to all lines so aren't considered.
155      line 0: class 3      line 0: class 2
156      line 1: class 4,5      line 1: class 3,4
157      line 2: class 6,7      line 2: class 5,6
     line 3: class 8,9  
158    
159      For classes with multiple members (eg class 2) we have vars to record if      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      there is at least one instance. e.g. hasnodes is true if we have at least
161      one instance of Nodes.      one instance of Nodes.
162      */      */
163      if (fs.empty())      if (fs.empty())
164          return false;          return false;
165      vector<int> hasclass(10);      vector<bool> hasclass(7, false);
166      vector<int> hasline(4);      vector<int> hasline(3, 0);
167      bool hasnodes=false;      bool hasnodes=false;
168      bool hasrednodes=false;      bool hasrednodes=false;
169      for (int 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: hasnodes=true; // no break is deliberate              case Nodes: hasnodes=true; // fall through
172              case DegreesOfFreedom:              case DegreesOfFreedom:
173                  hasclass[1]=1;                  hasclass[0]=true;
174                  break;                  break;
175              case ReducedNodes: hasrednodes=true; // no break is deliberate              case ReducedNodes: hasrednodes=true; // fall through
176              case ReducedDegreesOfFreedom:              case ReducedDegreesOfFreedom:
177                  hasclass[2]=1;                  hasclass[1]=true;
178                  break;                  break;
179              case Points:              case Points:
180                  hasline[0]=1;                  hasline[0]=1;
181                  hasclass[3]=1;                  hasclass[2]=true;
182                  break;                  break;
183              case Elements:              case Elements:
184                  hasclass[4]=1;                  hasclass[3]=true;
185                  hasline[1]=1;                  hasline[1]=1;
186                  break;                  break;
187              case ReducedElements:              case ReducedElements:
188                  hasclass[5]=1;                  hasclass[4]=true;
189                  hasline[1]=1;                  hasline[1]=1;
190                  break;                  break;
191              case FaceElements:              case FaceElements:
192                  hasclass[6]=1;                  hasclass[5]=true;
193                  hasline[2]=1;                  hasline[2]=1;
194                  break;                  break;
195              case ReducedFaceElements:              case ReducedFaceElements:
196                  hasclass[7]=1;                  hasclass[6]=true;
197                  hasline[2]=1;                  hasline[2]=1;
198                  break;                  break;
199              default:              default:
200                  return false;                  return false;
201          }          }
202      }      }
203      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int numLines=hasline[0]+hasline[1]+hasline[2];
204    
205      // fail if we have more than one leaf group      // fail if we have more than one leaf group
206      // = there are at least two branches we can't interpolate between      // = there are at least two branches we can't interpolate between
207      if (totlines>1)      if (numLines > 1)
208          return false;          return false;
209      else if (totlines==1) {      else if (numLines==1) {
210          // we have points          // 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) {
214              if (hasclass[5]==1)              if (hasclass[4])
215                  resultcode=ReducedElements;                  resultcode=ReducedElements;
216              else              else
217                  resultcode=Elements;                  resultcode=Elements;
218          } else if (hasline[2]==1) {          } else { // hasline[2]==1
219              if (hasclass[7]==1)              if (hasclass[6])
220                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
221              else              else
222                  resultcode=FaceElements;                  resultcode=FaceElements;
223          } else          }
224              throw RipleyException("Internal Ripley Error!");      } else { // numLines==0
225      } else { // totlines==0          if (hasclass[1])
226          if (hasclass[2]==1)              // something from class 1
             // something from class 2  
227              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228          else // something from class 1          else // something from class 0
229              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230      }      }
231      return true;      return true;
232  }  }
233    
234    bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235                                                  int fsType_target) const
236    {
237        if (fsType_target != Nodes &&
238                fsType_target != ReducedNodes &&
239                fsType_target != ReducedDegreesOfFreedom &&
240                fsType_target != DegreesOfFreedom &&
241                fsType_target != Elements &&
242                fsType_target != ReducedElements &&
243                fsType_target != FaceElements &&
244                fsType_target != ReducedFaceElements &&
245                fsType_target != Points) {
246            stringstream msg;
247            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
248                << fsType_target << " for " << getDescription();
249            throw RipleyException(msg.str());
250        }
251    
252        switch (fsType_source) {
253            case Nodes:
254            case DegreesOfFreedom:
255                return true;
256            case ReducedNodes:
257            case ReducedDegreesOfFreedom:
258                return (fsType_target != Nodes &&
259                        fsType_target != DegreesOfFreedom);
260            case Elements:
261                return (fsType_target==Elements ||
262                        fsType_target==ReducedElements);
263            case ReducedElements:
264                return (fsType_target==ReducedElements);
265            case FaceElements:
266                return (fsType_target==FaceElements ||
267                        fsType_target==ReducedFaceElements);
268            case ReducedFaceElements:
269                return (fsType_target==ReducedFaceElements);
270            case Points:
271                return (fsType_target==Points);
272    
273            default: {
274                stringstream msg;
275                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
276                    << fsType_source << " for " << getDescription();
277                throw RipleyException(msg.str());
278            }
279        }
280    }
281    
282  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
283                                         const escript::Data& in) const                                         const escript::Data& in) const
284  {  {
# Line 241  void RipleyDomain::interpolateOnDomain(e Line 295  void RipleyDomain::interpolateOnDomain(e
295          << " -> "          << " -> "
296          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
297    
298      switch (in.getFunctionSpace().getTypeCode()) {      const int inFS = in.getFunctionSpace().getTypeCode();
299          case Nodes:      const int outFS = target.getFunctionSpace().getTypeCode();
         case ReducedNodes:  
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     // FIXME  
                     copyNodalData(target, *const_cast<escript::Data*>(&in));  
                     break;  
   
                 case Elements:  
                     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);  
                     break;  
   
                 case ReducedElements:  
                     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);  
                     break;  
300    
301                  case FaceElements:      // simplest case: 1:1 copy
302                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);      if (inFS==outFS) {
303                      break;          copyData(target, *const_cast<escript::Data*>(&in));
304        // not allowed: reduced->non-reduced
305                  case ReducedFaceElements:      } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);              && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307                      break;          throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
308        } else if ((inFS==Elements && outFS==ReducedElements)
309                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
310            averageData(target, *const_cast<escript::Data*>(&in));
311        } else {
312            switch (inFS) {
313                case Nodes:
314                case ReducedNodes: //FIXME: reduced
315                    switch (outFS) {
316                        case DegreesOfFreedom:
317                        case ReducedDegreesOfFreedom: //FIXME: reduced
318                            if (getMPISize()==1)
319                                copyData(target, *const_cast<escript::Data*>(&in));
320                            else
321                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
322                            break;
323    
324                        case Nodes:
325                        case ReducedNodes: //FIXME: reduced
326                            copyData(target, *const_cast<escript::Data*>(&in));
327                            break;
328    
329                        case Elements:
330                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331                            break;
332    
333                        case ReducedElements:
334                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335                            break;
336    
337                        case FaceElements:
338                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339                            break;
340    
341                        case ReducedFaceElements:
342                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
343                            break;
344    
345                        default:
346                            throw RipleyException(msg.str());
347                    }
348                    break;
349    
350                  default:              case DegreesOfFreedom:
351                      throw RipleyException(msg.str());              case ReducedDegreesOfFreedom: //FIXME: reduced
352              }                  switch (outFS) {
353              break;                      case Nodes:
354          default:                      case ReducedNodes: //FIXME: reduced
355              throw RipleyException(msg.str());                          if (getMPISize()==1)
356                                copyData(target, *const_cast<escript::Data*>(&in));
357                            else
358                                dofToNodes(target, *const_cast<escript::Data*>(&in));
359                            break;
360    
361                        case DegreesOfFreedom:
362                        case ReducedDegreesOfFreedom: //FIXME: reduced
363                            copyData(target, *const_cast<escript::Data*>(&in));
364                            break;
365    
366                        case Elements:
367                        case ReducedElements:
368                            if (getMPISize()==1) {
369                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
370                            } else {
371                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
372                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
373                            }
374                            break;
375    
376                        case FaceElements:
377                        case ReducedFaceElements:
378                            if (getMPISize()==1) {
379                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
380                            } else {
381                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
382                                interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);
383                            }
384                            break;
385    
386                        default:
387                            throw RipleyException(msg.str());
388                    }
389                    break;
390                default:
391                    throw RipleyException(msg.str());
392            }
393      }      }
394  }  }
395    
# Line 331  bool RipleyDomain::isCellOriented(int fs Line 444  bool RipleyDomain::isCellOriented(int fs
444              break;              break;
445      }      }
446      stringstream msg;      stringstream msg;
447      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
448            << " on " << getDescription();
449      throw RipleyException(msg.str());      throw RipleyException(msg.str());
450  }  }
451    
# Line 352  bool RipleyDomain::canTag(int fsType) co Line 466  bool RipleyDomain::canTag(int fsType) co
466              break;              break;
467      }      }
468      stringstream msg;      stringstream msg;
469      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag(): Illegal function space type " << fsType << " on "
470            << getDescription();
471      throw RipleyException(msg.str());      throw RipleyException(msg.str());
472  }  }
473    
# Line 378  void RipleyDomain::setTags(const int fsT Line 493  void RipleyDomain::setTags(const int fsT
493              break;              break;
494          default: {          default: {
495              stringstream msg;              stringstream msg;
496              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags(): not implemented for "
497                    << functionSpaceTypeAsString(fsType);
498              throw RipleyException(msg.str());              throw RipleyException(msg.str());
499          }          }
500      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 531  int RipleyDomain::getTagFromSampleNo(int
531              break;              break;
532          default: {          default: {
533              stringstream msg;              stringstream msg;
534              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
535                    << functionSpaceTypeAsString(fsType);
536              throw RipleyException(msg.str());              throw RipleyException(msg.str());
537          }          }
538      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 552  int RipleyDomain::getNumberOfTagsInUse(i
552              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
553          default: {          default: {
554              stringstream msg;              stringstream msg;
555              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
556                    << functionSpaceTypeAsString(fsType);
557              throw RipleyException(msg.str());              throw RipleyException(msg.str());
558          }          }
559      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 572  const int* RipleyDomain::borrowListOfTag
572              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
573          default: {          default: {
574              stringstream msg;              stringstream msg;
575              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
576                    << functionSpaceTypeAsString(fsType);
577              throw RipleyException(msg.str());              throw RipleyException(msg.str());
578          }          }
579      }      }
# Line 502  escript::ASM_ptr RipleyDomain::newSystem Line 621  escript::ASM_ptr RipleyDomain::newSystem
621      // is the domain right?      // is the domain right?
622      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
623      if (row_domain!=*this)      if (row_domain!=*this)
624          throw RipleyException("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");
625      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626      if (col_domain!=*this)      if (col_domain!=*this)
627          throw RipleyException("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");
628      // is the function space type right?      // is the function space type right?
629      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630          reduceRowOrder=true;          reduceRowOrder=true;
631      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634          reduceColOrder=true;          reduceColOrder=true;
635      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637    
638      // generate matrix      // generate matrix
639      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 527  escript::ASM_ptr RipleyDomain::newSystem Line 646  escript::ASM_ptr RipleyDomain::newSystem
646      return sma;      return sma;
647  }  }
648    
649    void RipleyDomain::addPDEToSystem(
650            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651            const escript::Data& A, const escript::Data& B, const escript::Data& C,
652            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653            const escript::Data& d, const escript::Data& y,
654            const escript::Data& d_contact, const escript::Data& y_contact,
655            const escript::Data& d_dirac,const escript::Data& y_dirac) const
656    {
657        if (!d_contact.isEmpty() || !y_contact.isEmpty())
658            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
659    
660        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
661        if (!sma)
662            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
663    
664        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
666    
667        int fsType=UNKNOWN;
668        fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
669        fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
670        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
671        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
672        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
673        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
674        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
675        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
676        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
677        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
678        if (fsType==UNKNOWN)
679            return;
680    
681        const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
682    
683        //TODO: more input param checks (shape, etc)
684    
685        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
686    
687        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
688            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
689    
690        const int numEq=S->logical_row_block_size;
691        const int numComp=S->logical_col_block_size;
692    
693        if (numEq != numComp)
694            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
695        //TODO: more system matrix checks
696    
697        if (numEq==1)
698            if (reducedOrder)
699                assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y, d, y);
700            else
701                assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
702        else
703            if (reducedOrder)
704                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y, d, y);
705            else
706                assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
707    }
708    
709    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
710            const escript::Data& Y, const escript::Data& y,
711            const escript::Data& y_contact, const escript::Data& y_dirac) const
712    {
713        if (!y_contact.isEmpty())
714            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
715    
716        if (rhs.isEmpty()) {
717            if (!X.isEmpty() || !Y.isEmpty())
718                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
719            else
720                return;
721        }
722    
723        if (rhs.getDataPointSize() == 1)
724            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
725        else
726            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
727    }
728    
729  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
730  {  {
731      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
732  }  }
733    
734  //protected  //protected
735  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
736  {  {
737      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
738        const dim_t dpp = in.getNumDataPointsPerSample();
739        out.requireWrite();
740    #pragma omp parallel for
741        for (index_t i=0; i<in.getNumSamples(); i++) {
742            const double* src = in.getSampleDataRO(i);
743            double* dest = out.getSampleDataRW(i);
744            for (index_t c=0; c<numComp; c++) {
745                double res=0.;
746                for (index_t q=0; q<dpp; q++)
747                    res+=src[c+q*numComp];
748                *dest++ = res/dpp;
749            }
750        }
751    }
752    
753    //protected
754    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
755    {
756        const dim_t numComp = in.getDataPointSize();
757        out.requireWrite();
758  #pragma omp parallel for  #pragma omp parallel for
759      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
760          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 812  void RipleyDomain::updateTagsInUse(int f
812          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
813          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);          MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
814  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
815    
816            // if we found a new value add it to the tagsInUse vector
817          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
818              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
819              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 822  void RipleyDomain::updateTagsInUse(int f
822      }      }
823  }  }
824    
825    //protected
826    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
827            const IndexVector& index, const dim_t M, const dim_t N) const
828    {
829        // paso will manage the memory
830        index_t* indexC = MEMALLOC(index.size(), index_t);
831        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
832        copy(index.begin(), index.end(), indexC);
833        copy(ptr.begin(), ptr.end(), ptrC);
834        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
835    }
836    
837    //protected
838    Paso_Pattern* RipleyDomain::createMainPattern() const
839    {
840        IndexVector ptr(1,0);
841        IndexVector index;
842    
843        for (index_t i=0; i<getNumDOF(); i++) {
844            // add the DOF itself
845            index.push_back(i);
846            const dim_t num=insertNeighbourNodes(index, i);
847            ptr.push_back(ptr.back()+num+1);
848        }
849    
850        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
851    }
852    
853    //protected
854    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
855            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
856    {
857        IndexVector ptr(1,0);
858        IndexVector index;
859        for (index_t i=0; i<getNumDOF(); i++) {
860            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
861            ptr.push_back(ptr.back()+colIndices[i].size());
862        }
863    
864        const dim_t M=ptr.size()-1;
865        *colPattern=createPasoPattern(ptr, index, M, N);
866    
867        IndexVector rowPtr(1,0);
868        IndexVector rowIndex;
869        for (dim_t id=0; id<N; id++) {
870            dim_t n=0;
871            for (dim_t i=0; i<M; i++) {
872                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
873                    if (index[j]==id) {
874                        rowIndex.push_back(i);
875                        n++;
876                        break;
877                    }
878                }
879            }
880            rowPtr.push_back(rowPtr.back()+n);
881        }
882    
883        // M and N are now reversed
884        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
885    }
886    
887    //protected
888    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
889           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
890           dim_t num_Sol, const vector<double>& array) const
891    {
892        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
893        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
894        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
895        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
896    
897        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
898        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
899        double* mainBlock_val = mat->mainBlock->val;
900        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
901        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
902        double* col_coupleBlock_val = mat->col_coupleBlock->val;
903        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
904        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
905        double* row_coupleBlock_val = mat->row_coupleBlock->val;
906    
907        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
908            // down columns of array
909            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
910                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
911                //printf("row:%d\n", i_row);
912                // only look at the matrix rows stored on this processor
913                if (i_row < numMyRows) {
914                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
915                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
916                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
917                            if (i_col < numMyCols) {
918                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
919                                    if (mainBlock_index[k] == i_col) {
920                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
921                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
922                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
923                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
924                                                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())];
925                                            }
926                                        }
927                                        break;
928                                    }
929                                }
930                            } else {
931                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
932            //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
933                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
934                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
935                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
936                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
937                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
938                                                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())];
939                                            }
940                                        }
941                                        break;
942                                    }
943                                }
944                            }
945                        }
946                    }
947                } else {
948                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
949                        // across rows of array
950                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
951                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
952                            if (i_col < numMyCols) {
953                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
954                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
955                                {
956            //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
957                                    if (row_coupleBlock_index[k] == i_col) {
958                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
959                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
960                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
961                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
962                                            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())];
963                                            }
964                                        }
965                                        break;
966                                    }
967                                }
968                            }
969                        }
970                    }
971                }
972            }
973        }
974    }
975    
976  //  //
977  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
978  //  //
# Line 612  string RipleyDomain::getDescription() co Line 982  string RipleyDomain::getDescription() co
982      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
983  }  }
984    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
985  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
986  {  {
987      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 997  const int* RipleyDomain::borrowSampleRef
997      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
998  }  }
999    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
1000  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1001  {  {
1002      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 674  bool RipleyDomain::ownSample(int fsType, Line 1033  bool RipleyDomain::ownSample(int fsType,
1033      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1034  }  }
1035    
 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");  
 }  
   
1036  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1037          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
1038          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 692  void RipleyDomain::addPDEToLumpedSystem( Line 1040  void RipleyDomain::addPDEToLumpedSystem(
1040      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("addPDEToLumpedSystem() not implemented");
1041  }  }
1042    
 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");  
 }  
   
1043  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1044          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1045          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1110  dim_t RipleyDomain::getNumNodes() const
1110      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1111  }  }
1112    
1113    dim_t RipleyDomain::getNumDOF() const
1114    {
1115        throw RipleyException("getNumDOF() not implemented");
1116    }
1117    
1118    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1119    {
1120        throw RipleyException("insertNeighbourNodes() not implemented");
1121    }
1122    
1123  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1124  {  {
1125      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1126  }  }
1127    
1128    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1129            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1130            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1131            const escript::Data& d, const escript::Data& y) const
1132    {
1133        throw RipleyException("assemblePDESingle() not implemented");
1134    }
1135    
1136    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1137            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1138            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1139            const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1140    {
1141        throw RipleyException("assemblePDESingleReduced() not implemented");
1142    }
1143    
1144    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1145            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1146            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1147            const escript::Data& d, const escript::Data& y) const
1148    {
1149        throw RipleyException("assemblePDESystem() not implemented");
1150    }
1151    
1152    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1153            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1154            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1155            const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1156    {
1157        throw RipleyException("assemblePDESystemReduced() not implemented");
1158    }
1159    
1160  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1161  {  {
1162      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1167  void RipleyDomain::interpolateNodesOnFac
1167      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1168  }  }
1169    
1170    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1171    {
1172        throw RipleyException("nodesToDOF() not implemented");
1173    }
1174    
1175    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1176    {
1177        throw RipleyException("dofToNodes() not implemented");
1178    }
1179    
1180  } // end of namespace ripley  } // end of namespace ripley
1181    

Legend:
Removed from v.3722  
changed lines
  Added in v.3761

  ViewVC Help
Powered by ViewVC 1.1.26