/[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 3756 by caltinay, Fri Jan 6 02:35:19 2012 UTC
# 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 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 (fsType_target != Nodes &&
238                fsType_target != ReducedNodes &&
239                fsType_target != ReducedDegreesOfFreedom &&
240                fsType_target != DegreesOfFreedom &&
241                fsType_target != Elements &&
242                fsType_target != ReducedElements &&
243                fsType_target != FaceElements &&
244                fsType_target != ReducedFaceElements &&
245                fsType_target != Points) {
246            stringstream msg;
247            msg << "probeInterpolationOnDomain(): Invalid functionspace type "
248                << fsType_target << " for " << getDescription();
249            throw RipleyException(msg.str());
250        }
251    
252        switch (fsType_source) {
253            case Nodes:
254            case DegreesOfFreedom:
255                return true;
256            case ReducedNodes:
257            case ReducedDegreesOfFreedom:
258                return (fsType_target != Nodes &&
259                        fsType_target != DegreesOfFreedom);
260            case Elements:
261                return (fsType_target==Elements ||
262                        fsType_target==ReducedElements);
263            case ReducedElements:
264                return (fsType_target==ReducedElements);
265            case FaceElements:
266                return (fsType_target==FaceElements ||
267                        fsType_target==ReducedFaceElements);
268            case ReducedFaceElements:
269                return (fsType_target==ReducedFaceElements);
270            case Points:
271                return (fsType_target==Points);
272    
273            default: {
274                stringstream msg;
275                msg << "probeInterpolationOnDomain(): Invalid functionspace type "
276                    << fsType_source << " for " << getDescription();
277                throw RipleyException(msg.str());
278            }
279        }
280    }
281    
282  void RipleyDomain::interpolateOnDomain(escript::Data& target,  void RipleyDomain::interpolateOnDomain(escript::Data& target,
283                                         const escript::Data& in) const                                         const escript::Data& in) const
284  {  {
# Line 241  void RipleyDomain::interpolateOnDomain(e Line 295  void RipleyDomain::interpolateOnDomain(e
295          << " -> "          << " -> "
296          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());          << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
297    
298      switch (in.getFunctionSpace().getTypeCode()) {      const int inFS = in.getFunctionSpace().getTypeCode();
299          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;  
300    
301                  case ReducedFaceElements:      // simplest case: 1:1 copy
302                      interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);      if (inFS==outFS) {
303                      break;          copyData(target, *const_cast<escript::Data*>(&in));
304        // not allowed: reduced->non-reduced
305        } else if ((inFS==ReducedNodes || inFS==ReducedDegreesOfFreedom)
306                && (outFS==Nodes || outFS==DegreesOfFreedom)) {
307            throw RipleyException("interpolateOnDomain(): Cannot interpolate reduced data to non-reduced data.");
308        } else if ((inFS==Elements && outFS==ReducedElements)
309                || (inFS==FaceElements && outFS==ReducedFaceElements)) {
310            averageData(target, *const_cast<escript::Data*>(&in));
311        } else {
312            switch (inFS) {
313                case Nodes:
314                case ReducedNodes: //FIXME: reduced
315                    switch (outFS) {
316                        case DegreesOfFreedom:
317                        case ReducedDegreesOfFreedom: //FIXME: reduced
318                            if (getMPISize()==1)
319                                copyData(target, *const_cast<escript::Data*>(&in));
320                            else
321                                nodesToDOF(target,*const_cast<escript::Data*>(&in));
322                            break;
323    
324                        case Nodes:
325                        case ReducedNodes: //FIXME: reduced
326                            copyData(target, *const_cast<escript::Data*>(&in));
327                            break;
328    
329                        case Elements:
330                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
331                            break;
332    
333                        case ReducedElements:
334                            interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
335                            break;
336    
337                        case FaceElements:
338                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
339                            break;
340    
341                        case ReducedFaceElements:
342                            interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
343                            break;
344    
345                        default:
346                            throw RipleyException(msg.str());
347                    }
348                    break;
349    
350                  default:              case DegreesOfFreedom:
351                      throw RipleyException(msg.str());              case ReducedDegreesOfFreedom: //FIXME: reduced
352              }                  switch (outFS) {
353              break;                      case Nodes:
354          default:                      case ReducedNodes: //FIXME: reduced
355              throw RipleyException(msg.str());                          if (getMPISize()==1)
356                                copyData(target, *const_cast<escript::Data*>(&in));
357                            else
358                                dofToNodes(target, *const_cast<escript::Data*>(&in));
359                            break;
360    
361                        case DegreesOfFreedom:
362                        case ReducedDegreesOfFreedom: //FIXME: reduced
363                            copyData(target, *const_cast<escript::Data*>(&in));
364                            break;
365    
366                        case Elements:
367                        case ReducedElements:
368                            if (getMPISize()==1) {
369                                interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), outFS==ReducedElements);
370                            } else {
371                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
372                                interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
373                            }
374                            break;
375    
376                        case FaceElements:
377                        case ReducedFaceElements:
378                            if (getMPISize()==1) {
379                                interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), outFS==ReducedFaceElements);
380                            } else {
381                                escript::Data contIn=escript::Data(in, continuousFunction(*this));
382                                interpolateNodesOnElements(target, contIn, outFS==ReducedFaceElements);
383                            }
384                            break;
385    
386                        default:
387                            throw RipleyException(msg.str());
388                    }
389                    break;
390                default:
391                    throw RipleyException(msg.str());
392            }
393      }      }
394  }  }
395    
# Line 331  bool RipleyDomain::isCellOriented(int fs Line 444  bool RipleyDomain::isCellOriented(int fs
444              break;              break;
445      }      }
446      stringstream msg;      stringstream msg;
447      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType
448            << " on " << getDescription();
449      throw RipleyException(msg.str());      throw RipleyException(msg.str());
450  }  }
451    
# Line 352  bool RipleyDomain::canTag(int fsType) co Line 466  bool RipleyDomain::canTag(int fsType) co
466              break;              break;
467      }      }
468      stringstream msg;      stringstream msg;
469      msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();      msg << "canTag(): Illegal function space type " << fsType << " on "
470            << getDescription();
471      throw RipleyException(msg.str());      throw RipleyException(msg.str());
472  }  }
473    
# Line 378  void RipleyDomain::setTags(const int fsT Line 493  void RipleyDomain::setTags(const int fsT
493              break;              break;
494          default: {          default: {
495              stringstream msg;              stringstream msg;
496              msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "setTags(): not implemented for "
497                    << functionSpaceTypeAsString(fsType);
498              throw RipleyException(msg.str());              throw RipleyException(msg.str());
499          }          }
500      }      }
# Line 415  int RipleyDomain::getTagFromSampleNo(int Line 531  int RipleyDomain::getTagFromSampleNo(int
531              break;              break;
532          default: {          default: {
533              stringstream msg;              stringstream msg;
534              msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getTagFromSampleNo(): not implemented for "
535                    << functionSpaceTypeAsString(fsType);
536              throw RipleyException(msg.str());              throw RipleyException(msg.str());
537          }          }
538      }      }
# Line 435  int RipleyDomain::getNumberOfTagsInUse(i Line 552  int RipleyDomain::getNumberOfTagsInUse(i
552              return m_faceTagsInUse.size();              return m_faceTagsInUse.size();
553          default: {          default: {
554              stringstream msg;              stringstream msg;
555              msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "getNumberOfTagsInUse(): not implemented for "
556                    << functionSpaceTypeAsString(fsType);
557              throw RipleyException(msg.str());              throw RipleyException(msg.str());
558          }          }
559      }      }
# Line 454  const int* RipleyDomain::borrowListOfTag Line 572  const int* RipleyDomain::borrowListOfTag
572              return &m_faceTagsInUse[0];              return &m_faceTagsInUse[0];
573          default: {          default: {
574              stringstream msg;              stringstream msg;
575              msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);              msg << "borrowListOfTagsInUse(): not implemented for "
576                    << functionSpaceTypeAsString(fsType);
577              throw RipleyException(msg.str());              throw RipleyException(msg.str());
578          }          }
579      }      }
# Line 502  escript::ASM_ptr RipleyDomain::newSystem Line 621  escript::ASM_ptr RipleyDomain::newSystem
621      // is the domain right?      // is the domain right?
622      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));      const RipleyDomain& row_domain=dynamic_cast<const RipleyDomain&>(*(row_functionspace.getDomain()));
623      if (row_domain!=*this)      if (row_domain!=*this)
624          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");
625      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));      const RipleyDomain& col_domain=dynamic_cast<const RipleyDomain&>(*(column_functionspace.getDomain()));
626      if (col_domain!=*this)      if (col_domain!=*this)
627          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");
628      // is the function space type right?      // is the function space type right?
629      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (row_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
630          reduceRowOrder=true;          reduceRowOrder=true;
631      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (row_functionspace.getTypeCode()!=DegreesOfFreedom)
632          throw RipleyException("Illegal function space type for system matrix rows");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix rows");
633      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)      if (column_functionspace.getTypeCode()==ReducedDegreesOfFreedom)
634          reduceColOrder=true;          reduceColOrder=true;
635      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)      else if (column_functionspace.getTypeCode()!=DegreesOfFreedom)
636          throw RipleyException("Illegal function space type for system matrix columns");          throw RipleyException("newSystemMatrix(): Illegal function space type for system matrix columns");
637    
638      // generate matrix      // generate matrix
639      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);      Paso_SystemMatrixPattern* pattern=getPattern(reduceRowOrder, reduceColOrder);
# Line 527  escript::ASM_ptr RipleyDomain::newSystem Line 646  escript::ASM_ptr RipleyDomain::newSystem
646      return sma;      return sma;
647  }  }
648    
649    void RipleyDomain::addPDEToSystem(
650            escript::AbstractSystemMatrix& mat, escript::Data& rhs,
651            const escript::Data& A, const escript::Data& B, const escript::Data& C,
652            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
653            const escript::Data& d, const escript::Data& y,
654            const escript::Data& d_contact, const escript::Data& y_contact,
655            const escript::Data& d_dirac,const escript::Data& y_dirac) const
656    {
657        if (!d_contact.isEmpty() || !y_contact.isEmpty())
658            throw RipleyException("addPDEToSystem(): Ripley does not support contact elements");
659    
660        paso::SystemMatrixAdapter* sma=dynamic_cast<paso::SystemMatrixAdapter*>(&mat);
661        if (!sma)
662            throw RipleyException("addPDEToSystem(): Ripley only accepts Paso system matrices");
663    
664        if (rhs.isEmpty() && (!X.isEmpty() || !Y.isEmpty()))
665            throw RipleyException("addPDEToSystem(): Right hand side coefficients are provided but no right hand side vector given");
666    
667        //TODO: more input param checks (shape, function space etc)
668    
669        Paso_SystemMatrix* S = sma->getPaso_SystemMatrix();
670    
671        if (!rhs.isEmpty() && rhs.getDataPointSize() != S->logical_row_block_size)
672            throw RipleyException("addPDEToSystem(): Matrix row block size and number of components of right hand side don't match");
673    
674        const int numEq=S->logical_row_block_size;
675        const int numComp=S->logical_col_block_size;
676    
677        if (numEq != numComp)
678            throw RipleyException("addPDEToSystem(): Number of equations and number of solutions don't match");
679        //TODO: more system matrix checks
680    
681        if (numEq==1)
682            assemblePDESingle(S, rhs, A, B, C, D, X, Y, d, y);
683        else
684            assemblePDESystem(S, rhs, A, B, C, D, X, Y, d, y);
685    }
686    
687  void RipleyDomain::setNewX(const escript::Data& arg)  void RipleyDomain::setNewX(const escript::Data& arg)
688  {  {
689      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
690  }  }
691    
692  //protected  //protected
693  void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const  void RipleyDomain::averageData(escript::Data& out, escript::Data& in) const
694  {  {
695      const dim_t numComp = in.getDataPointSize();      const dim_t numComp = in.getDataPointSize();
696        const dim_t dpp = in.getNumDataPointsPerSample();
697        out.requireWrite();
698    #pragma omp parallel for
699        for (index_t i=0; i<in.getNumSamples(); i++) {
700            const double* src = in.getSampleDataRO(i);
701            double* dest = out.getSampleDataRW(i);
702            for (index_t c=0; c<numComp; c++) {
703                double res=0.;
704                for (index_t q=0; q<dpp; q++)
705                    res+=src[c+q*numComp];
706                *dest++ = res/dpp;
707            }
708        }
709    }
710    
711    //protected
712    void RipleyDomain::copyData(escript::Data& out, escript::Data& in) const
713    {
714        const dim_t numComp = in.getDataPointSize();
715        out.requireWrite();
716  #pragma omp parallel for  #pragma omp parallel for
717      for (index_t i=0; i<in.getNumSamples(); i++) {      for (index_t i=0; i<in.getNumSamples(); i++) {
718          const double* src = in.getSampleDataRO(i);          const double* src = in.getSampleDataRO(i);
# Line 593  void RipleyDomain::updateTagsInUse(int f Line 770  void RipleyDomain::updateTagsInUse(int f
770          local_minFoundValue = minFoundValue;          local_minFoundValue = minFoundValue;
771          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);
772  #endif  #endif
         // if we found a new value add it to the tagsInUse vector  
773    
774            // if we found a new value add it to the tagsInUse vector
775          if (minFoundValue < INDEX_T_MAX) {          if (minFoundValue < INDEX_T_MAX) {
776              tagsInUse->push_back(minFoundValue);              tagsInUse->push_back(minFoundValue);
777              lastFoundValue = minFoundValue;              lastFoundValue = minFoundValue;
# Line 603  void RipleyDomain::updateTagsInUse(int f Line 780  void RipleyDomain::updateTagsInUse(int f
780      }      }
781  }  }
782    
783    //protected
784    Paso_Pattern* RipleyDomain::createPasoPattern(const IndexVector& ptr,
785            const IndexVector& index, const dim_t M, const dim_t N) const
786    {
787        // paso will manage the memory
788        index_t* indexC = MEMALLOC(index.size(), index_t);
789        index_t* ptrC = MEMALLOC(ptr.size(), index_t);
790        copy(index.begin(), index.end(), indexC);
791        copy(ptr.begin(), ptr.end(), ptrC);
792        return Paso_Pattern_alloc(MATRIX_FORMAT_DEFAULT, M, N, ptrC, indexC);
793    }
794    
795    //protected
796    Paso_Pattern* RipleyDomain::createMainPattern() const
797    {
798        IndexVector ptr(1,0);
799        IndexVector index;
800    
801        for (index_t i=0; i<getNumDOF(); i++) {
802            // add the DOF itself
803            index.push_back(i);
804            const dim_t num=insertNeighbourNodes(index, i);
805            ptr.push_back(ptr.back()+num+1);
806        }
807    
808        return createPasoPattern(ptr, index, ptr.size()-1, ptr.size()-1);
809    }
810    
811    //protected
812    void RipleyDomain::createCouplePatterns(const vector<IndexVector>& colIndices,
813            const dim_t N, Paso_Pattern** colPattern, Paso_Pattern** rowPattern) const
814    {
815        IndexVector ptr(1,0);
816        IndexVector index;
817        for (index_t i=0; i<getNumDOF(); i++) {
818            index.insert(index.end(), colIndices[i].begin(), colIndices[i].end());
819            ptr.push_back(ptr.back()+colIndices[i].size());
820        }
821    
822        const dim_t M=ptr.size()-1;
823        *colPattern=createPasoPattern(ptr, index, M, N);
824    
825        IndexVector rowPtr(1,0);
826        IndexVector rowIndex;
827        for (dim_t id=0; id<N; id++) {
828            dim_t n=0;
829            for (dim_t i=0; i<M; i++) {
830                for (dim_t j=ptr[i]; j<ptr[i+1]; j++) {
831                    if (index[j]==id) {
832                        rowIndex.push_back(i);
833                        n++;
834                        break;
835                    }
836                }
837            }
838            rowPtr.push_back(rowPtr.back()+n);
839        }
840    
841        // M and N are now reversed
842        *rowPattern=createPasoPattern(rowPtr, rowIndex, N, M);
843    }
844    
845    //protected
846    void RipleyDomain::addToSystemMatrix(Paso_SystemMatrix* mat,
847           const IndexVector& nodes_Eq, dim_t num_Eq, const IndexVector& nodes_Sol,
848           dim_t num_Sol, const vector<double>& array) const
849    {
850        const dim_t numMyCols = mat->pattern->mainPattern->numInput;
851        const dim_t numMyRows = mat->pattern->mainPattern->numOutput;
852    
853        const index_t* mainBlock_ptr = mat->mainBlock->pattern->ptr;
854        const index_t* mainBlock_index = mat->mainBlock->pattern->index;
855        double* mainBlock_val = mat->mainBlock->val;
856        const index_t* col_coupleBlock_ptr = mat->col_coupleBlock->pattern->ptr;
857        const index_t* col_coupleBlock_index = mat->col_coupleBlock->pattern->index;
858        double* col_coupleBlock_val = mat->col_coupleBlock->val;
859        const index_t* row_coupleBlock_ptr = mat->row_coupleBlock->pattern->ptr;
860        const index_t* row_coupleBlock_index = mat->row_coupleBlock->pattern->index;
861        double* row_coupleBlock_val = mat->row_coupleBlock->val;
862    
863        for (dim_t k_Eq = 0; k_Eq < nodes_Eq.size(); ++k_Eq) {
864            // down columns of array
865            const dim_t j_Eq = nodes_Eq[k_Eq];
866            const dim_t i_row = j_Eq;
867    //printf("row:%d\n", i_row);
868            // only look at the matrix rows stored on this processor
869            if (i_row < numMyRows) {
870                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
871                    const dim_t i_col = nodes_Sol[k_Sol];
872                    if (i_col < numMyCols) {
873                        for (dim_t k = mainBlock_ptr[i_row]; k < mainBlock_ptr[i_row + 1]; ++k) {
874                            if (mainBlock_index[k] == i_col) {
875                                mainBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
876                                break;
877                            }
878                        }
879                    } else {
880                        for (dim_t k = col_coupleBlock_ptr[i_row]; k < col_coupleBlock_ptr[i_row + 1]; ++k) {
881    //cout << "col:" << i_col-numMyCols << " colIdx:" << col_coupleBlock_index[k] << endl;
882                            if (col_coupleBlock_index[k] == i_col - numMyCols) {
883                                col_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
884                                break;
885                            }
886                        }
887                    }
888                }
889            } else {
890                for (dim_t k_Sol = 0; k_Sol < nodes_Sol.size(); ++k_Sol) {
891                    // across rows of array
892                    const dim_t i_col = nodes_Sol[k_Sol];
893                    if (i_col < numMyCols) {
894                        for (dim_t k = row_coupleBlock_ptr[i_row - numMyRows];
895                             k < row_coupleBlock_ptr[i_row - numMyRows + 1]; ++k)
896                        {
897    //cout << "col:" << i_col << " rowIdx:" << row_coupleBlock_index[k] << endl;
898                            if (row_coupleBlock_index[k] == i_col) {
899                                row_coupleBlock_val[k] += array[INDEX2(k_Eq, k_Sol, nodes_Eq.size())];
900                                break;
901                            }
902                        }
903                    }
904                }
905            }
906        }
907    }
908    
909  //  //
910  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
911  //  //
# Line 612  string RipleyDomain::getDescription() co Line 915  string RipleyDomain::getDescription() co
915      throw RipleyException("getDescription() not implemented");      throw RipleyException("getDescription() not implemented");
916  }  }
917    
 bool RipleyDomain::operator==(const AbstractDomain& other) const  
 {  
     throw RipleyException("operator==() not implemented");  
 }  
   
918  void RipleyDomain::write(const string& filename) const  void RipleyDomain::write(const string& filename) const
919  {  {
920      throw RipleyException("write() not implemented");      throw RipleyException("write() not implemented");
# Line 632  const int* RipleyDomain::borrowSampleRef Line 930  const int* RipleyDomain::borrowSampleRef
930      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
931  }  }
932    
 bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  
                                              int fsType_target) const  
 {  
     throw RipleyException("probeInterpolationOnDomain() not implemented");  
 }  
   
933  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const  void RipleyDomain::interpolateACross(escript::Data& target, const escript::Data& source) const
934  {  {
935      throw RipleyException("interpolateACross() not implemented");      throw RipleyException("interpolateACross() not implemented");
# Line 674  bool RipleyDomain::ownSample(int fsType, Line 966  bool RipleyDomain::ownSample(int fsType,
966      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
967  }  }
968    
 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");  
 }  
   
969  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,  void RipleyDomain::addPDEToLumpedSystem(escript::Data& mat,
970          const escript::Data& D, const escript::Data& d,          const escript::Data& D, const escript::Data& d,
971          const escript::Data& d_dirac, const bool useHRZ) const          const escript::Data& d_dirac, const bool useHRZ) const
# Line 769  dim_t RipleyDomain::getNumNodes() const Line 1050  dim_t RipleyDomain::getNumNodes() const
1050      throw RipleyException("getNumNodes() not implemented");      throw RipleyException("getNumNodes() not implemented");
1051  }  }
1052    
1053    dim_t RipleyDomain::getNumDOF() const
1054    {
1055        throw RipleyException("getNumDOF() not implemented");
1056    }
1057    
1058    dim_t RipleyDomain::insertNeighbourNodes(IndexVector& index, index_t node) const
1059    {
1060        throw RipleyException("insertNeighbourNodes() not implemented");
1061    }
1062    
1063  void RipleyDomain::assembleCoordinates(escript::Data& arg) const  void RipleyDomain::assembleCoordinates(escript::Data& arg) const
1064  {  {
1065      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
1066  }  }
1067    
1068    void RipleyDomain::assemblePDESingle(Paso_SystemMatrix* mat, escript::Data& rhs,
1069            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1070            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1071            const escript::Data& d, const escript::Data& y) const
1072    {
1073        throw RipleyException("assemblePDESingle() not implemented");
1074    }
1075    
1076    void RipleyDomain::assemblePDESystem(Paso_SystemMatrix* mat, escript::Data& rhs,
1077            const escript::Data& A, const escript::Data& B, const escript::Data& C,
1078            const escript::Data& D, const escript::Data& X, const escript::Data& Y,
1079            const escript::Data& d, const escript::Data& y) const
1080    {
1081        throw RipleyException("assemblePDESystem() not implemented");
1082    }
1083    
1084  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const  void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
1085  {  {
1086      throw RipleyException("interpolateNodesOnElements() not implemented");      throw RipleyException("interpolateNodesOnElements() not implemented");
# Line 784  void RipleyDomain::interpolateNodesOnFac Line 1091  void RipleyDomain::interpolateNodesOnFac
1091      throw RipleyException("interpolateNodesOnFaces() not implemented");      throw RipleyException("interpolateNodesOnFaces() not implemented");
1092  }  }
1093    
1094    void RipleyDomain::nodesToDOF(escript::Data& out, escript::Data& in) const
1095    {
1096        throw RipleyException("nodesToDOF() not implemented");
1097    }
1098    
1099    void RipleyDomain::dofToNodes(escript::Data& out, escript::Data& in) const
1100    {
1101        throw RipleyException("dofToNodes() not implemented");
1102    }
1103    
1104  } // end of namespace ripley  } // end of namespace ripley
1105    

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

  ViewVC Help
Powered by ViewVC 1.1.26