/[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 3697 by caltinay, Mon Nov 28 04:52:00 2011 UTC revision 3722 by caltinay, Wed Dec 7 05:53:22 2011 UTC
# Line 89  pair<int,int> RipleyDomain::getDataShape Line 89  pair<int,int> RipleyDomain::getDataShape
89      const int ptsPerSample = (m_numDim==2 ? 4 : 8);      const int ptsPerSample = (m_numDim==2 ? 4 : 8);
90      switch (fsType) {      switch (fsType) {
91          case Nodes:          case Nodes:
92            case ReducedNodes:
93          case DegreesOfFreedom:          case DegreesOfFreedom:
94            case ReducedDegreesOfFreedom:
95              return pair<int,int>(1, getNumNodes());              return pair<int,int>(1, getNumNodes());
96          case Elements:          case Elements:
97              return pair<int,int>(ptsPerSample, getNumElements());              return pair<int,int>(ptsPerSample, getNumElements());
98          case FaceElements:          case FaceElements:
99              return pair<int,int>(ptsPerSample/2, getNumFaceElements());              return pair<int,int>(ptsPerSample/2, getNumFaceElements());
             /*  
         case Points:  
             return pair<int,int>(1, getNumNodes());  
         case ReducedNodes:  
         case ReducedDegreesOfFreedom:  
             return pair<int,int>(1, getNumReducedNodes());  
100          case ReducedElements:          case ReducedElements:
101              return pair<int,int>(1, getNumReducedElements());              return pair<int,int>(1, getNumElements());
102          case ReducedFaceElements:          case ReducedFaceElements:
103              return pair<int,int>(1, getNumReducedFaceElements());              return pair<int,int>(1, getNumFaceElements());
104                /*
105            case Points:
106                return pair<int,int>(1, getNumPoints());
107              */              */
108          default:          default:
109              break;              break;
110      }      }
111    
112      stringstream msg;      stringstream msg;
113      msg << "Invalid function space type " << fsType << " for "      msg << "getDataShape(): Unsupported function space type "
114          << getDescription();          << functionSpaceTypeAsString(fsType) << " for " << getDescription();
115      throw RipleyException(msg.str());      throw RipleyException(msg.str());
116  }  }
117    
 int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const  
 {  
     throw RipleyException("getTagFromSampleNo() not implemented");  
 }  
   
118  string RipleyDomain::showTagNames() const  string RipleyDomain::showTagNames() const
119  {  {
120      stringstream ret;      stringstream ret;
# Line 231  bool RipleyDomain::commonFunctionSpace(c Line 225  bool RipleyDomain::commonFunctionSpace(c
225      return true;      return true;
226  }  }
227    
228    void RipleyDomain::interpolateOnDomain(escript::Data& target,
229                                           const escript::Data& in) const
230    {
231        const RipleyDomain& inDomain=dynamic_cast<const RipleyDomain&>(*(in.getFunctionSpace().getDomain()));
232        const RipleyDomain& targetDomain=dynamic_cast<const RipleyDomain&>(*(target.getFunctionSpace().getDomain()));
233        if (inDomain != *this)
234            throw RipleyException("Illegal domain of interpolant");
235        if (targetDomain != *this)
236            throw RipleyException("Illegal domain of interpolation target");
237    
238        stringstream msg;
239        msg << "interpolateOnDomain() not implemented for function space "
240            << functionSpaceTypeAsString(in.getFunctionSpace().getTypeCode())
241            << " -> "
242            << functionSpaceTypeAsString(target.getFunctionSpace().getTypeCode());
243    
244        switch (in.getFunctionSpace().getTypeCode()) {
245            case Nodes:
246            case ReducedNodes:
247            case DegreesOfFreedom:
248            case ReducedDegreesOfFreedom:
249                switch (target.getFunctionSpace().getTypeCode()) {
250                    case Nodes:
251                    case ReducedNodes:
252                    case DegreesOfFreedom:
253                    case ReducedDegreesOfFreedom:
254                        // FIXME
255                        copyNodalData(target, *const_cast<escript::Data*>(&in));
256                        break;
257    
258                    case Elements:
259                        interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), false);
260                        break;
261    
262                    case ReducedElements:
263                        interpolateNodesOnElements(target, *const_cast<escript::Data*>(&in), true);
264                        break;
265    
266                    case FaceElements:
267                        interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), false);
268                        break;
269    
270                    case ReducedFaceElements:
271                        interpolateNodesOnFaces(target, *const_cast<escript::Data*>(&in), true);
272                        break;
273    
274                    default:
275                        throw RipleyException(msg.str());
276                }
277                break;
278            default:
279                throw RipleyException(msg.str());
280        }
281    }
282    
283  escript::Data RipleyDomain::getX() const  escript::Data RipleyDomain::getX() const
284  {  {
285      return escript::continuousFunction(*this).getX();      return escript::continuousFunction(*this).getX();
# Line 282  bool RipleyDomain::isCellOriented(int fs Line 331  bool RipleyDomain::isCellOriented(int fs
331              break;              break;
332      }      }
333      stringstream msg;      stringstream msg;
334      msg << "Illegal function space type " << fsType << " on " << getDescription();      msg << "isCellOriented(): Illegal function space type " << fsType << " on " << getDescription();
335        throw RipleyException(msg.str());
336    }
337    
338    bool RipleyDomain::canTag(int fsType) const
339    {
340        switch(fsType) {
341            case Nodes:
342            case Elements:
343            case ReducedElements:
344            case FaceElements:
345            case ReducedFaceElements:
346                return true;
347            case DegreesOfFreedom:
348            case ReducedDegreesOfFreedom:
349            case Points:
350                return false;
351            default:
352                break;
353        }
354        stringstream msg;
355        msg << "canTag(): Illegal function space type " << fsType << " on " << getDescription();
356      throw RipleyException(msg.str());      throw RipleyException(msg.str());
357  }  }
358    
359    void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& cMask) const
360    {
361        IndexVector* target=NULL;
362        dim_t num=0;
363    
364        switch(fsType) {
365            case Nodes:
366                num=getNumNodes();
367                target=&m_nodeTags;
368                break;
369            case Elements:
370            case ReducedElements:
371                num=getNumElements();
372                target=&m_elementTags;
373                break;
374            case FaceElements:
375            case ReducedFaceElements:
376                num=getNumFaceElements();
377                target=&m_faceTags;
378                break;
379            default: {
380                stringstream msg;
381                msg << "setTags(): not implemented for " << functionSpaceTypeAsString(fsType);
382                throw RipleyException(msg.str());
383            }
384        }
385        if (target->size() != num) {
386            target->assign(num, -1);
387        }
388    
389        escript::Data& mask=*const_cast<escript::Data*>(&cMask);
390    #pragma omp parallel for
391        for (index_t i=0; i<num; i++) {
392            if (mask.getSampleDataRO(i)[0] > 0) {
393                (*target)[i]=newTag;
394            }
395        }
396        updateTagsInUse(fsType);
397    }
398    
399    int RipleyDomain::getTagFromSampleNo(int fsType, int sampleNo) const
400    {
401        switch(fsType) {
402            case Nodes:
403                if (m_nodeTags.size() > sampleNo)
404                    return m_nodeTags[sampleNo];
405                break;
406            case Elements:
407            case ReducedElements:
408                if (m_elementTags.size() > sampleNo)
409                    return m_elementTags[sampleNo];
410                break;
411            case FaceElements:
412            case ReducedFaceElements:
413                if (m_faceTags.size() > sampleNo)
414                    return m_faceTags[sampleNo];
415                break;
416            default: {
417                stringstream msg;
418                msg << "getTagFromSampleNo(): not implemented for " << functionSpaceTypeAsString(fsType);
419                throw RipleyException(msg.str());
420            }
421        }
422        return -1;
423    }
424    
425    int RipleyDomain::getNumberOfTagsInUse(int fsType) const
426    {
427        switch(fsType) {
428            case Nodes:
429                return m_nodeTagsInUse.size();
430            case Elements:
431            case ReducedElements:
432                return m_elementTagsInUse.size();
433            case FaceElements:
434            case ReducedFaceElements:
435                return m_faceTagsInUse.size();
436            default: {
437                stringstream msg;
438                msg << "getNumberOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);
439                throw RipleyException(msg.str());
440            }
441        }
442    }
443    
444    const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const
445    {
446        switch(fsType) {
447            case Nodes:
448                return &m_nodeTagsInUse[0];
449            case Elements:
450            case ReducedElements:
451                return &m_elementTagsInUse[0];
452            case FaceElements:
453            case ReducedFaceElements:
454                return &m_faceTagsInUse[0];
455            default: {
456                stringstream msg;
457                msg << "borrowListOfTagsInUse(): not implemented for " << functionSpaceTypeAsString(fsType);
458                throw RipleyException(msg.str());
459            }
460        }
461    }
462    
463  void RipleyDomain::Print_Mesh_Info(const bool full) const  void RipleyDomain::Print_Mesh_Info(const bool full) const
464  {  {
465      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "      cout << "Print_Mesh_Info for " << getDescription() << " running on CPU "
# Line 358  void RipleyDomain::setNewX(const escript Line 532  void RipleyDomain::setNewX(const escript
532      throw RipleyException("setNewX(): Operation not supported");      throw RipleyException("setNewX(): Operation not supported");
533  }  }
534    
535    //protected
536    void RipleyDomain::copyNodalData(escript::Data& out, escript::Data& in) const
537    {
538        const dim_t numComp = in.getDataPointSize();
539    #pragma omp parallel for
540        for (index_t i=0; i<in.getNumSamples(); i++) {
541            const double* src = in.getSampleDataRO(i);
542            copy(src, src+numComp, out.getSampleDataRW(i));
543        }
544    }
545    
546    //protected
547    void RipleyDomain::updateTagsInUse(int fsType) const
548    {
549        IndexVector* tagsInUse=NULL;
550        const IndexVector* tags=NULL;
551        switch(fsType) {
552            case Nodes:
553                tags=&m_nodeTags;
554                tagsInUse=&m_nodeTagsInUse;
555                break;
556            case Elements:
557            case ReducedElements:
558                tags=&m_elementTags;
559                tagsInUse=&m_elementTagsInUse;
560                break;
561            case FaceElements:
562            case ReducedFaceElements:
563                tags=&m_faceTags;
564                tagsInUse=&m_faceTagsInUse;
565                break;
566            default:
567                return;
568        }
569    
570        // gather global unique tag values from tags into tagsInUse
571        tagsInUse->clear();
572        index_t lastFoundValue = INDEX_T_MIN, minFoundValue, local_minFoundValue;
573    
574        while (true) {
575            // find smallest value bigger than lastFoundValue
576            minFoundValue = INDEX_T_MAX;
577    #pragma omp parallel private(local_minFoundValue)
578            {
579                local_minFoundValue = minFoundValue;
580    #pragma omp for schedule(static) nowait
581                for (size_t i = 0; i < tags->size(); i++) {
582                    const index_t v = (*tags)[i];
583                    if ((v > lastFoundValue) && (v < local_minFoundValue))
584                        local_minFoundValue = v;
585                }
586    #pragma omp critical
587                {
588                    if (local_minFoundValue < minFoundValue)
589                        minFoundValue = local_minFoundValue;
590                }
591            }
592    #ifdef ESYS_MPI
593            local_minFoundValue = minFoundValue;
594            MPI_Allreduce(&local_minFoundValue, &minFoundValue, 1, MPI_INT, MPI_MIN, m_mpiInfo->comm);
595    #endif
596            // if we found a new value add it to the tagsInUse vector
597    
598            if (minFoundValue < INDEX_T_MAX) {
599                tagsInUse->push_back(minFoundValue);
600                lastFoundValue = minFoundValue;
601            } else
602                break;
603        }
604    }
605    
606  //  //
607  // the methods that follow have to be implemented by the subclasses  // the methods that follow have to be implemented by the subclasses
608  //  //
# Line 387  const int* RipleyDomain::borrowSampleRef Line 632  const int* RipleyDomain::borrowSampleRef
632      throw RipleyException("borrowSampleReferenceIDs() not implemented");      throw RipleyException("borrowSampleReferenceIDs() not implemented");
633  }  }
634    
 void RipleyDomain::interpolateOnDomain(escript::Data& target,  
                                        const escript::Data& in) const  
 {  
     throw RipleyException("interpolateOnDomain() not implemented");  
 }  
   
635  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,  bool RipleyDomain::probeInterpolationOnDomain(int fsType_source,
636                                               int fsType_target) const                                               int fsType_target) const
637  {  {
# Line 435  bool RipleyDomain::ownSample(int fsType, Line 674  bool RipleyDomain::ownSample(int fsType,
674      throw RipleyException("ownSample() not implemented");      throw RipleyException("ownSample() not implemented");
675  }  }
676    
 void RipleyDomain::setTags(const int fsType, const int newTag, const escript::Data& mask) const  
 {  
     throw RipleyException("setTags() not implemented");  
 }  
   
 int RipleyDomain::getNumberOfTagsInUse(int fsType) const  
 {  
     throw RipleyException("getNumberOfTagsInUse() not implemented");  
 }  
   
 const int* RipleyDomain::borrowListOfTagsInUse(int fsType) const  
 {  
     throw RipleyException("borrowListOfTagsInUse() not implemented");  
 }  
   
 bool RipleyDomain::canTag(int fsType) const  
 {  
     throw RipleyException("canTag() not implemented");  
 }  
   
   
677  void RipleyDomain::addPDEToSystem(  void RipleyDomain::addPDEToSystem(
678          escript::AbstractSystemMatrix& mat, escript::Data& rhs,          escript::AbstractSystemMatrix& mat, escript::Data& rhs,
679          const escript::Data& A, const escript::Data& B, const escript::Data& C,          const escript::Data& A, const escript::Data& B, const escript::Data& C,
# Line 556  void RipleyDomain::assembleCoordinates(e Line 774  void RipleyDomain::assembleCoordinates(e
774      throw RipleyException("assembleCoordinates() not implemented");      throw RipleyException("assembleCoordinates() not implemented");
775  }  }
776    
777    void RipleyDomain::interpolateNodesOnElements(escript::Data& out, escript::Data& in, bool reduced) const
778    {
779        throw RipleyException("interpolateNodesOnElements() not implemented");
780    }
781    
782    void RipleyDomain::interpolateNodesOnFaces(escript::Data& out, escript::Data& in, bool reduced) const
783    {
784        throw RipleyException("interpolateNodesOnFaces() not implemented");
785    }
786    
787    
788  } // end of namespace ripley  } // end of namespace ripley
789    

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

  ViewVC Help
Powered by ViewVC 1.1.26