/[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 3691 by caltinay, Wed Nov 23 23:07:37 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: //FIXME: reduced
104              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
105          case DegreesOfFreedom:          case DegreesOfFreedom:
106              return pair<int,int>(1, getNumNodes());          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:
111              return pair<int,int>(ptsPerSample/2, getNumFaceElements());              return pair<int,int>(ptsPerSample/2, getNumFaceElements());
             /*  
         case ReducedNodes:  
             return pair<int,int>(1, getNumReducedNodes());  
112          case ReducedElements:          case ReducedElements:
113              return pair<int,int>(1, getNumReducedElements());              return pair<int,int>(1, getNumElements());
114          case ReducedFaceElements:          case ReducedFaceElements:
115              return pair<int,int>(1, getNumReducedFaceElements());              return pair<int,int>(1, getNumFaceElements());
116          case Points:          case Points:
117              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, 0); //FIXME: dirac
         case ReducedDegreesOfFreedom:  
             return pair<int,int>(1, getNumReducedNodes());  
             */  
118          default:          default:
119              break;              break;
120      }      }
121    
122      stringstream msg;      stringstream msg;
123      msg << "Invalid function space type " << fsType << " for "      msg << "getDataShape(): Unsupported function space type " << fsType
124          << getDescription();          << " for " << getDescription();
125      throw RipleyException(msg.str());      throw RipleyException(msg.str());
126  }  }
127    
 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const  
 {  
     throw RipleyException("getTagFromSampleNo() not implemented");  
 }  
   
128  string RipleyDomain::showTagNames() const  string RipleyDomain::showTagNames() const
129  {  {
130      stringstream ret;      stringstream ret;
# Line 139  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,
273                                           const escript::Data& in) const
274    {
275        const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
276        const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
277        if (inDomain != *this)
278            throw RipleyException("Illegal domain of interpolant");
279        if (targetDomain != *this)
280            throw RipleyException("Illegal domain of interpolation target");
281    
282        stringstream msg;
283        msg << "interpolateOnDomain() not implemented for function space "
284            << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
285            << " -> "
286            << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
287    
288        const int inFS = in.getFunctionSpace().getTypeCode();
289        const int outFS = target.getFunctionSpace().getTypeCode();
290    
291        // simplest case: 1:1 copy
292        if (inFS==outFS) {
293            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                case DegreesOfFreedom:
344                case ReducedDegreesOfFreedom: //FIXME: reduced
345                    switch (outFS) {
346                        case Nodes:
347                        case ReducedNodes: //FIXME: reduced
348                            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    
391  escript::Data RipleyDomain::getX() const  escript::Data RipleyDomain::getX() const
392  {  {
393      return escript::continuousFunction(*this).getX();      return escript::continuousFunction(*this).getX();
# Line 267  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 284  bool RipleyDomain::isCellOriented(int fs Line 513  bool RipleyDomain::isCellOriented(int fs
513              break;              break;
514      }      }
515      stringstream msg;      stringstream msg;
516      msg << "Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
517            << " on " << getDescription();
518        throw RipleyException(msg.str());
519    }
520    
521    bool RipleyDomain::canTag(int fsType) const
522    {
523        switch(fsType) {
524            case Nodes:
525            case Elements:
526            case ReducedElements:
527            case FaceElements:
528            case ReducedFaceElements:
529                return true;
530            case DegreesOfFreedom:
531            case ReducedDegreesOfFreedom:
532            case Points:
533            case ReducedNodes:
534                return false;
535            default:
536                break;
537        }
538        stringstream msg;
539        msg << "canTag(): Illegal function space type " << fsType << " on "
540            << getDescription();
541      throw RipleyException(msg.str());      throw RipleyException(msg.str());
542  }  }
543    
544    void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
545    {
546        IndexVector* target=NULL;
547        dim_t num=0;
548    
549        switch(fsType) {
550            case Nodes:
551                num=getNumNodes();
552                target=&m_nodeTags;
553                break;
554            case Elements:
555            case ReducedElements:
556                num=getNumElements();
557                target=&m_elementTags;
558                break;
559            case FaceElements:
560            case ReducedFaceElements:
561                num=getNumFaceElements();
562                target=&m_faceTags;
563                break;
564            default: {
565                stringstream msg;
566                msg << "setTags(): not implemented for "
567                    << functionSpaceTypeAsString(fsType);
568                throw RipleyException(msg.str());
569            }
570        }
571        if (target->size() != num) {
572            target->assign(num, -1);
573        }
574    
575        escript::Data& mask=*const_cast<escript::Data*>(&cMask);
576    #pragma omp parallel for
577        for (index_t i=0; i<num; i++) {
578            if (mask.getSampleDataRO(i)[0] > 0) {
579                (*target)[i]=newTag;
580            }
581        }
582        updateTagsInUse(fsType);
583    }
584    
585    int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
586    {
587        switch(fsType) {
588            case Nodes:
589                if (m_nodeTags.size() > sampleNo)
590                    return m_nodeTags[sampleNo];
591                break;
592            case Elements:
593            case ReducedElements:
594                if (m_elementTags.size() > sampleNo)
595                    return m_elementTags[sampleNo];
596                break;
597            case FaceElements:
598            case ReducedFaceElements:
599                if (m_faceTags.size() > sampleNo)
600                    return m_faceTags[sampleNo];
601                break;
602            default: {
603                stringstream msg;
604                msg << "getTagFromSampleNo(): not implemented for "
605                    << functionSpaceTypeAsString(fsType);
606                throw RipleyException(msg.str());
607            }
608        }
609        return -1;
610    }
611    
612    int RipleyDomain::getNumberOfTagsInUse(int fsType) const
613    {
614        switch(fsType) {
615            case Nodes:
616                return m_nodeTagsInUse.size();
617            case Elements:
618            case ReducedElements:
619                return m_elementTagsInUse.size();
620            case FaceElements:
621            case ReducedFaceElements:
622                return m_faceTagsInUse.size();
623            default: {
624                stringstream msg;
625                msg << "getNumberOfTagsInUse(): not implemented for "
626                    << functionSpaceTypeAsString(fsType);
627                throw RipleyException(msg.str());
628            }
629        }
630    }
631    
632    const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
633    {
634        switch(fsType) {
635            case Nodes:
636                return &m_nodeTagsInUse[0];
637            case Elements:
638            case ReducedElements:
639                return &m_elementTagsInUse[0];
640            case FaceElements:
641            case ReducedFaceElements:
642                return &m_faceTagsInUse[0];
643            default: {
644                stringstream msg;
645                msg << "borrowListOfTagsInUse(): not implemented for "
646                    << functionSpaceTypeAsString(fsType);
647                throw RipleyException(msg.str());
648            }
649        }
650    }
651    
652  void RipleyDomain::Print_Mesh_Info(const bool full) const  void RipleyDomain::Print_Mesh_Info(const bool full) const
653  {  {
654      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
# Line 329  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  // the methods that follow have to be implemented by the subclasses          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  string RipleyDomain::getDescription() const          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      throw RipleyException("getDescription() not implemented");      if (!d_contact.isEmpty() || !y_contact.isEmpty())
727  }          throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
728    
729  bool RipleyDomain::operator==(const AbstractDomain& other) const      paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
730  {      if (!sma)
731      throw RipleyException("operator==() not implemented");          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::write(const string& filename) const  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      throw RipleyException("write() not implemented");      if (!y_contact.isEmpty())
790  }          throw RipleyException("addPDEToRHS(): Ripley does not support contact elements");
791    
792  void RipleyDomain::dump(const string& filename) const      if (rhs.isEmpty()) {
793  {          if (!X.isEmpty() || !Y.isEmpty())
794      throw RipleyException("dump() not implemented");              throw RipleyException("addPDEToRHS(): Right hand side coefficients are provided but no right hand side vector given");
795  }          else
796                return;
797        }
798    
799  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const      if (rhs.getDataPointSize() == 1) {
800  {          assemblePDESingle(NULL, rhs, escript::Data(), escript::Data(), escript::Data(), escript::Data(), X, Y);
801      throw RipleyException("borrowSampleReferenceIDs() not implemented");          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)
# Line 389  void RipleyDomain::setNewX(const escript Line 810  void RipleyDomain::setNewX(const escript
810      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
811  }  }
812    
813  void RipleyDomain::interpolateOnDomain(escript::Data& target,  //protected
814                                         const escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
815  {  {
816      throw RipleyException("interpolateOnDomain() not implemented");      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  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  //protected
833                                               int fsType_target) const  void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
834  {  {
835      throw RipleyException("probeInterpolationOnDomain() not implemented");      const dim_t numComp = in.getDataPointSize();
836        out.requireWrite();
837    #pragma omp parallel for
838        for (index_t i=0; i<in.getNumSamples(); i++) {
839            const double* src = in.getSampleDataRO(i);
840            copy(src, src+numComp, out.getSampleDataRW(i));
841        }
842  }  }
843    
844  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  //protected
845    void RipleyDomain::updateTagsInUse(int fsType) const
846  {  {
847      throw RipleyException("interpolateACross() not implemented");      IndexVector* tagsInUse=NULL;
848  }      const IndexVector* tags=NULL;
849        switch(fsType) {
850            case Nodes:
851                tags=&m_nodeTags;
852                tagsInUse=&m_nodeTagsInUse;
853                break;
854            case Elements:
855            case ReducedElements:
856                tags=&m_elementTags;
857                tagsInUse=&m_elementTagsInUse;
858                break;
859            case FaceElements:
860            case ReducedFaceElements:
861                tags=&m_faceTags;
862                tagsInUse=&m_faceTagsInUse;
863                break;
864            default:
865                return;
866        }
867    
868  bool RipleyDomain::probeInterpolationACross(int fsType_source,      // gather global unique tag values from tags into tagsInUse
869          const escript::AbstractDomain&, int fsType_target) const      tagsInUse->clear();
870  {      index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
871      throw RipleyException("probeInterpolationACross() not implemented");  
872        while (true) {
873            // find smallest value bigger than lastFoundValue
874            minFoundValue = INDEX_T_MAX;
875    #pragma omp parallel private(local_minFoundValue)
876            {
877                local_minFoundValue = minFoundValue;
878    #pragma omp for schedule(static) nowait
879                for (size_t i = 0; i < tags->size(); i++) {
880                    const index_t v = (*tags)[i];
881                    if ((v > lastFoundValue) && (v < local_minFoundValue))
882                        local_minFoundValue = v;
883                }
884    #pragma omp critical
885                {
886                    if (local_minFoundValue < minFoundValue)
887                        minFoundValue = local_minFoundValue;
888                }
889            }
890    #ifdef ESYS_MPI
891            local_minFoundValue = minFoundValue;
892            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
893    #endif
894    
895            // if we found a new value add it to the tagsInUse vector
896            if (minFoundValue < INDEX_T_MAX) {
897                tagsInUse->push_back(minFoundValue);
898                lastFoundValue = minFoundValue;
899            } else
900                break;
901        }
902  }  }
903    
904  void RipleyDomain::setToNormal(escript::Data& normal) const  //protected
905  {  Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
906      throw RipleyException("setToNormal() not implemented");          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  void RipleyDomain::setToSize(escript::Data& size) const  //protected
933  {  void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
934      throw RipleyException("setToSize() not implemented");          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  void RipleyDomain::setToGradient(escript::Data& grad, const escript::Data& arg) const  //protected
967  {  void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
968      throw RipleyException("setToGradient() not implemented");         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  void RipleyDomain::setToIntegrals(vector<double>& integrals, const escript::Data& arg) const  //
1053    // the methods that follow have to be implemented by the subclasses
1054    //
1055    
1056    string RipleyDomain::getDescription() const
1057  {  {
1058      throw RipleyException("setToIntegrals() not implemented");      throw RipleyException("getDescription() not implemented");
1059  }  }
1060    
1061  bool RipleyDomain::ownSample(int fsType, index_t id) const  void RipleyDomain::write(const string& filename) const
1062  {  {
1063      throw RipleyException("ownSample() not implemented");      throw RipleyException("write() not implemented");
1064  }  }
1065    
1066  void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const  void RipleyDomain::dump(const string& filename) const
1067  {  {
1068      throw RipleyException("setTags() not implemented");      throw RipleyException("dump() not implemented");
1069  }  }
1070    
1071  int RipleyDomain::getNumberOfTagsInUse(int fsType) const  const int* RipleyDomain::borrowSampleReferenceIDs(int fsType) const
1072  {  {
1073      throw RipleyException("getNumberOfTagsInUse() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
1074  }  }
1075    
1076  const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
1077  {  {
1078      throw RipleyException("borrowListOfTagsInUse() not implemented");      throw RipleyException("interpolateACross() not implemented");
1079  }  }
1080    
1081  bool RipleyDomain::canTag(int fsType) const  bool RipleyDomain::probeInterpolationACross(int fsType_source,
1082            const escript::AbstractDomain&, int fsType_target) const
1083  {  {
1084      throw RipleyException("canTag() not implemented");      throw RipleyException("probeInterpolationACross() not implemented");
1085  }  }
1086    
1087    void RipleyDomain::setToNormal(escript::Data& normal) const
 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  
1088  {  {
1089      throw RipleyException("addPDEToSystem() not implemented");      throw RipleyException("setToNormal() not implemented");
1090  }  }
1091    
1092  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::setToSize(escript::Data& size) const
         const escript::Data& D, const escript::Data& d,  
         const escript::Data& d_dirac, const bool useHRZ) const  
1093  {  {
1094      throw RipleyException("addPDEToLumpedSystem() not implemented");      throw RipleyException("setToSize() not implemented");
1095  }  }
1096    
1097  void RipleyDomain::addPDEToRHS(escript::Data& rhs, const escript::Data& X,  bool RipleyDomain::ownSample(int fsType, index_t id) const
         const escript::Data& Y, const escript::Data& y,  
         const escript::Data& y_contact, const escript::Data& y_dirac) const  
1098  {  {
1099      throw RipleyException("addPDEToRHS() not implemented");      throw RipleyException("ownSample() not implemented");
1100  }  }
1101    
1102  void RipleyDomain::addPDEToTransportProblem(  void RipleyDomain::addPDEToTransportProblem(
# Line 513  dim_t RipleyDomain::getNumDataPointsGlob Line 1129  dim_t RipleyDomain::getNumDataPointsGlob
1129      throw RipleyException("getNumDataPointsGlobal() not implemented");      throw RipleyException("getNumDataPointsGlobal() not implemented");
1130  }  }
1131    
1132    IndexVector RipleyDomain::getNumNodesPerDim() const
1133    {
1134        throw RipleyException("getNumNodesPerDim() not implemented");
1135    }
1136    
1137    IndexVector RipleyDomain::getNumElementsPerDim() const
1138    {
1139        throw RipleyException("getNumElementsPerDim() not implemented");
1140    }
1141    
1142    IndexVector RipleyDomain::getNumFacesPerBoundary() const
1143    {
1144        throw RipleyException("getNumFacesPerBoundary() not implemented");
1145    }
1146    
1147    IndexVector RipleyDomain::getNodeDistribution() const
1148    {
1149        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
1158    {
1159        throw RipleyException("getFirstCoordAndSpacing() not implemented");
1160    }
1161    
1162  dim_t RipleyDomain::getNumFaceElements() const  dim_t RipleyDomain::getNumFaceElements() const
1163  {  {
1164      throw RipleyException("getNumFaceElements() not implemented");      throw RipleyException("getNumFaceElements() not implemented");
# Line 528  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
1257    {
1258        throw RipleyException("interpolateNodesOnElements() not implemented");
1259    }
1260    
1261    void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
1262    {
1263        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.3691  
changed lines
  Added in v.3785

  ViewVC Help
Powered by ViewVC 1.1.26