/[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 3777 by caltinay, Thu Jan 19 06:17:38 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            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(): Illegal 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(): Illegal 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(): not implemented for "
567                    << functionSpaceTypeAsString(fsType);
568              throw RipleyException(msg.str());              throw RipleyException(msg.str());
569          }          }
570      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 601  int RipleyDomain::getTagFromSampleNo(int
601              break;              break;
602          default: {          default: {
603              stringstream msg;              stringstream msg;
604              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
605                    << functionSpaceTypeAsString(fsType);
606              throw RipleyException(msg.str());              throw RipleyException(msg.str());
607          }          }
608      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 622  int RipleyDomain::getNumberOfTagsInUse(i
622              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
623          default: {          default: {
624              stringstream msg;              stringstream msg;
625              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
626                    << functionSpaceTypeAsString(fsType);
627              throw RipleyException(msg.str());              throw RipleyException(msg.str());
628          }          }
629      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 642  const int* RipleyDomain::borrowListOfTag
642              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
643          default: {          default: {
644              stringstream msg;              stringstream msg;
645              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
646                    << functionSpaceTypeAsString(fsType);
647              throw RipleyException(msg.str());              throw RipleyException(msg.str());
648          }          }
649      }      }
# Line 501  escript::ASM_ptr RipleyDomain::newSystem Line 690  escript::ASM_ptr RipleyDomain::newSystem
690      bool reduceColOrder=false;      bool reduceColOrder=false;
691      // is the domain right?      // is the domain right?
692      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
693      if (row_domain!=*this)      if (row_domain != *this)
694          throw RipleyException("Domain of row function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
695      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
696      if (col_domain!=*this)      if (col_domain != *this)
697          throw RipleyException("Domain of column function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
698      // is the function space type right?      // is the function space type right?
699      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
700          reduceRowOrder=true;          reduceRowOrder=true;
701      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
702          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
703      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
704          reduceColOrder=true;          reduceColOrder=true;
705      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
706          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
707    
708      // generate matrix      // generate matrix
709      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 527  escript::ASM_ptr RipleyDomain::newSystem Line 716  escript::ASM_ptr RipleyDomain::newSystem
716      return sma;      return sma;
717  }  }
718    
719    void RipleyDomain::addPDEToSystem(
720            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
721            const escript::Data& A, const escript::Data& B, const escript::Data& C,
722            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
723            const escript::Data& d, const escript::Data& y,
724            const escript::Data& d_contact, const escript::Data& y_contact,
725            const escript::Data& d_dirac,const escript::Data& y_dirac) const
726    {
727        if (!d_contact.isEmpty() || !y_contact.isEmpty())
728            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
729    
730        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
731        if (!sma)
732            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
733    
734        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
735            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
736    
737        int fsType=UNKNOWN;
738        fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
739        fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
740        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
741        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
742        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
743        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
744        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
745        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
746        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
747        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
748        if (fsType==UNKNOWN)
749            return;
750    
751        const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
752    
753        //TODO: more input param checks (shape, etc)
754    
755        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
756    
757        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
758            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
759    
760        const int numEq=S->logical_row_block_size;
761        const int numComp=S->logical_col_block_size;
762    
763        if (numEq != numComp)
764            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
765        //TODO: more system matrix checks
766    
767        if (numEq==1) {
768            if (reducedOrder) {
769                assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
770                assemblePDEBoundarySingleReduced(S, rhs, d, y);
771            } else {
772                assemblePDESingle(S, rhs, A, B, C, D, X, Y);
773                assemblePDEBoundarySingle(S, rhs, d, y);
774            }
775        } else {
776            if (reducedOrder) {
777                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
778                assemblePDEBoundarySystemReduced(S, rhs, d, y);
779            } else {
780                assemblePDESystem(S, rhs, A, B, C, D, X, Y);
781                assemblePDEBoundarySystem(S, rhs, d, y);
782            }
783        }
784    }
785    
786    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
787            const escript::Data& Y, const escript::Data& y,
788            const escript::Data& y_contact, const escript::Data& y_dirac) const
789    {
790        if (!y_contact.isEmpty())
791            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
792    
793        if (rhs.isEmpty()) {
794            if (!X.isEmpty() || !Y.isEmpty())
795                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
796            else
797                return;
798        }
799    
800        if (rhs.getDataPointSize() == 1) {
801            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
802            assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
803        } else {
804            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
805            assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
806        }
807    }
808    
809  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
810  {  {
811      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
812  }  }
813    
814  //protected  //protected
815  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
816    {
817        const dim_t numComp = in.getDataPointSize();
818        const dim_t dpp = in.getNumDataPointsPerSample();
819        out.requireWrite();
820    #pragma omp parallel for
821        for (index_t i=0; i<in.getNumSamples(); i++) {
822            const double* src = in.getSampleDataRO(i);
823            double* dest = out.getSampleDataRW(i);
824            for (index_t c=0; c<numComp; c++) {
825                double res=0.;
826                for (index_t q=0; q<dpp; q++)
827                    res+=src[c+q*numComp];
828                *dest++ = res/dpp;
829            }
830        }
831    }
832    
833    //protected
834    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
835  {  {
836      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
837        out.requireWrite();
838  #pragma omp parallel for  #pragma omp parallel for
839      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
840          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 892  void RipleyDomain::updateTagsInUse(int f
892          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
893          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);
894  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
895    
896            // if we found a new value add it to the tagsInUse vector
897          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
898              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
899              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 902  void RipleyDomain::updateTagsInUse(int f
902      }      }
903  }  }
904    
905    //protected
906    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
907            const IndexVector& index, const dim_t M, const dim_t N) const
908    {
909        // paso will manage the memory
910        index_t* indexC = MEMALLOC(index.size(), index_t);
911        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
912        copy(index.begin(), index.end(), indexC);
913        copy(ptr.begin(), ptr.end(), ptrC);
914        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
915    }
916    
917    //protected
918    Paso_Pattern* RipleyDomain::createMainPattern() const
919    {
920        IndexVector ptr(1,0);
921        IndexVector index;
922    
923        for (index_t i=0; i<getNumDOF(); i++) {
924            // add the DOF itself
925            index.push_back(i);
926            const dim_t num=insertNeighbourNodes(index, i);
927            ptr.push_back(ptr.back()+num+1);
928        }
929    
930        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
931    }
932    
933    //protected
934    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
935            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
936    {
937        IndexVector ptr(1,0);
938        IndexVector index;
939        for (index_t i=0; i<getNumDOF(); i++) {
940            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
941            ptr.push_back(ptr.back()+colIndices[i].size());
942        }
943    
944        const dim_t M=ptr.size()-1;
945        *colPattern=createPasoPattern(ptr, index, M, N);
946    
947        IndexVector rowPtr(1,0);
948        IndexVector rowIndex;
949        for (dim_t id=0; id<N; id++) {
950            dim_t n=0;
951            for (dim_t i=0; i<M; i++) {
952                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
953                    if (index[j]==id) {
954                        rowIndex.push_back(i);
955                        n++;
956                        break;
957                    }
958                }
959            }
960            rowPtr.push_back(rowPtr.back()+n);
961        }
962    
963        // M and N are now reversed
964        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
965    }
966    
967    //protected
968    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
969           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
970           dim_t num_Sol, const vector<double>& array) const
971    {
972        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
973        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
974        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
975        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
976    
977        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
978        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
979        double* mainBlock_val = mat->mainBlock->val;
980        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
981        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
982        double* col_coupleBlock_val = mat->col_coupleBlock->val;
983        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
984        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
985        double* row_coupleBlock_val = mat->row_coupleBlock->val;
986    
987        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
988            // down columns of array
989            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
990                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
991                // only look at the matrix rows stored on this processor
992                if (i_row < numMyRows) {
993                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
994                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
995                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
996                            if (i_col < numMyCols) {
997                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
998                                    if (mainBlock_index[k] == i_col) {
999                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1000                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1001                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1002                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1003                                                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())];
1004                                            }
1005                                        }
1006                                        break;
1007                                    }
1008                                }
1009                            } else {
1010                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1011                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1012                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1013                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1014                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1015                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1016                                                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())];
1017                                            }
1018                                        }
1019                                        break;
1020                                    }
1021                                }
1022                            }
1023                        }
1024                    }
1025                } else {
1026                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1027                        // across rows of array
1028                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1029                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1030                            if (i_col < numMyCols) {
1031                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1032                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1033                                {
1034                                    if (row_coupleBlock_index[k] == i_col) {
1035                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1036                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1037                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1038                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1039                                            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())];
1040                                            }
1041                                        }
1042                                        break;
1043                                    }
1044                                }
1045                            }
1046                        }
1047                    }
1048                }
1049            }
1050        }
1051    }
1052    
1053  //  //
1054  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
1055  //  //
# Line 612  string RipleyDomain::getDescription() co Line 1059  string RipleyDomain::getDescription() co
1059      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
1060  }  }
1061    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
1062  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
1063  {  {
1064      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 1074  const int* RipleyDomain::borrowSampleRef
1074      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1075  }  }
1076    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
1077  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1078  {  {
1079      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 659  void RipleyDomain::setToSize(escript::Da Line 1095  void RipleyDomain::setToSize(escript::Da
1095      throw RipleyException("setToSize() not implemented");      throw RipleyException("setToSize() not implemented");
1096  }  }
1097    
 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");  
 }  
   
1098  bool RipleyDomain::ownSample(int fsType, index_t id) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1099  {  {
1100      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1101  }  }
1102    
 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");  
 }  
   
1103  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1104          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
1105          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 692  void RipleyDomain::addPDEToLumpedSystem( Line 1107  void RipleyDomain::addPDEToLumpedSystem(
1107      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("addPDEToLumpedSystem() not implemented");
1108  }  }
1109    
 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");  
 }  
   
1110  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1111          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1112          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 749  IndexVector RipleyDomain::getNodeDistrib Line 1157  IndexVector RipleyDomain::getNodeDistrib
1157      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1158  }  }
1159    
1160    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1161    {
1162        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1163    }
1164    
1165  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1166  {  {
1167      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1182  dim_t RipleyDomain::getNumNodes() const
1182      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1183  }  }
1184    
1185    dim_t RipleyDomain::getNumDOF() const
1186    {
1187        throw RipleyException("getNumDOF() not implemented");
1188    }
1189    
1190    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1191    {
1192        throw RipleyException("insertNeighbourNodes() not implemented");
1193    }
1194    
1195  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1196  {  {
1197      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1198  }  }
1199    
1200    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1201    {
1202        throw RipleyException("assembleGradient() not implemented");
1203    }
1204    
1205    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1206    {
1207        throw RipleyException("assembleIntegrate() not implemented");
1208    }
1209    
1210    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1211            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1212            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1213    {
1214        throw RipleyException("assemblePDESingle() not implemented");
1215    }
1216    
1217    void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1218          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1219    {
1220        throw RipleyException("assemblePDEBoundarySingle() not implemented");
1221    }
1222    
1223    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1224            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1225            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1226            const escript::Data& Y) const
1227    {
1228        throw RipleyException("assemblePDESingleReduced() not implemented");
1229    }
1230    
1231    void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1232          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1233    {
1234        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1235    }
1236    
1237    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1238            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1239            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1240    {
1241        throw RipleyException("assemblePDESystem() not implemented");
1242    }
1243    
1244    void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1245          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1246    {
1247        throw RipleyException("assemblePDEBoundarySystem() not implemented");
1248    }
1249    
1250    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1251            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1252            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1253            const escript::Data& Y) const
1254    {
1255        throw RipleyException("assemblePDESystemReduced() not implemented");
1256    }
1257    
1258    void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1259          escript::Data& rhs, const escript::Data& d, const escript::Data& y) const
1260    {
1261        throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1262    }
1263    
1264  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1265  {  {
1266      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1271  void RipleyDomain::interpolateNodesOnFac
1271      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1272  }  }
1273    
1274    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1275    {
1276        throw RipleyException("nodesToDOF() not implemented");
1277    }
1278    
1279    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1280    {
1281        throw RipleyException("dofToNodes() not implemented");
1282    }
1283    
1284  } // end of namespace ripley  } // end of namespace ripley
1285    

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

  ViewVC Help
Powered by ViewVC 1.1.26