/[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 3776 by caltinay, Thu Jan 19 03:48:35 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 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 72  string RipleyDomain::functionSpaceTypeAs Line 83  string RipleyDomain::functionSpaceTypeAs
83          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85          case Nodes: return "Ripley_Nodes";          case Nodes: return "Ripley_Nodes";
86          case ReducedNodes: return "Ripley_Reduced_Nodes";          case ReducedNodes: return "Ripley_ReducedNodes";
87          case Elements: return "Ripley_Elements";          case Elements: return "Ripley_Elements";
88          case ReducedElements: return "Ripley_Reduced_Elements";          case ReducedElements: return "Ripley_ReducedElements";
89          case FaceElements: return "Ripley_Face_Elements";          case FaceElements: return "Ripley_FaceElements";
90          case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";          case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91          case Points: return "Ripley_Points";          case Points: return "Ripley_Points";
92          default:          default:
93              break;              break;
# Line 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 (!isValidFunctionSpaceType(fsType_target)) {
238            stringstream msg;
239            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
240                << fsType_target << " for " << getDescription();
241            throw RipleyException(msg.str());
242        }
243    
244        switch (fsType_source) {
245            case Nodes:
246            case DegreesOfFreedom:
247                return true;
248            case ReducedNodes:
249            case ReducedDegreesOfFreedom:
250                return (fsType_target != Nodes &&
251                        fsType_target != DegreesOfFreedom);
252            case Elements:
253                return (fsType_target==Elements ||
254                        fsType_target==ReducedElements);
255            case FaceElements:
256                return (fsType_target==FaceElements ||
257                        fsType_target==ReducedFaceElements);
258            case ReducedElements:
259            case ReducedFaceElements:
260            case Points:
261                return (fsType_target==fsType_source);
262    
263            default: {
264                stringstream msg;
265                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
266                    << fsType_source << " for " << getDescription();
267                throw RipleyException(msg.str());
268            }
269        }
270    }
271    
272  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
273                                         const escript::Data& in) const                                         const escript::Data& in) const
274  {  {
# Line 241  void RipleyDomain::interpolateOnDomain(e Line 285  void RipleyDomain::interpolateOnDomain(e
285          << " -> "          << " -> "
286          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
287    
288      switch (in.getFunctionSpace().getTypeCode()) {      const int inFS = in.getFunctionSpace().getTypeCode();
289          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;  
   
                 case FaceElements:  
                     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);  
                     break;  
290    
291                  case ReducedFaceElements:      // simplest case: 1:1 copy
292                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);      if (inFS==outFS) {
293                      break;          copyData(target, *const_cast<escript::Data*>(&in));
294        // not allowed: reduced->non-reduced
295        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
296                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
297            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
298        } else if ((inFS==Elements && outFS==ReducedElements)
299                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
300            averageData(target, *const_cast<escript::Data*>(&in));
301        } else {
302            switch (inFS) {
303                case Nodes:
304                case ReducedNodes: //FIXME: reduced
305                    switch (outFS) {
306                        case DegreesOfFreedom:
307                        case ReducedDegreesOfFreedom: //FIXME: reduced
308                            if (getMPISize()==1)
309                                copyData(target, *const_cast<escript::Data*>(&in));
310                            else
311                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
312                            break;
313    
314                        case Nodes:
315                        case ReducedNodes: //FIXME: reduced
316                            copyData(target, *const_cast<escript::Data*>(&in));
317                            break;
318    
319                        case Elements:
320                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
321                            break;
322    
323                        case ReducedElements:
324                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
325                            break;
326    
327                        case FaceElements:
328                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
329                            break;
330    
331                        case ReducedFaceElements:
332                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
333                            break;
334    
335                        default:
336                            throw RipleyException(msg.str());
337                    }
338                    break;
339    
340                  default:              case DegreesOfFreedom:
341                      throw RipleyException(msg.str());              case ReducedDegreesOfFreedom: //FIXME: reduced
342              }                  switch (outFS) {
343              break;                      case Nodes:
344          default:                      case ReducedNodes: //FIXME: reduced
345              throw RipleyException(msg.str());                          if (getMPISize()==1)
346                                copyData(target, *const_cast<escript::Data*>(&in));
347                            else
348                                dofToNodes(target, *const_cast<escript::Data*>(&in));
349                            break;
350    
351                        case DegreesOfFreedom:
352                        case ReducedDegreesOfFreedom: //FIXME: reduced
353                            copyData(target, *const_cast<escript::Data*>(&in));
354                            break;
355    
356                        case Elements:
357                        case ReducedElements:
358                            if (getMPISize()==1) {
359                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
360                            } else {
361                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
362                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
363                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
364                            }
365                            break;
366    
367                        case FaceElements:
368                        case ReducedFaceElements:
369                            if (getMPISize()==1) {
370                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
371                            } else {
372                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
373                                         escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
374                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
375                            }
376                            break;
377    
378                        default:
379                            throw RipleyException(msg.str());
380                    }
381                    break;
382                default:
383                    throw RipleyException(msg.str());
384            }
385      }      }
386  }  }
387    
# Line 314  void RipleyDomain::setToX(escript::Data& Line 419  void RipleyDomain::setToX(escript::Data&
419      }      }
420  }  }
421    
422    void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
423    {
424        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
425                *(arg.getFunctionSpace().getDomain()));
426        if (argDomain != *this)
427            throw RipleyException("setToGradient: Illegal domain of gradient argument");
428        const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
429                *(grad.getFunctionSpace().getDomain()));
430        if (gradDomain != *this)
431            throw RipleyException("setToGradient: Illegal domain of gradient");
432    
433        switch (grad.getFunctionSpace().getTypeCode()) {
434            case Elements:
435            case ReducedElements:
436            case FaceElements:
437            case ReducedFaceElements:
438                break;
439            default: {
440                stringstream msg;
441                msg << "setToGradient(): not supported for "
442                    << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
443                throw RipleyException(msg.str());
444            }
445        }
446    
447        if (getMPISize()>1) {
448            if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
449                escript::Data contArg(arg, escript::continuousFunction(*this));
450                assembleGradient(grad, contArg);
451            } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
452                escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
453                assembleGradient(grad, contArg);
454            } else {
455                assembleGradient(grad, *const_cast<escript::Data*>(&arg));
456            }
457        } else {
458            assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459        }
460    }
461    
462    void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
463    {
464        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
465                *(arg.getFunctionSpace().getDomain()));
466        if (argDomain != *this)
467            throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
468    
469        switch (arg.getFunctionSpace().getTypeCode()) {
470            case Nodes:
471            case ReducedNodes:
472            case DegreesOfFreedom:
473            case ReducedDegreesOfFreedom:
474                {
475                    escript::Data funcArg(arg, escript::function(*this));
476                    assembleIntegrate(integrals, funcArg);
477                }
478                break;
479            case Elements:
480            case ReducedElements:
481            case FaceElements:
482            case ReducedFaceElements:
483                assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
484                break;
485            default: {
486                stringstream msg;
487                msg << "setToIntegrals(): not supported for "
488                    << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
489                throw RipleyException(msg.str());
490            }
491        }
492    
493    }
494    
495  bool RipleyDomain::isCellOriented(int fsType) const  bool RipleyDomain::isCellOriented(int fsType) const
496  {  {
497      switch(fsType) {      switch(fsType) {
498          case Nodes:          case Nodes:
499            case ReducedNodes:
500          case DegreesOfFreedom:          case DegreesOfFreedom:
501          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
502              return false;              return false;
# Line 331  bool RipleyDomain::isCellOriented(int fs Line 510  bool RipleyDomain::isCellOriented(int fs
510              break;              break;
511      }      }
512      stringstream msg;      stringstream msg;
513      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
514            << " on " << getDescription();
515      throw RipleyException(msg.str());      throw RipleyException(msg.str());
516  }  }
517    
# Line 347  bool RipleyDomain::canTag(int fsType) co Line 527  bool RipleyDomain::canTag(int fsType) co
527          case DegreesOfFreedom:          case DegreesOfFreedom:
528          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
529          case Points:          case Points:
530            case ReducedNodes:
531              return false;              return false;
532          default:          default:
533              break;              break;
534      }      }
535      stringstream msg;      stringstream msg;
536      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag(): Illegal function space type " << fsType << " on "
537            << getDescription();
538      throw RipleyException(msg.str());      throw RipleyException(msg.str());
539  }  }
540    
# Line 378  void RipleyDomain::setTags(const int fsT Line 560  void RipleyDomain::setTags(const int fsT
560              break;              break;
561          default: {          default: {
562              stringstream msg;              stringstream msg;
563              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags(): not implemented for "
564                    << functionSpaceTypeAsString(fsType);
565              throw RipleyException(msg.str());              throw RipleyException(msg.str());
566          }          }
567      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 598  int RipleyDomain::getTagFromSampleNo(int
598              break;              break;
599          default: {          default: {
600              stringstream msg;              stringstream msg;
601              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
602                    << functionSpaceTypeAsString(fsType);
603              throw RipleyException(msg.str());              throw RipleyException(msg.str());
604          }          }
605      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 619  int RipleyDomain::getNumberOfTagsInUse(i
619              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
620          default: {          default: {
621              stringstream msg;              stringstream msg;
622              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
623                    << functionSpaceTypeAsString(fsType);
624              throw RipleyException(msg.str());              throw RipleyException(msg.str());
625          }          }
626      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 639  const int* RipleyDomain::borrowListOfTag
639              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
640          default: {          default: {
641              stringstream msg;              stringstream msg;
642              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
643                    << functionSpaceTypeAsString(fsType);
644              throw RipleyException(msg.str());              throw RipleyException(msg.str());
645          }          }
646      }      }
# Line 501  escript::ASM_ptr RipleyDomain::newSystem Line 687  escript::ASM_ptr RipleyDomain::newSystem
687      bool reduceColOrder=false;      bool reduceColOrder=false;
688      // is the domain right?      // is the domain right?
689      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
690      if (row_domain!=*this)      if (row_domain != *this)
691          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");
692      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
693      if (col_domain!=*this)      if (col_domain != *this)
694          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");
695      // is the function space type right?      // is the function space type right?
696      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
697          reduceRowOrder=true;          reduceRowOrder=true;
698      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
699          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
700      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
701          reduceColOrder=true;          reduceColOrder=true;
702      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
703          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
704    
705      // generate matrix      // generate matrix
706      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 527  escript::ASM_ptr RipleyDomain::newSystem Line 713  escript::ASM_ptr RipleyDomain::newSystem
713      return sma;      return sma;
714  }  }
715    
716    void RipleyDomain::addPDEToSystem(
717            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
718            const escript::Data& A, const escript::Data& B, const escript::Data& C,
719            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
720            const escript::Data& d, const escript::Data& y,
721            const escript::Data& d_contact, const escript::Data& y_contact,
722            const escript::Data& d_dirac,const escript::Data& y_dirac) const
723    {
724        if (!d_contact.isEmpty() || !y_contact.isEmpty())
725            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
726    
727        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
728        if (!sma)
729            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
730    
731        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
732            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
733    
734        int fsType=UNKNOWN;
735        fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
736        fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
737        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
738        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
739        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
740        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
741        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
742        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
743        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
744        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
745        if (fsType==UNKNOWN)
746            return;
747    
748        const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
749    
750        //TODO: more input param checks (shape, etc)
751    
752        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
753    
754        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
755            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
756    
757        const int numEq=S->logical_row_block_size;
758        const int numComp=S->logical_col_block_size;
759    
760        if (numEq != numComp)
761            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
762        //TODO: more system matrix checks
763    
764        if (numEq==1) {
765            if (reducedOrder) {
766                assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
767                assemblePDEBoundarySingleReduced(S, rhs, d, y);
768            } else {
769                assemblePDESingle(S, rhs, A, B, C, D, X, Y);
770                assemblePDEBoundarySingle(S, rhs, d, y);
771            }
772        } else {
773            if (reducedOrder) {
774                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
775                assemblePDEBoundarySystemReduced(S, rhs, d, y);
776            } else {
777                assemblePDESystem(S, rhs, A, B, C, D, X, Y);
778                assemblePDEBoundarySystem(S, rhs, d, y);
779            }
780        }
781    }
782    
783    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
784            const escript::Data& Y, const escript::Data& y,
785            const escript::Data& y_contact, const escript::Data& y_dirac) const
786    {
787        if (!y_contact.isEmpty())
788            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
789    
790        if (rhs.isEmpty()) {
791            if (!X.isEmpty() || !Y.isEmpty())
792                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
793            else
794                return;
795        }
796    
797        if (rhs.getDataPointSize() == 1) {
798            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
799            assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
800        } else {
801            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
802            assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
803        }
804    }
805    
806  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
807  {  {
808      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
809  }  }
810    
811  //protected  //protected
812  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
813    {
814        const dim_t numComp = in.getDataPointSize();
815        const dim_t dpp = in.getNumDataPointsPerSample();
816        out.requireWrite();
817    #pragma omp parallel for
818        for (index_t i=0; i<in.getNumSamples(); i++) {
819            const double* src = in.getSampleDataRO(i);
820            double* dest = out.getSampleDataRW(i);
821            for (index_t c=0; c<numComp; c++) {
822                double res=0.;
823                for (index_t q=0; q<dpp; q++)
824                    res+=src[c+q*numComp];
825                *dest++ = res/dpp;
826            }
827        }
828    }
829    
830    //protected
831    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
832  {  {
833      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
834        out.requireWrite();
835  #pragma omp parallel for  #pragma omp parallel for
836      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
837          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 889  void RipleyDomain::updateTagsInUse(int f
889          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
890          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);
891  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
892    
893            // if we found a new value add it to the tagsInUse vector
894          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
895              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
896              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 899  void RipleyDomain::updateTagsInUse(int f
899      }      }
900  }  }
901    
902    //protected
903    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
904            const IndexVector& index, const dim_t M, const dim_t N) const
905    {
906        // paso will manage the memory
907        index_t* indexC = MEMALLOC(index.size(), index_t);
908        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
909        copy(index.begin(), index.end(), indexC);
910        copy(ptr.begin(), ptr.end(), ptrC);
911        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
912    }
913    
914    //protected
915    Paso_Pattern* RipleyDomain::createMainPattern() const
916    {
917        IndexVector ptr(1,0);
918        IndexVector index;
919    
920        for (index_t i=0; i<getNumDOF(); i++) {
921            // add the DOF itself
922            index.push_back(i);
923            const dim_t num=insertNeighbourNodes(index, i);
924            ptr.push_back(ptr.back()+num+1);
925        }
926    
927        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
928    }
929    
930    //protected
931    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
932            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
933    {
934        IndexVector ptr(1,0);
935        IndexVector index;
936        for (index_t i=0; i<getNumDOF(); i++) {
937            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
938            ptr.push_back(ptr.back()+colIndices[i].size());
939        }
940    
941        const dim_t M=ptr.size()-1;
942        *colPattern=createPasoPattern(ptr, index, M, N);
943    
944        IndexVector rowPtr(1,0);
945        IndexVector rowIndex;
946        for (dim_t id=0; id<N; id++) {
947            dim_t n=0;
948            for (dim_t i=0; i<M; i++) {
949                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
950                    if (index[j]==id) {
951                        rowIndex.push_back(i);
952                        n++;
953                        break;
954                    }
955                }
956            }
957            rowPtr.push_back(rowPtr.back()+n);
958        }
959    
960        // M and N are now reversed
961        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
962    }
963    
964    //protected
965    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
966           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
967           dim_t num_Sol, const vector<double>& array) const
968    {
969        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
970        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
971        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
972        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
973    
974        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
975        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
976        double* mainBlock_val = mat->mainBlock->val;
977        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
978        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
979        double* col_coupleBlock_val = mat->col_coupleBlock->val;
980        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
981        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
982        double* row_coupleBlock_val = mat->row_coupleBlock->val;
983    
984        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
985            // down columns of array
986            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
987                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
988                // only look at the matrix rows stored on this processor
989                if (i_row < numMyRows) {
990                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
991                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
992                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
993                            if (i_col < numMyCols) {
994                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
995                                    if (mainBlock_index[k] == i_col) {
996                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
997                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
998                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
999                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1000                                                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())];
1001                                            }
1002                                        }
1003                                        break;
1004                                    }
1005                                }
1006                            } else {
1007                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1008                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1009                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1010                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1011                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1012                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1013                                                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())];
1014                                            }
1015                                        }
1016                                        break;
1017                                    }
1018                                }
1019                            }
1020                        }
1021                    }
1022                } else {
1023                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1024                        // across rows of array
1025                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1026                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1027                            if (i_col < numMyCols) {
1028                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1029                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1030                                {
1031                                    if (row_coupleBlock_index[k] == i_col) {
1032                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1033                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1034                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1035                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1036                                            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())];
1037                                            }
1038                                        }
1039                                        break;
1040                                    }
1041                                }
1042                            }
1043                        }
1044                    }
1045                }
1046            }
1047        }
1048    }
1049    
1050  //  //
1051  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
1052  //  //
# Line 612  string RipleyDomain::getDescription() co Line 1056  string RipleyDomain::getDescription() co
1056      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
1057  }  }
1058    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
1059  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
1060  {  {
1061      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 1071  const int* RipleyDomain::borrowSampleRef
1071      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1072  }  }
1073    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
1074  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1075  {  {
1076      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 659  void RipleyDomain::setToSize(escript::Da Line 1092  void RipleyDomain::setToSize(escript::Da
1092      throw RipleyException("setToSize() not implemented");      throw RipleyException("setToSize() not implemented");
1093  }  }
1094    
 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");  
 }  
   
1095  bool RipleyDomain::ownSample(int fsType, index_t id) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1096  {  {
1097      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1098  }  }
1099    
 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");  
 }  
   
1100  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1101          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
1102          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 692  void RipleyDomain::addPDEToLumpedSystem( Line 1104  void RipleyDomain::addPDEToLumpedSystem(
1104      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("addPDEToLumpedSystem() not implemented");
1105  }  }
1106    
 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");  
 }  
   
1107  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1108          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1109          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 749  IndexVector RipleyDomain::getNodeDistrib Line 1154  IndexVector RipleyDomain::getNodeDistrib
1154      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1155  }  }
1156    
1157    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1158    {
1159        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1160    }
1161    
1162  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1163  {  {
1164      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1179  dim_t RipleyDomain::getNumNodes() const
1179      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1180  }  }
1181    
1182    dim_t RipleyDomain::getNumDOF() const
1183    {
1184        throw RipleyException("getNumDOF() not implemented");
1185    }
1186    
1187    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1188    {
1189        throw RipleyException("insertNeighbourNodes() not implemented");
1190    }
1191    
1192  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1193  {  {
1194      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1195  }  }
1196    
1197    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1198    {
1199        throw RipleyException("assembleGradient() not implemented");
1200    }
1201    
1202    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1203    {
1204        throw RipleyException("assembleIntegrate() not implemented");
1205    }
1206    
1207    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1208            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1209            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1210    {
1211        throw RipleyException("assemblePDESingle() not implemented");
1212    }
1213    
1214    void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1215          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1216    {
1217        throw RipleyException("assemblePDEBoundarySingle() not implemented");
1218    }
1219    
1220    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1221            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1222            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1223            const escript::Data& Y) const
1224    {
1225        throw RipleyException("assemblePDESingleReduced() not implemented");
1226    }
1227    
1228    void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1229          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1230    {
1231        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1232    }
1233    
1234    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1235            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1236            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1237    {
1238        throw RipleyException("assemblePDESystem() not implemented");
1239    }
1240    
1241    void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1242          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1243    {
1244        throw RipleyException("assemblePDEBoundarySystem() not implemented");
1245    }
1246    
1247    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1248            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1249            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1250            const escript::Data& Y) const
1251    {
1252        throw RipleyException("assemblePDESystemReduced() not implemented");
1253    }
1254    
1255    void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1256          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1257    {
1258        throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1259    }
1260    
1261  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1262  {  {
1263      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1268  void RipleyDomain::interpolateNodesOnFac
1268      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1269  }  }
1270    
1271    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1272    {
1273        throw RipleyException("nodesToDOF() not implemented");
1274    }
1275    
1276    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1277    {
1278        throw RipleyException("dofToNodes() not implemented");
1279    }
1280    
1281  } // end of namespace ripley  } // end of namespace ripley
1282    

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

  ViewVC Help
Powered by ViewVC 1.1.26