/[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 3766 by caltinay, Thu Jan 12 08:17:49 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:
         case ReducedNodes:  
         case DegreesOfFreedom:  
         case ReducedDegreesOfFreedom:  
103              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
104            case ReducedNodes: //FIXME: reduced
105                if (getCoarseMesh())
106                    return getCoarseMesh()->getDataShape(Nodes);
107                break;
108            case DegreesOfFreedom:
109            case ReducedDegreesOfFreedom: //FIXME: reduced
110                return pair<int,int>(1, getNumDOF());
111          case Elements:          case Elements:
112              return pair<int,int>(ptsPerSample, getNumElements());              return pair<int,int>(ptsPerSample, getNumElements());
113          case FaceElements:          case FaceElements:
# Line 101  pair<int,int> RipleyDomain::getDataShape Line 116  pair<int,int> RipleyDomain::getDataShape
116              return pair<int,int>(1, getNumElements());              return pair<int,int>(1, getNumElements());
117          case ReducedFaceElements:          case ReducedFaceElements:
118              return pair<int,int>(1, getNumFaceElements());              return pair<int,int>(1, getNumFaceElements());
             /*  
119          case Points:          case Points:
120              return pair<int,int>(1, getNumPoints());              return pair<int,int>(1, 0); //FIXME: dirac
             */  
121          default:          default:
122              break;              break;
123      }      }
124    
125      stringstream msg;      stringstream msg;
126      msg << "getDataShape(): Unsupported function space type "      msg << "getDataShape(): Unsupported function space type " << fsType
127          << functionSpaceTypeAsString(fsType) << " for " << getDescription();          << " for " << getDescription();
128      throw RipleyException(msg.str());      throw RipleyException(msg.str());
129  }  }
130    
# Line 131  bool RipleyDomain::commonFunctionSpace(c Line 144  bool RipleyDomain::commonFunctionSpace(c
144     /*     /*
145      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
146      interpolated back and forth):      interpolated back and forth):
147      class 1: DOF <-> Nodes      class 0: DOF <-> Nodes
148      class 2: ReducedDOF <-> ReducedNodes      class 1: ReducedDOF <-> ReducedNodes
149      class 3: Points      class 2: Points
150      class 4: Elements      class 3: Elements
151      class 5: ReducedElements      class 4: ReducedElements
152      class 6: FaceElements      class 5: FaceElements
153      class 7: ReducedFaceElements      class 6: ReducedFaceElements
     class 8: ContactElementZero <-> ContactElementOne  
     class 9: ReducedContactElementZero <-> ReducedContactElementOne  
154    
155      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
156      between lines.      between lines.
157      class 1 and 2 belong to all lines so aren't considered.      class 0 and 1 belong to all lines so aren't considered.
158      line 0: class 3      line 0: class 2
159      line 1: class 4,5      line 1: class 3,4
160      line 2: class 6,7      line 2: class 5,6
     line 3: class 8,9  
161    
162      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
163      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
164      one instance of Nodes.      one instance of Nodes.
165      */      */
166      if (fs.empty())      if (fs.empty())
167          return false;          return false;
168      vector<int> hasclass(10);      vector<bool> hasclass(7, false);
169      vector<int> hasline(4);      vector<int> hasline(3, 0);
170      bool hasnodes=false;      bool hasnodes=false;
171      bool hasrednodes=false;      bool hasrednodes=false;
172      for (int i=0;i<fs.size();++i) {      for (size_t i=0; i<fs.size(); ++i) {
173          switch (fs[i]) {          switch (fs[i]) {
174              case Nodes: hasnodes=true; // no break is deliberate              case Nodes: hasnodes=true; // fall through
175              case DegreesOfFreedom:              case DegreesOfFreedom:
176                  hasclass[1]=1;                  hasclass[0]=true;
177                  break;                  break;
178              case ReducedNodes: hasrednodes=true; // no break is deliberate              case ReducedNodes: hasrednodes=true; // fall through
179              case ReducedDegreesOfFreedom:              case ReducedDegreesOfFreedom:
180                  hasclass[2]=1;                  hasclass[1]=true;
181                  break;                  break;
182              case Points:              case Points:
183                  hasline[0]=1;                  hasline[0]=1;
184                  hasclass[3]=1;                  hasclass[2]=true;
185                  break;                  break;
186              case Elements:              case Elements:
187                  hasclass[4]=1;                  hasclass[3]=true;
188                  hasline[1]=1;                  hasline[1]=1;
189                  break;                  break;
190              case ReducedElements:              case ReducedElements:
191                  hasclass[5]=1;                  hasclass[4]=true;
192                  hasline[1]=1;                  hasline[1]=1;
193                  break;                  break;
194              case FaceElements:              case FaceElements:
195                  hasclass[6]=1;                  hasclass[5]=true;
196                  hasline[2]=1;                  hasline[2]=1;
197                  break;                  break;
198              case ReducedFaceElements:              case ReducedFaceElements:
199                  hasclass[7]=1;                  hasclass[6]=true;
200                  hasline[2]=1;                  hasline[2]=1;
201                  break;                  break;
202              default:              default:
203                  return false;                  return false;
204          }          }
205      }      }
206      int totlines=hasline[0]+hasline[1]+hasline[2]+hasline[3];      int numLines=hasline[0]+hasline[1]+hasline[2];
207    
208      // fail if we have more than one leaf group      // fail if we have more than one leaf group
209      // = there are at least two branches we can't interpolate between      // = there are at least two branches we can't interpolate between
210      if (totlines>1)      if (numLines > 1)
211          return false;          return false;
212      else if (totlines==1) {      else if (numLines==1) {
213          // we have points          // we have points
214          if (hasline[0]==1)          if (hasline[0]==1)
215              resultcode=Points;              resultcode=Points;
216          else if (hasline[1]==1) {          else if (hasline[1]==1) {
217              if (hasclass[5]==1)              if (hasclass[4])
218                  resultcode=ReducedElements;                  resultcode=ReducedElements;
219              else              else
220                  resultcode=Elements;                  resultcode=Elements;
221          } else if (hasline[2]==1) {          } else { // hasline[2]==1
222              if (hasclass[7]==1)              if (hasclass[6])
223                  resultcode=ReducedFaceElements;                  resultcode=ReducedFaceElements;
224              else              else
225                  resultcode=FaceElements;                  resultcode=FaceElements;
226          } else          }
227              throw RipleyException("Internal Ripley Error!");      } else { // numLines==0
228      } else { // totlines==0          if (hasclass[1])
229          if (hasclass[2]==1)              // something from class 1
             // something from class 2  
230              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);              resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
231          else // something from class 1          else // something from class 0
232              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);              resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
233      }      }
234      return true;      return true;
235  }  }
236    
237    bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
238                                                  int fsType_target) const
239    {
240        if (!isValidFunctionSpaceType(fsType_target)) {
241            stringstream msg;
242            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
243                << fsType_target << " for " << getDescription();
244            throw RipleyException(msg.str());
245        }
246    
247        switch (fsType_source) {
248            case Nodes:
249            case DegreesOfFreedom:
250                return true;
251            case ReducedNodes:
252            case ReducedDegreesOfFreedom:
253                return (fsType_target != Nodes &&
254                        fsType_target != DegreesOfFreedom);
255            case Elements:
256                return (fsType_target==Elements ||
257                        fsType_target==ReducedElements);
258            case FaceElements:
259                return (fsType_target==FaceElements ||
260                        fsType_target==ReducedFaceElements);
261            case ReducedElements:
262            case ReducedFaceElements:
263            case Points:
264                return (fsType_target==fsType_source);
265    
266            default: {
267                stringstream msg;
268                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
269                    << fsType_source << " for " << getDescription();
270                throw RipleyException(msg.str());
271            }
272        }
273    }
274    
275  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
276                                         const escript::Data& in) const                                         const escript::Data& in) const
277  {  {
# Line 241  void RipleyDomain::interpolateOnDomain(e Line 288  void RipleyDomain::interpolateOnDomain(e
288          << " -> "          << " -> "
289          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
290    
291      switch (in.getFunctionSpace().getTypeCode()) {      const int inFS = in.getFunctionSpace().getTypeCode();
292          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;  
293    
294                  case ReducedFaceElements:      // simplest case: 1:1 copy
295                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);      if (inFS==outFS) {
296                      break;          copyData(target, *const_cast<escript::Data*>(&in));
297        // not allowed: reduced->non-reduced
298        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
299                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
300            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
301        } else if ((inFS==Elements && outFS==ReducedElements)
302                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
303            averageData(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                            copyData(target, *const_cast<escript::Data*>(&in));
319                        case ReducedNodes: //FIXME: reduced
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, d, y);
770            else
771                assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
772        else
773            if (reducedOrder)
774                assemblePDESystemReduced(S, rhs, A, B, C, D, X, Y, d, y);
775            else
776                assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
777    }
778    
779    void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,
780            const escript::Data& Y, const escript::Data& y,
781            const escript::Data& y_contact, const escript::Data& y_dirac) const
782    {
783        if (!y_contact.isEmpty())
784            throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
785    
786        if (rhs.isEmpty()) {
787            if (!X.isEmpty() || !Y.isEmpty())
788                throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
789            else
790                return;
791        }
792    
793        if (rhs.getDataPointSize() == 1)
794            assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
795        else
796            assemblePDESystem(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y, escript::Data(), y);
797    }
798    
799  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
800  {  {
801      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
802  }  }
803    
804  //protected  //protected
805  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
806    {
807        const dim_t numComp = in.getDataPointSize();
808        const dim_t dpp = in.getNumDataPointsPerSample();
809        out.requireWrite();
810    #pragma omp parallel for
811        for (index_t i=0; i<in.getNumSamples(); i++) {
812            const double* src = in.getSampleDataRO(i);
813            double* dest = out.getSampleDataRW(i);
814            for (index_t c=0; c<numComp; c++) {
815                double res=0.;
816                for (index_t q=0; q<dpp; q++)
817                    res+=src[c+q*numComp];
818                *dest++ = res/dpp;
819            }
820        }
821    }
822    
823    //protected
824    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
825  {  {
826      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
827        out.requireWrite();
828  #pragma omp parallel for  #pragma omp parallel for
829      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
830          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 882  void RipleyDomain::updateTagsInUse(int f
882          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
883          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);
884  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
885    
886            // if we found a new value add it to the tagsInUse vector
887          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
888              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
889              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 892  void RipleyDomain::updateTagsInUse(int f
892      }      }
893  }  }
894    
895    //protected
896    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
897            const IndexVector& index, const dim_t M, const dim_t N) const
898    {
899        // paso will manage the memory
900        index_t* indexC = MEMALLOC(index.size(), index_t);
901        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
902        copy(index.begin(), index.end(), indexC);
903        copy(ptr.begin(), ptr.end(), ptrC);
904        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
905    }
906    
907    //protected
908    Paso_Pattern* RipleyDomain::createMainPattern() const
909    {
910        IndexVector ptr(1,0);
911        IndexVector index;
912    
913        for (index_t i=0; i<getNumDOF(); i++) {
914            // add the DOF itself
915            index.push_back(i);
916            const dim_t num=insertNeighbourNodes(index, i);
917            ptr.push_back(ptr.back()+num+1);
918        }
919    
920        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
921    }
922    
923    //protected
924    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
925            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
926    {
927        IndexVector ptr(1,0);
928        IndexVector index;
929        for (index_t i=0; i<getNumDOF(); i++) {
930            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
931            ptr.push_back(ptr.back()+colIndices[i].size());
932        }
933    
934        const dim_t M=ptr.size()-1;
935        *colPattern=createPasoPattern(ptr, index, M, N);
936    
937        IndexVector rowPtr(1,0);
938        IndexVector rowIndex;
939        for (dim_t id=0; id<N; id++) {
940            dim_t n=0;
941            for (dim_t i=0; i<M; i++) {
942                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
943                    if (index[j]==id) {
944                        rowIndex.push_back(i);
945                        n++;
946                        break;
947                    }
948                }
949            }
950            rowPtr.push_back(rowPtr.back()+n);
951        }
952    
953        // M and N are now reversed
954        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
955    }
956    
957    //protected
958    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
959           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
960           dim_t num_Sol, const vector<double>& array) const
961    {
962        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
963        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
964        const dim_t numSubblocks_Eq = num_Eq / mat->row_block_size;
965        const dim_t numSubblocks_Sol = num_Sol / mat->col_block_size;
966    
967        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
968        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
969        double* mainBlock_val = mat->mainBlock->val;
970        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
971        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
972        double* col_coupleBlock_val = mat->col_coupleBlock->val;
973        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
974        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
975        double* row_coupleBlock_val = mat->row_coupleBlock->val;
976    
977        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
978            // down columns of array
979            for (dim_t l_row = 0; l_row < numSubblocks_Eq; ++l_row) {
980                const dim_t i_row = nodes_Eq[k_Eq]*numSubblocks_Eq+l_row;
981                // only look at the matrix rows stored on this processor
982                if (i_row < numMyRows) {
983                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
984                        for (dim_t l_col = 0; l_col < numSubblocks_Sol; ++l_col) {
985                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
986                            if (i_col < numMyCols) {
987                                for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
988                                    if (mainBlock_index[k] == i_col) {
989                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
990                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
991                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
992                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
993                                                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())];
994                                            }
995                                        }
996                                        break;
997                                    }
998                                }
999                            } else {
1000                                for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
1001                                    if (col_coupleBlock_index[k] == i_col - numMyCols) {
1002                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1003                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1004                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1005                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1006                                                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())];
1007                                            }
1008                                        }
1009                                        break;
1010                                    }
1011                                }
1012                            }
1013                        }
1014                    }
1015                } else {
1016                    for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
1017                        // across rows of array
1018                        for (dim_t l_col=0; l_col<numSubblocks_Sol; ++l_col) {
1019                            const dim_t i_col = nodes_Sol[k_Sol]*numSubblocks_Sol+l_col;
1020                            if (i_col < numMyCols) {
1021                                for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
1022                                     k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
1023                                {
1024                                    if (row_coupleBlock_index[k] == i_col) {
1025                                        for (dim_t ic=0; ic<mat->col_block_size; ++ic) {
1026                                            const dim_t i_Sol=ic+mat->col_block_size*l_col;
1027                                            for (dim_t ir=0; ir<mat->row_block_size; ++ir) {
1028                                                const dim_t i_Eq=ir+mat->row_block_size*l_row;
1029                                            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())];
1030                                            }
1031                                        }
1032                                        break;
1033                                    }
1034                                }
1035                            }
1036                        }
1037                    }
1038                }
1039            }
1040        }
1041    }
1042    
1043  //  //
1044  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
1045  //  //
# Line 612  string RipleyDomain::getDescription() co Line 1049  string RipleyDomain::getDescription() co
1049      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
1050  }  }
1051    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
1052  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
1053  {  {
1054      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 1064  const int* RipleyDomain::borrowSampleRef
1064      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1065  }  }
1066    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
1067  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1068  {  {
1069      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 659  void RipleyDomain::setToSize(escript::Da Line 1085  void RipleyDomain::setToSize(escript::Da
1085      throw RipleyException("setToSize() not implemented");      throw RipleyException("setToSize() not implemented");
1086  }  }
1087    
 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");  
 }  
   
1088  bool RipleyDomain::ownSample(int fsType, index_t id) const  bool RipleyDomain::ownSample(int fsType, index_t id) const
1089  {  {
1090      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
1091  }  }
1092    
 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");  
 }  
   
1093  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
1094          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
1095          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 692  void RipleyDomain::addPDEToLumpedSystem( Line 1097  void RipleyDomain::addPDEToLumpedSystem(
1097      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("addPDEToLumpedSystem() not implemented");
1098  }  }
1099    
 void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,  
         const escript::Data& Y, const escript::Data& y,  
         const escript::Data& y_contact, const escript::Data& y_dirac) const  
 {  
     throw RipleyException("addPDEToRHS() not implemented");  
 }  
   
1100  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
1101          escript::AbstractTransportProblem& tp,          escript::AbstractTransportProblem& tp,
1102          escript::Data& source, const escript::Data& M,          escript::Data& source, const escript::Data& M,
# Line 749  IndexVector RipleyDomain::getNodeDistrib Line 1147  IndexVector RipleyDomain::getNodeDistrib
1147      throw RipleyException("getNodeDistribution() not implemented");      throw RipleyException("getNodeDistribution() not implemented");
1148  }  }
1149    
1150    IndexVector RipleyDomain::getNumSubdivisionsPerDim() const
1151    {
1152        throw RipleyException("getNumSubdivisionsPerDim() not implemented");
1153    }
1154    
1155  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const  pair<double,double> RipleyDomain::getFirstCoordAndSpacing(dim_t dim) const
1156  {  {
1157      throw RipleyException("getFirstCoordAndSpacing() not implemented");      throw RipleyException("getFirstCoordAndSpacing() not implemented");
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1172  dim_t RipleyDomain::getNumNodes() const
1172      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1173  }  }
1174    
1175    dim_t RipleyDomain::getNumDOF() const
1176    {
1177        throw RipleyException("getNumDOF() not implemented");
1178    }
1179    
1180    escript::Domain_ptr RipleyDomain::getCoarseMesh() const
1181    {
1182        throw RipleyException("getCoarseMesh() not implemented");
1183    }
1184    
1185    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1186    {
1187        throw RipleyException("insertNeighbourNodes() not implemented");
1188    }
1189    
1190  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1191  {  {
1192      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1193  }  }
1194    
1195    void RipleyDomain::assembleGradient(escript::Data& out, escript::Data& in) const
1196    {
1197        throw RipleyException("assembleGradient() not implemented");
1198    }
1199    
1200    void RipleyDomain::assembleIntegrate(vector<double>& integrals, escript::Data& arg) const
1201    {
1202        throw RipleyException("assembleIntegrate() not implemented");
1203    }
1204    
1205    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1206            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1207            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1208            const escript::Data& d, const escript::Data& y) const
1209    {
1210        throw RipleyException("assemblePDESingle() not implemented");
1211    }
1212    
1213    void RipleyDomain::assemblePDESingleReduced(Paso_SystemMatrix* mat,
1214            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1215            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1216            const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1217    {
1218        throw RipleyException("assemblePDESingleReduced() not implemented");
1219    }
1220    
1221    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1222            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1223            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1224            const escript::Data& d, const escript::Data& y) const
1225    {
1226        throw RipleyException("assemblePDESystem() not implemented");
1227    }
1228    
1229    void RipleyDomain::assemblePDESystemReduced(Paso_SystemMatrix* mat,
1230            escript::Data& rhs, const escript::Data& A, const escript::Data& B,
1231            const escript::Data& C, const escript::Data& D, const escript::Data& X,
1232            const escript::Data& Y, const escript::Data& d, const escript::Data& y) const
1233    {
1234        throw RipleyException("assemblePDESystemReduced() not implemented");
1235    }
1236    
1237  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1238  {  {
1239      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1244  void RipleyDomain::interpolateNodesOnFac
1244      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1245  }  }
1246    
1247    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1248    {
1249        throw RipleyException("nodesToDOF() not implemented");
1250    }
1251    
1252    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1253    {
1254        throw RipleyException("dofToNodes() not implemented");
1255    }
1256    
1257  } // end of namespace ripley  } // end of namespace ripley
1258    

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

  ViewVC Help
Powered by ViewVC 1.1.26