/[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 3773 by caltinay, Wed Jan 18 04:27:53 2012 UTC
# Line 1  Line 1 
1    
2  /*******************************************************  /*******************************************************
3  *  *
4  * Copyright (c) 2003-2011 by University of Queensland  * Copyright (c) 2003-2012 by University of Queensland
5  * Earth Systems Science Computational Center (ESSCC)  * Earth Systems Science Computational Center (ESSCC)
6  * http://www.uq.edu.au/esscc  * http://www.uq.edu.au/esscc
7  *  *
# Line 47  RipleyDomain::~RipleyDomain() Line 47  RipleyDomain::~RipleyDomain()
47      Esys_MPIInfo_free(m_mpiInfo);      Esys_MPIInfo_free(m_mpiInfo);
48  }  }
49    
50    bool RipleyDomain::operator==(const AbstractDomain& other) const
51    {
52        const RipleyDomain* o=dynamic_cast<const RipleyDomain*>(&other);
53        if (o) {
54            return (m_tagMap==o->m_tagMap && m_nodeTags==o->m_nodeTags
55                    && m_elementTags==o->m_elementTags
56                    && m_faceTags==o->m_faceTags);
57        }
58        return false;
59    }
60    
61  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const  bool RipleyDomain::isValidFunctionSpaceType(int fsType) const
62  {  {
63      switch (fsType) {      switch (fsType) {
# Line 72  string RipleyDomain::functionSpaceTypeAs Line 83  string RipleyDomain::functionSpaceTypeAs
83          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";          case DegreesOfFreedom: return "Ripley_DegreesOfFreedom";
84          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";          case ReducedDegreesOfFreedom: return "Ripley_ReducedDegreesOfFreedom";
85          case Nodes: return "Ripley_Nodes";          case Nodes: return "Ripley_Nodes";
86          case ReducedNodes: return "Ripley_Reduced_Nodes";          case ReducedNodes: return "Ripley_ReducedNodes";
87          case Elements: return "Ripley_Elements";          case Elements: return "Ripley_Elements";
88          case ReducedElements: return "Ripley_Reduced_Elements";          case ReducedElements: return "Ripley_ReducedElements";
89          case FaceElements: return "Ripley_Face_Elements";          case FaceElements: return "Ripley_FaceElements";
90          case ReducedFaceElements: return "Ripley_Reduced_Face_Elements";          case ReducedFaceElements: return "Ripley_ReducedFaceElements";
91          case Points: return "Ripley_Points";          case Points: return "Ripley_Points";
92          default:          default:
93              break;              break;
# Line 89  pair<int,int> RipleyDomain::getDataShape Line 100  pair<int,int> RipleyDomain::getDataShape
100      const int ptsPerSample = (m_numDim==2 ? 4 : 8);      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
101      switch (fsType) {      switch (fsType) {
102          case Nodes:          case Nodes:
103          case ReducedNodes:          case ReducedNodes: //FIXME: reduced
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
104              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
105            case DegreesOfFreedom:
106            case ReducedDegreesOfFreedom: //FIXME: reduced
107                return pair<int,int>(1, getNumDOF());
108          case Elements:          case Elements:
109              return pair<int,int>(ptsPerSample, getNumElements());              return pair<int,int>(ptsPerSample, getNumElements());
110          case FaceElements:          case FaceElements:
# Line 101  pair<int,int> RipleyDomain::getDataShape Line 113  pair<int,int> RipleyDomain::getDataShape
113              return pair<int,int>(1, getNumElements());              return pair<int,int>(1, getNumElements());
114          case ReducedFaceElements:          case ReducedFaceElements:
115              return pair<int,int>(1, getNumFaceElements());              return pair<int,int>(1, getNumFaceElements());
             /*  
116          case Points:          case Points:
117              return pair<int,int>(1, getNumPoints());              return pair<int,int>(1, 0); //FIXME: dirac
             */  
118          default:          default:
119              break;              break;
120      }      }
121    
122      stringstream msg;      stringstream msg;
123      msg << "getDataShape(): Unsupported function space type "      msg << "getDataShape(): Unsupported function space type " << fsType
124          << functionSpaceTypeAsString(fsType) << " for " << getDescription();          << " for " << getDescription();
125      throw RipleyException(msg.str());      throw RipleyException(msg.str());
126  }  }
127    
# Line 131  bool RipleyDomain::commonFunctionSpace(c Line 141  bool RipleyDomain::commonFunctionSpace(c
141     /*     /*
142      The idea is to use equivalence classes (i.e. types which can be      The idea is to use equivalence classes (i.e. types which can be
143      interpolated back and forth):      interpolated back and forth):
144      class 1: DOF <-> Nodes      class 0: DOF <-> Nodes
145      class 2: ReducedDOF <-> ReducedNodes      class 1: ReducedDOF <-> ReducedNodes
146      class 3: Points      class 2: Points
147      class 4: Elements      class 3: Elements
148      class 5: ReducedElements      class 4: ReducedElements
149      class 6: FaceElements      class 5: FaceElements
150      class 7: ReducedFaceElements      class 6: ReducedFaceElements
     class 8: ContactElementZero <-> ContactElementOne  
     class 9: ReducedContactElementZero <-> ReducedContactElementOne  
151    
152      There is also a set of lines. Interpolation is possible down a line but not      There is also a set of lines. Interpolation is possible down a line but not
153      between lines.      between lines.
154      class 1 and 2 belong to all lines so aren't considered.      class 0 and 1 belong to all lines so aren't considered.
155      line 0: class 3      line 0: class 2
156      line 1: class 4,5      line 1: class 3,4
157      line 2: class 6,7      line 2: class 5,6
     line 3: class 8,9  
158    
159      For classes with multiple members (eg class 2) we have vars to record if      For classes with multiple members (eg class 1) we have vars to record if
160      there is at least one instance. e.g. hasnodes is true if we have at least      there is at least one instance. e.g. hasnodes is true if we have at least
161      one instance of Nodes.      one instance of Nodes.
162      */      */
163      if (fs.empty())      if (fs.empty())
164          return false;          return false;
165      vector<int> hasclass(10);      vector<bool> hasclass(7, false);
166      vector<int> hasline(4);      vector<int> hasline(3, 0);
167      bool hasnodes=false;      bool hasnodes=false;
168      bool hasrednodes=false;      bool hasrednodes=false;
169      for (int i=0;i<fs.size();++i) {      for (size_t i=0; i<fs.size(); ++i) {
170          switch (fs[i]) {          switch (fs[i]) {
171              case Nodes: hasnodes=true; // no break is deliberate              case Nodes: hasnodes=true; // fall through
172              case DegreesOfFreedom:              case DegreesOfFreedom:
173                  hasclass[1]=1;                  hasclass[0]=true;
174                  break;                  break;
175              case ReducedNodes: hasrednodes=true; // no break is deliberate              case ReducedNodes: hasrednodes=true; // fall through
176              case ReducedDegreesOfFreedom:              case ReducedDegreesOfFreedom:
177                  hasclass[2]=1;                  hasclass[1]=true;
178                  break;                  break;
179              case Points:              case Points:
180                  hasline[0]=1;                  hasline[0]=1;
181                  hasclass[3]=1;                  hasclass[2]=true;
182                  break;                  break;
183              case Elements:              case Elements:
184                  hasclass[4]=1;                  hasclass[3]=true;
185                  hasline[1]=1;                  hasline[1]=1;
186                  break;                  break;
187              case ReducedElements:              case ReducedElements:
188                  hasclass[5]=1;                  hasclass[4]=true;
189                  hasline[1]=1;                  hasline[1]=1;
190                  break;                  break;
191              case FaceElements:              case FaceElements:
192                  hasclass[6]=1;                  hasclass[5]=true;
193                  hasline[2]=1;                  hasline[2]=1;
194                  break;                  break;
195              case ReducedFaceElements:              case ReducedFaceElements:
196                  hasclass[7]=1;                  hasclass[6]=true;
197                  hasline[2]=1;                  hasline[2]=1;
198                  break;                  break;
199              default:              default:
200                  return false;                  return false;
201          }          }
202      }      }
203      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int numLines=hasline[0]+hasline[1]+hasline[2];
204    
205      // fail if we have more than one leaf group      // fail if we have more than one leaf group
206      // = there are at least two branches we can't interpolate between      // = there are at least two branches we can't interpolate between
207      if (totlines>1)      if (numLines > 1)
208          return false;          return false;
209      else if (totlines==1) {      else if (numLines==1) {
210          // we have points          // we have points
211          if (hasline[0]==1)          if (hasline[0]==1)
212              resultcode=Points;              resultcode=Points;
213          else if (hasline[1]==1) {          else if (hasline[1]==1) {
214              if (hasclass[5]==1)              if (hasclass[4])
215                  resultcode=ReducedElements;                  resultcode=ReducedElements;
216              else              else
217                  resultcode=Elements;                  resultcode=Elements;
218          } else if (hasline[2]==1) {          } else { // hasline[2]==1
219              if (hasclass[7]==1)              if (hasclass[6])
220                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
221              else              else
222                  resultcode=FaceElements;                  resultcode=FaceElements;
223          } else          }
224              throw RipleyException("Internal Ripley Error!");      } else { // numLines==0
225      } else { // totlines==0          if (hasclass[1])
226          if (hasclass[2]==1)              // something from class 1
             // something from class 2  
227              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
228          else // something from class 1          else // something from class 0
229              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
230      }      }
231      return true;      return true;
232  }  }
233    
234    bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
235                                                  int fsType_target) const
236    {
237        if (!isValidFunctionSpaceType(fsType_target)) {
238            stringstream msg;
239            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
240                << fsType_target << " for " << getDescription();
241            throw RipleyException(msg.str());
242        }
243    
244        switch (fsType_source) {
245            case Nodes:
246            case DegreesOfFreedom:
247                return true;
248            case ReducedNodes:
249            case ReducedDegreesOfFreedom:
250                return (fsType_target != Nodes &&
251                        fsType_target != DegreesOfFreedom);
252            case Elements:
253                return (fsType_target==Elements ||
254                        fsType_target==ReducedElements);
255            case FaceElements:
256                return (fsType_target==FaceElements ||
257                        fsType_target==ReducedFaceElements);
258            case ReducedElements:
259            case ReducedFaceElements:
260            case Points:
261                return (fsType_target==fsType_source);
262    
263            default: {
264                stringstream msg;
265                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
266                    << fsType_source << " for " << getDescription();
267                throw RipleyException(msg.str());
268            }
269        }
270    }
271    
272  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
273                                         const escript::Data& in) const                                         const escript::Data& in) const
274  {  {
# Line 241  void RipleyDomain::interpolateOnDomain(e Line 285  void RipleyDomain::interpolateOnDomain(e
285          << " -> "          << " -> "
286          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
287    
288      switch (in.getFunctionSpace().getTypeCode()) {      const int inFS = in.getFunctionSpace().getTypeCode();
289          case Nodes:      const int outFS = target.getFunctionSpace().getTypeCode();
         case ReducedNodes:  
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
             switch (target.getFunctionSpace().getTypeCode()) {  
                 case Nodes:  
                 case ReducedNodes:  
                 case DegreesOfFreedom:  
                 case ReducedDegreesOfFreedom:  
                     // FIXME  
                     copyNodalData(target, *const_cast<escript::Data*>(&in));  
                     break;  
   
                 case Elements:  
                     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);  
                     break;  
   
                 case ReducedElements:  
                     interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);  
                     break;  
   
                 case FaceElements:  
                     interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);  
                     break;  
290    
291                  case ReducedFaceElements:      // simplest case: 1:1 copy
292                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);      if (inFS==outFS) {
293                      break;          copyData(target, *const_cast<escript::Data*>(&in));
294        // not allowed: reduced->non-reduced
295        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
296                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
297            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
298        } else if ((inFS==Elements && outFS==ReducedElements)
299                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
300            averageData(target, *const_cast<escript::Data*>(&in));
301        } else {
302            switch (inFS) {
303                case Nodes:
304                case ReducedNodes: //FIXME: reduced
305                    switch (outFS) {
306                        case DegreesOfFreedom:
307                        case ReducedDegreesOfFreedom: //FIXME: reduced
308                            if (getMPISize()==1)
309                                copyData(target, *const_cast<escript::Data*>(&in));
310                            else
311                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
312                            break;
313    
314                        case Nodes:
315                            copyData(target, *const_cast<escript::Data*>(&in));
316                        case ReducedNodes: //FIXME: reduced
317                            break;
318    
319                        case Elements:
320                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
321                            break;
322    
323                        case ReducedElements:
324                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
325                            break;
326    
327                        case FaceElements:
328                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
329                            break;
330    
331                        case ReducedFaceElements:
332                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
333                            break;
334    
335                        default:
336                            throw RipleyException(msg.str());
337                    }
338                    break;
339    
340                  default:              case DegreesOfFreedom:
341                      throw RipleyException(msg.str());              case ReducedDegreesOfFreedom: //FIXME: reduced
342              }                  switch (outFS) {
343              break;                      case Nodes:
344          default:                      case ReducedNodes: //FIXME: reduced
345              throw RipleyException(msg.str());                          if (getMPISize()==1)
346                                copyData(target, *const_cast<escript::Data*>(&in));
347                            else
348                                dofToNodes(target, *const_cast<escript::Data*>(&in));
349                            break;
350    
351                        case DegreesOfFreedom:
352                        case ReducedDegreesOfFreedom: //FIXME: reduced
353                            copyData(target, *const_cast<escript::Data*>(&in));
354                            break;
355    
356                        case Elements:
357                        case ReducedElements:
358                            if (getMPISize()==1) {
359                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
360                            } else {
361                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
362                                            escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
363                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
364                            }
365                            break;
366    
367                        case FaceElements:
368                        case ReducedFaceElements:
369                            if (getMPISize()==1) {
370                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
371                            } else {
372                                escript::Data contIn(in, (inFS==DegreesOfFreedom ?
373                                         escript::continuousFunction(*this) : escript::reducedContinuousFunction(*this)));
374                                interpolateNodesOnFaces(target, contIn, outFS==ReducedFaceElements);
375                            }
376                            break;
377    
378                        default:
379                            throw RipleyException(msg.str());
380                    }
381                    break;
382                default:
383                    throw RipleyException(msg.str());
384            }
385      }      }
386  }  }
387    
# Line 314  void RipleyDomain::setToX(escript::Data& Line 419  void RipleyDomain::setToX(escript::Data&
419      }      }
420  }  }
421    
422    void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const
423    {
424        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
425                *(arg.getFunctionSpace().getDomain()));
426        if (argDomain != *this)
427            throw RipleyException("setToGradient: Illegal domain of gradient argument");
428        const RipleyDomain& gradDomain=dynamic_cast<const RipleyDomain&>(
429                *(grad.getFunctionSpace().getDomain()));
430        if (gradDomain != *this)
431            throw RipleyException("setToGradient: Illegal domain of gradient");
432    
433        switch (grad.getFunctionSpace().getTypeCode()) {
434            case Elements:
435            case ReducedElements:
436            case FaceElements:
437            case ReducedFaceElements:
438                break;
439            default: {
440                stringstream msg;
441                msg << "setToGradient(): not supported for "
442                    << functionSpaceTypeAsString(grad.getFunctionSpace().getTypeCode());
443                throw RipleyException(msg.str());
444            }
445        }
446    
447        if (getMPISize()>1) {
448            if (arg.getFunctionSpace().getTypeCode()==DegreesOfFreedom) {
449                escript::Data contArg(arg, escript::continuousFunction(*this));
450                assembleGradient(grad, contArg);
451            } else if (arg.getFunctionSpace().getTypeCode()==ReducedDegreesOfFreedom) {
452                escript::Data contArg(arg, escript::reducedContinuousFunction(*this));
453                assembleGradient(grad, contArg);
454            } else {
455                assembleGradient(grad, *const_cast<escript::Data*>(&arg));
456            }
457        } else {
458            assembleGradient(grad, *const_cast<escript::Data*>(&arg));
459        }
460    }
461    
462    void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const
463    {
464        const RipleyDomain& argDomain=dynamic_cast<const RipleyDomain&>(
465                *(arg.getFunctionSpace().getDomain()));
466        if (argDomain != *this)
467            throw RipleyException("setToIntegrals: Illegal domain of integration kernel");
468    
469        switch (arg.getFunctionSpace().getTypeCode()) {
470            case Nodes:
471            case ReducedNodes:
472            case DegreesOfFreedom:
473            case ReducedDegreesOfFreedom:
474                {
475                    escript::Data funcArg(arg, escript::function(*this));
476                    assembleIntegrate(integrals, funcArg);
477                }
478                break;
479            case Elements:
480            case ReducedElements:
481            case FaceElements:
482            case ReducedFaceElements:
483                assembleIntegrate(integrals, *const_cast<escript::Data*>(&arg));
484                break;
485            default: {
486                stringstream msg;
487                msg << "setToIntegrals(): not supported for "
488                    << functionSpaceTypeAsString(arg.getFunctionSpace().getTypeCode());
489                throw RipleyException(msg.str());
490            }
491        }
492    
493    }
494    
495  bool RipleyDomain::isCellOriented(int fsType) const  bool RipleyDomain::isCellOriented(int fsType) const
496  {  {
497      switch(fsType) {      switch(fsType) {
498          case Nodes:          case Nodes:
499            case ReducedNodes:
500          case DegreesOfFreedom:          case DegreesOfFreedom:
501          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
502              return false;              return false;
# Line 331  bool RipleyDomain::isCellOriented(int fs Line 510  bool RipleyDomain::isCellOriented(int fs
510              break;              break;
511      }      }
512      stringstream msg;      stringstream msg;
513      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
514            << " on " << getDescription();
515      throw RipleyException(msg.str());      throw RipleyException(msg.str());
516  }  }
517    
# Line 347  bool RipleyDomain::canTag(int fsType) co Line 527  bool RipleyDomain::canTag(int fsType) co
527          case DegreesOfFreedom:          case DegreesOfFreedom:
528          case ReducedDegreesOfFreedom:          case ReducedDegreesOfFreedom:
529          case Points:          case Points:
530            case ReducedNodes:
531              return false;              return false;
532          default:          default:
533              break;              break;
534      }      }
535      stringstream msg;      stringstream msg;
536      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag(): Illegal function space type " << fsType << " on "
537            << getDescription();
538      throw RipleyException(msg.str());      throw RipleyException(msg.str());
539  }  }
540    
# Line 378  void RipleyDomain::setTags(const int fsT Line 560  void RipleyDomain::setTags(const int fsT
560              break;              break;
561          default: {          default: {
562              stringstream msg;              stringstream msg;
563              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags(): not implemented for "
564                    << functionSpaceTypeAsString(fsType);
565              throw RipleyException(msg.str());              throw RipleyException(msg.str());
566          }          }
567      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 598  int RipleyDomain::getTagFromSampleNo(int
598              break;              break;
599          default: {          default: {
600              stringstream msg;              stringstream msg;
601              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
602                    << functionSpaceTypeAsString(fsType);
603              throw RipleyException(msg.str());              throw RipleyException(msg.str());
604          }          }
605      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 619  int RipleyDomain::getNumberOfTagsInUse(i
619              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
620          default: {          default: {
621              stringstream msg;              stringstream msg;
622              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
623                    << functionSpaceTypeAsString(fsType);
624              throw RipleyException(msg.str());              throw RipleyException(msg.str());
625          }          }
626      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 639  const int* RipleyDomain::borrowListOfTag
639              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
640          default: {          default: {
641              stringstream msg;              stringstream msg;
642              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
643                    << functionSpaceTypeAsString(fsType);
644              throw RipleyException(msg.str());              throw RipleyException(msg.str());
645          }          }
646      }      }
# Line 501  escript::ASM_ptr RipleyDomain::newSystem Line 687  escript::ASM_ptr RipleyDomain::newSystem
687      bool reduceColOrder=false;      bool reduceColOrder=false;
688      // is the domain right?      // is the domain right?
689      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
690      if (row_domain!=*this)      if (row_domain != *this)
691          throw RipleyException("Domain of row function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of row function space does not match the domain of matrix generator");
692      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
693      if (col_domain!=*this)      if (col_domain != *this)
694          throw RipleyException("Domain of column function space does not match the domain of matrix generator");          throw RipleyException("newSystemMatrix(): Domain of column function space does not match the domain of matrix generator");
695      // is the function space type right?      // is the function space type right?
696      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
697          reduceRowOrder=true;          reduceRowOrder=true;
698      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
699          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
700      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
701          reduceColOrder=true;          reduceColOrder=true;
702      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
703          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
704    
705      // generate matrix      // generate matrix
706      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 527  escript::ASM_ptr RipleyDomain::newSystem Line 713  escript::ASM_ptr RipleyDomain::newSystem
713      return sma;      return sma;
714  }  }
715    
716    void RipleyDomain::addPDEToSystem(
717            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
718            const escript::Data& A, const escript::Data& B, const escript::Data& C,
719            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
720            const escript::Data& d, const escript::Data& y,
721            const escript::Data& d_contact, const escript::Data& y_contact,
722            const escript::Data& d_dirac,const escript::Data& y_dirac) const
723    {
724        if (!d_contact.isEmpty() || !y_contact.isEmpty())
725            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
726    
727        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
728        if (!sma)
729            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
730    
731        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
732            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
733    
734        int fsType=UNKNOWN;
735        fsType=(A.isEmpty() ? fsType : A.getFunctionSpace().getTypeCode());
736        fsType=(B.isEmpty() ? fsType : B.getFunctionSpace().getTypeCode());
737        fsType=(C.isEmpty() ? fsType : C.getFunctionSpace().getTypeCode());
738        fsType=(D.isEmpty() ? fsType : D.getFunctionSpace().getTypeCode());
739        fsType=(X.isEmpty() ? fsType : X.getFunctionSpace().getTypeCode());
740        fsType=(Y.isEmpty() ? fsType : Y.getFunctionSpace().getTypeCode());
741        fsType=(d.isEmpty() ? fsType : d.getFunctionSpace().getTypeCode());
742        fsType=(y.isEmpty() ? fsType : y.getFunctionSpace().getTypeCode());
743        fsType=(d_dirac.isEmpty() ? fsType : d_dirac.getFunctionSpace().getTypeCode());
744        fsType=(y_dirac.isEmpty() ? fsType : y_dirac.getFunctionSpace().getTypeCode());
745        if (fsType==UNKNOWN)
746            return;
747    
748        const bool reducedOrder=(fsType==ReducedElements || fsType==ReducedFaceElements);
749    
750        //TODO: more input param checks (shape, etc)
751    
752        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
753    
754        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
755            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
756    
757        const int numEq=S->logical_row_block_size;
758        const int numComp=S->logical_col_block_size;
759    
760        if (numEq != numComp)
761            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
762        //TODO: more system matrix checks
763    
764        if (numEq==1) {
765            if (reducedOrder) {
766                assemblePDESingleReduced(S, rhs, A, B, C, D, X, Y);
767                assemblePDEBoundarySingleReduced(S, rhs, escript::Data(),
768                        escript::Data(), escript::Data(), d, escript::Data(), y);
769            } else {
770                assemblePDESingle(S, rhs, A, B, C, D, X, Y);
771                assemblePDEBoundarySingle(S, rhs, escript::Data(),
772                        escript::Data(), escript::Data(), d, escript::Data(), y);
773            }
774        } else {
775            if (reducedOrder) {
776                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y);
777                assemblePDEBoundarySystemReduced(S, rhs, escript::Data(),
778                        escript::Data(), escript::Data(), d, escript::Data(), y);
779            } else {
780                assemblePDESystem(S, rhs, A, B, C, D, X, Y);
781                assemblePDEBoundarySystem(S, rhs, escript::Data(),
782                        escript::Data(), escript::Data(), d, escript::Data(), y);
783            }
784        }
785    }
786    
787    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
788            const escript::Data& Y, const escript::Data& y,
789            const escript::Data& y_contact, const escript::Data& y_dirac) const
790    {
791        if (!y_contact.isEmpty())
792            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
793    
794        if (rhs.isEmpty()) {
795            if (!X.isEmpty() || !Y.isEmpty())
796                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
797            else
798                return;
799        }
800    
801        if (rhs.getDataPointSize() == 1) {
802            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
803            assemblePDEBoundarySingle(NULL, rhs, escript::Data(), y);
804        } else {
805            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
806            assemblePDEBoundarySystem(NULL, rhs, escript::Data(), y);
807        }
808    }
809    
810  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
811  {  {
812      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
813  }  }
814    
815  //protected  //protected
816  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
817    {
818        const dim_t numComp = in.getDataPointSize();
819        const dim_t dpp = in.getNumDataPointsPerSample();
820        out.requireWrite();
821    #pragma omp parallel for
822        for (index_t i=0; i<in.getNumSamples(); i++) {
823            const double* src = in.getSampleDataRO(i);
824            double* dest = out.getSampleDataRW(i);
825            for (index_t c=0; c<numComp; c++) {
826                double res=0.;
827                for (index_t q=0; q<dpp; q++)
828                    res+=src[c+q*numComp];
829                *dest++ = res/dpp;
830            }
831        }
832    }
833    
834    //protected
835    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
836  {  {
837      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
838        out.requireWrite();
839  #pragma omp parallel for  #pragma omp parallel for
840      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
841          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 893  void RipleyDomain::updateTagsInUse(int f
893          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
894          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);
895  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
896    
897            // if we found a new value add it to the tagsInUse vector
898          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
899              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
900              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 903  void RipleyDomain::updateTagsInUse(int f
903      }      }
904  }  }
905    
906    //protected
907    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
908            const IndexVector& index, const dim_t M, const dim_t N) const
909    {
910        // paso will manage the memory
911        index_t* indexC = MEMALLOC(index.size(), index_t);
912        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
913        copy(index.begin(), index.end(), indexC);
914        copy(ptr.begin(), ptr.end(), ptrC);
915        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
916    }
917    
918    //protected
919    Paso_Pattern* RipleyDomain::createMainPattern() const
920    {
921        IndexVector ptr(1,0);
922        IndexVector index;
923    
924        for (index_t i=0; i<getNumDOF(); i++) {
925            // add the DOF itself
926            index.push_back(i);
927            const dim_t num=insertNeighbourNodes(index, i);
928            ptr.push_back(ptr.back()+num+1);
929        }
930    
931        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
932    }
933    
934    //protected
935    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
936            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
937    {
938        IndexVector ptr(1,0);
939        IndexVector index;
940        for (index_t i=0; i<getNumDOF(); i++) {
941            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
942            ptr.push_back(ptr.back()+colIndices[i].size());
943        }
944    
945        const dim_t M=ptr.size()-1;
946        *colPattern=createPasoPattern(ptr, index, M, N);
947    
948        IndexVector rowPtr(1,0);
949        IndexVector rowIndex;
950        for (dim_t id=0; id<N; id++) {
951            dim_t n=0;
952            for (dim_t i=0; i<M; i++) {
953                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
954                    if (index[j]==id) {
955                        rowIndex.push_back(i);
956                        n++;
957                        break;
958                    }
959                }
960            }
961            rowPtr.push_back(rowPtr.back()+n);
962        }
963    
964        // M and N are now reversed
965        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
966    }
967    
968    //protected
969    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
970           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
971           dim_t num_Sol, const vector<double>& array) const
972    {
973        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
974        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
975        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
976        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
977    
978        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
979        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
980        double* mainBlock_val = mat->mainBlock->val;
981        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
982        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
983        double* col_coupleBlock_val = mat->col_coupleBlock->val;
984        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
985        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
986        double* row_coupleBlock_val = mat->row_coupleBlock->val;
987    
988        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
989            // down columns of array
990            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
991                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
992                // only look at the matrix rows stored on this processor
993                if (i_row < numMyRows) {
994                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
995                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
996                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
997                            if (i_col < numMyCols) {
998                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
999                                    if (mainBlock_index[k] == i_col) {
1000                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1001                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1002                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1003                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1004                                                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())];
1005                                            }
1006                                        }
1007                                        break;
1008                                    }
1009                                }
1010                            } else {
1011                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1012                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1013                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1014                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1015                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1016                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1017                                                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())];
1018                                            }
1019                                        }
1020                                        break;
1021                                    }
1022                                }
1023                            }
1024                        }
1025                    }
1026                } else {
1027                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1028                        // across rows of array
1029                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1030                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1031                            if (i_col < numMyCols) {
1032                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1033                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1034                                {
1035                                    if (row_coupleBlock_index[k] == i_col) {
1036                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1037                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1038                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1039                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1040                                            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())];
1041                                            }
1042                                        }
1043                                        break;
1044                                    }
1045                                }
1046                            }
1047                        }
1048                    }
1049                }
1050            }
1051        }
1052    }
1053    
1054  //  //
1055  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
1056  //  //
# Line 612  string RipleyDomain::getDescription() co Line 1060  string RipleyDomain::getDescription() co
1060      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
1061  }  }
1062    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
1063  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
1064  {  {
1065      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 1075  const int* RipleyDomain::borrowSampleRef
1075      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1076  }  }
1077    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
1078  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1079  {  {
1080      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 659  void RipleyDomain::setToSize(escript::Da Line 1096  void RipleyDomain::setToSize(escript::Da
1096      throw RipleyException("setToSize() not implemented");      throw RipleyException("setToSize() not implemented");
1097  }  }
1098    
 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");  
 }  
   
1099  bool RipleyDomain::ownSample(int fsType, index_t id) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1100  {  {
1101      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1102  }  }
1103    
 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");  
 }  
   
1104  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1105          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
1106          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 692  void RipleyDomain::addPDEToLumpedSystem( Line 1108  void RipleyDomain::addPDEToLumpedSystem(
1108      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("addPDEToLumpedSystem() not implemented");
1109  }  }
1110    
 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");  
 }  
   
1111  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1112          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1113          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 749  IndexVector RipleyDomain::getNodeDistrib Line 1158  IndexVector RipleyDomain::getNodeDistrib
1158      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1159  }  }
1160    
1161    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1162    {
1163        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1164    }
1165    
1166  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1167  {  {
1168      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1183  dim_t RipleyDomain::getNumNodes() const
1183      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1184  }  }
1185    
1186    dim_t RipleyDomain::getNumDOF() const
1187    {
1188        throw RipleyException("getNumDOF() not implemented");
1189    }
1190    
1191    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1192    {
1193        throw RipleyException("insertNeighbourNodes() not implemented");
1194    }
1195    
1196  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1197  {  {
1198      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1199  }  }
1200    
1201    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1202    {
1203        throw RipleyException("assembleGradient() not implemented");
1204    }
1205    
1206    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1207    {
1208        throw RipleyException("assembleIntegrate() not implemented");
1209    }
1210    
1211    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1212            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1213            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1214    {
1215        throw RipleyException("assemblePDESingle() not implemented");
1216    }
1217    
1218    void RipleyDomain::assemblePDEBoundarySingle(Paso_SystemMatrix* mat,
1219            escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1220            const escript::Data& c, const escript::Data& d, const escript::Data& x,
1221            const escript::Data& y) const
1222    {
1223        throw RipleyException("assemblePDEBoundarySingle() not implemented");
1224    }
1225    
1226    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1227            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1228            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1229            const escript::Data& Y) const
1230    {
1231        throw RipleyException("assemblePDESingleReduced() not implemented");
1232    }
1233    
1234    void RipleyDomain::assemblePDEBoundarySingleReduced(Paso_SystemMatrix* mat,
1235            escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1236            const escript::Data& c, const escript::Data& d, const escript::Data& x,
1237            const escript::Data& y) const
1238    {
1239        throw RipleyException("assemblePDEBoundarySingleReduced() not implemented");
1240    }
1241    
1242    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1243            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1244            const escript::Data& D, const escript::Data& X, const escript::Data& Y) const
1245    {
1246        throw RipleyException("assemblePDESystem() not implemented");
1247    }
1248    
1249    void RipleyDomain::assemblePDEBoundarySystem(Paso_SystemMatrix* mat,
1250            escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1251            const escript::Data& c, const escript::Data& d, const escript::Data& x,
1252            const escript::Data& y) const
1253    {
1254        throw RipleyException("assemblePDEBoundarySystem() not implemented");
1255    }
1256    
1257    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1258            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1259            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1260            const escript::Data& Y) const
1261    {
1262        throw RipleyException("assemblePDESystemReduced() not implemented");
1263    }
1264    
1265    void RipleyDomain::assemblePDEBoundarySystemReduced(Paso_SystemMatrix* mat,
1266            escript::Data& rhs, const escript::Data& a, const escript::Data& b,
1267            const escript::Data& c, const escript::Data& d, const escript::Data& x,
1268            const escript::Data& y) const
1269    {
1270        throw RipleyException("assemblePDEBoundarySystemReduced() not implemented");
1271    }
1272    
1273  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1274  {  {
1275      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1280  void RipleyDomain::interpolateNodesOnFac
1280      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1281  }  }
1282    
1283    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1284    {
1285        throw RipleyException("nodesToDOF() not implemented");
1286    }
1287    
1288    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1289    {
1290        throw RipleyException("dofToNodes() not implemented");
1291    }
1292    
1293  } // end of namespace ripley  } // end of namespace ripley
1294    

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

  ViewVC Help
Powered by ViewVC 1.1.26