/[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

branches/ripleygmg_from_3668/ripley/src/RipleyDomain.cpp revision 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC trunk/ripley/src/RipleyDomain.cpp revision 3793 by gross, Wed Feb 1 07:39:43 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: Invalid 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 function space 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 function space 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            if (in.actsExpanded())
301                averageData(target, *const_cast<escript::Data*>(&in));
302            else
303                copyData(target, *const_cast<escript::Data*>(&in));
304        } else {
305            switch (inFS) {
306                case Nodes:
307                case ReducedNodes: //FIXME: reduced
308                    switch (outFS) {
309                        case DegreesOfFreedom:
310                        case ReducedDegreesOfFreedom: //FIXME: reduced
311                            if (getMPISize()==1)
312                                copyData(target, *const_cast<escript::Data*>(&in));
313                            else
314                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
315                            break;
316    
317                        case Nodes:
318                        case ReducedNodes: //FIXME: reduced
319                            copyData(target, *const_cast<escript::Data*>(&in));
320                            break;
321    
322                        case Elements:
323                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
324                            break;
325    
326                        case ReducedElements:
327                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
328                            break;
329    
330                        case FaceElements:
331                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
332                            break;
333    
334                        case ReducedFaceElements:
335                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
336                            break;
337    
338                        default:
339                            throw RipleyException(msg.str());
340                    }
341                    break;
342    
343                  default:              case DegreesOfFreedom:
344                      throw RipleyException(msg.str());              case ReducedDegreesOfFreedom: //FIXME: reduced
345              }                  switch (outFS) {
346              break;                      case Nodes:
347          default:                      case ReducedNodes: //FIXME: reduced
348              throw RipleyException(msg.str());                          if (getMPISize()==1)
349                                copyData(target, *const_cast<escript::Data*>(&in));
350                            else
351                                dofToNodes(target, *const_cast<escript::Data*>(&in));
352                            break;
353    
354                        case DegreesOfFreedom:
355                        case ReducedDegreesOfFreedom: //FIXME: reduced
356                            copyData(target, *const_cast<escript::Data*>(&in));
357                            break;
358    
359                        case Elements:
360                        case ReducedElements:
361                            if (getMPISize()==1) {
362                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
363                            } else {
364                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
365                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
366                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
367                            }
368                            break;
369    
370                        case FaceElements:
371                        case ReducedFaceElements:
372                            if (getMPISize()==1) {
373                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
374                            } else {
375                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
376                                         escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
377                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
378                            }
379                            break;
380    
381                        default:
382                            throw RipleyException(msg.str());
383                    }
384                    break;
385                default:
386                    throw RipleyException(msg.str());
387            }
388      }      }
389  }  }
390    
# Line 314  void RipleyDomain::setToX(escript::Data& Line 422  void RipleyDomain::setToX(escript::Data&
422      }      }
423  }  }
424    
425    void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
426    {
427        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
428                *(arg.getFunctionSpace().getDomain()));
429        if (argDomain != *this)
430            throw RipleyException("setToGradient: Illegal domain of gradient argument");
431        const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
432                *(grad.getFunctionSpace().getDomain()));
433        if (gradDomain != *this)
434            throw RipleyException("setToGradient: Illegal domain of gradient");
435    
436        switch (grad.getFunctionSpace().getTypeCode()) {
437            case Elements:
438            case ReducedElements:
439            case FaceElements:
440            case ReducedFaceElements:
441                break;
442            default: {
443                stringstream msg;
444                msg << "setToGradient: not supported for "
445                    << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
446                throw RipleyException(msg.str());
447            }
448        }
449    
450        if (getMPISize()>1) {
451            if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
452                escript::Data contArg(arg, escript::continuousFunction(*this));
453                assembleGradient(grad, contArg);
454            } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
455                escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
456                assembleGradient(grad, contArg);
457            } else {
458                assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459            }
460        } else {
461            assembleGradient(grad, *const_cast<escript::Data*>(&arg));
462        }
463    }
464    
465    void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
466    {
467        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
468                *(arg.getFunctionSpace().getDomain()));
469        if (argDomain != *this)
470            throw RipleyException("setToIntegrals: illegal domain of integration kernel");
471    
472        switch (arg.getFunctionSpace().getTypeCode()) {
473            case Nodes:
474            case ReducedNodes:
475            case DegreesOfFreedom:
476            case ReducedDegreesOfFreedom:
477                {
478                    escript::Data funcArg(arg, escript::function(*this));
479                    assembleIntegrate(integrals, funcArg);
480                }
481                break;
482            case Elements:
483            case ReducedElements:
484            case FaceElements:
485            case ReducedFaceElements:
486                assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
487                break;
488            default: {
489                stringstream msg;
490                msg << "setToIntegrals: not supported for "
491                    << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
492                throw RipleyException(msg.str());
493            }
494        }
495    
496    }
497    
498  bool RipleyDomain::isCellOriented(int fsType) const  bool RipleyDomain::isCellOriented(int fsType) const
499  {  {
500      switch(fsType) {      switch(fsType) {
501          case Nodes:          case Nodes:
502            case ReducedNodes:
503          case DegreesOfFreedom:          case DegreesOfFreedom:
504          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
505              return false;              return false;
# Line 331  bool RipleyDomain::isCellOriented(int fs Line 513  bool RipleyDomain::isCellOriented(int fs
513              break;              break;
514      }      }
515      stringstream msg;      stringstream msg;
516      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented: invalid function space type " << fsType
517            << " on " << getDescription();
518      throw RipleyException(msg.str());      throw RipleyException(msg.str());
519  }  }
520    
# Line 347  bool RipleyDomain::canTag(int fsType) co Line 530  bool RipleyDomain::canTag(int fsType) co
530          case DegreesOfFreedom:          case DegreesOfFreedom:
531          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
532          case Points:          case Points:
533            case ReducedNodes:
534              return false;              return false;
535          default:          default:
536              break;              break;
537      }      }
538      stringstream msg;      stringstream msg;
539      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag: invalid function space type " << fsType << " on "
540            << getDescription();
541      throw RipleyException(msg.str());      throw RipleyException(msg.str());
542  }  }
543    
# Line 378  void RipleyDomain::setTags(const int fsT Line 563  void RipleyDomain::setTags(const int fsT
563              break;              break;
564          default: {          default: {
565              stringstream msg;              stringstream msg;
566              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags: invalid function space type " << fsType;
567              throw RipleyException(msg.str());              throw RipleyException(msg.str());
568          }          }
569      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 600  int RipleyDomain::getTagFromSampleNo(int
600              break;              break;
601          default: {          default: {
602              stringstream msg;              stringstream msg;
603              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo: invalid function space type " << fsType;
604              throw RipleyException(msg.str());              throw RipleyException(msg.str());
605          }          }
606      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 620  int RipleyDomain::getNumberOfTagsInUse(i
620              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
621          default: {          default: {
622              stringstream msg;              stringstream msg;
623              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse: invalid function space type "
624                    << fsType;
625              throw RipleyException(msg.str());              throw RipleyException(msg.str());
626          }          }
627      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 640  const int* RipleyDomain::borrowListOfTag
640              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
641          default: {          default: {
642              stringstream msg;              stringstream msg;
643              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse: invalid function space type "
644                    << fsType;
645              throw RipleyException(msg.str());              throw RipleyException(msg.str());
646          }          }
647      }      }
# Line 501  escript::ASM_ptr RipleyDomain::newSystem Line 688  escript::ASM_ptr RipleyDomain::newSystem
688      bool reduceColOrder=false;      bool reduceColOrder=false;
689      // is the domain right?      // is the domain right?
690      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
691      if (row_domain!=*this)      if (row_domain != *this)
692          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");
693      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
694      if (col_domain!=*this)      if (col_domain != *this)
695          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");
696      // is the function space type right?      // is the function space type right?
697      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
698          reduceRowOrder=true;          reduceRowOrder=true;
699      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
700          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix: illegal function space type for system matrix rows");
701      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
702          reduceColOrder=true;          reduceColOrder=true;
703      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
704          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix: illegal function space type for system matrix columns");
705    
706      // generate matrix      // generate matrix
707      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
708      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,      Paso_SystemMatrix* matrix = Paso_SystemMatrix_alloc(type, pattern,
709              row_blocksize, column_blocksize, FALSE);              row_blocksize, column_blocksize, FALSE);
710      paso::checkPasoError();      paso::checkPasoError();
     Paso_SystemMatrixPattern_free(pattern);  
711      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,      escript::ASM_ptr sma(new SystemMatrixAdapter(matrix, row_blocksize,
712                  row_functionspace, column_blocksize, column_functionspace));                  row_functionspace, column_blocksize, column_functionspace));
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
812    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  //protected
831  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  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");  
 }  
   
 void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  
         const escript::Data& D, const escript::Data& d,  
         const escript::Data& d_dirac, const bool useHRZ) const  
 {  
     throw RipleyException("addPDEToLumpedSystem() not implemented");  
 }  
   
 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,  
         const escript::Data& Y, const escript::Data& y,  
         const escript::Data& y_contact, const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToRHS() not implemented");  
 }  
   
1100  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1101          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1102          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 711  void RipleyDomain::addPDEToTransportProb Line 1109  void RipleyDomain::addPDEToTransportProb
1109      throw RipleyException("addPDEToTransportProblem() not implemented");      throw RipleyException("addPDEToTransportProblem() not implemented");
1110  }  }
1111    
1112  escript::ATP_ptr RipleyDomain::newTransportProblem(const bool useBackwardEuler,  escript::ATP_ptr RipleyDomain::newTransportProblem(
1113          const int blocksize, const escript::FunctionSpace& functionspace,          const int blocksize, const escript::FunctionSpace& functionspace,
1114          const int type) const          const int type) const
1115  {  {
# Line 749  IndexVector RipleyDomain::getNodeDistrib Line 1147  IndexVector RipleyDomain::getNodeDistrib
1147      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1148  }  }
1149    
1150    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1151    {
1152        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1153    }
1154    
1155  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1156  {  {
1157      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1172  dim_t RipleyDomain::getNumNodes() const
1172      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1173  }  }
1174    
1175    dim_t RipleyDomain::getNumDOF() const
1176    {
1177        throw RipleyException("getNumDOF() not implemented");
1178    }
1179    
1180    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1181    {
1182        throw RipleyException("insertNeighbourNodes() not implemented");
1183    }
1184    
1185  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1186  {  {
1187      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1188  }  }
1189    
1190    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1191    {
1192        throw RipleyException("assembleGradient() not implemented");
1193    }
1194    
1195    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1196    {
1197        throw RipleyException("assembleIntegrate() not implemented");
1198    }
1199    
1200    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1201            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1202            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1203    {
1204        throw RipleyException("assemblePDESingle() not implemented");
1205    }
1206    
1207    void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1208          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1209    {
1210        throw RipleyException("assemblePDEBoundarySingle() not implemented");
1211    }
1212    
1213    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1214            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1215            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1216            const escript::Data& Y) const
1217    {
1218        throw RipleyException("assemblePDESingleReduced() not implemented");
1219    }
1220    
1221    void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1222          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1223    {
1224        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1225    }
1226    
1227    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1228            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1229            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1230    {
1231        throw RipleyException("assemblePDESystem() not implemented");
1232    }
1233    
1234    void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1235          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1236    {
1237        throw RipleyException("assemblePDEBoundarySystem() not implemented");
1238    }
1239    
1240    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1241            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1242            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1243            const escript::Data& Y) const
1244    {
1245        throw RipleyException("assemblePDESystemReduced() not implemented");
1246    }
1247    
1248    void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1249          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1250    {
1251        throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1252    }
1253    
1254  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1255  {  {
1256      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1261  void RipleyDomain::interpolateNodesOnFac
1261      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1262  }  }
1263    
1264    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1265    {
1266        throw RipleyException("nodesToDOF() not implemented");
1267    }
1268    
1269    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1270    {
1271        throw RipleyException("dofToNodes() not implemented");
1272    }
1273    
1274  } // end of namespace ripley  } // end of namespace ripley
1275    

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

  ViewVC Help
Powered by ViewVC 1.1.26