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

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

  ViewVC Help
Powered by ViewVC 1.1.26