/[escript]/trunk/speckley/src/SpeckleyDomain.cpp
ViewVC logotype

Diff of /trunk/speckley/src/SpeckleyDomain.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 5477 by jfenwick, Sat Feb 14 00:25:03 2015 UTC revision 5478 by sshaw, Wed Feb 18 01:57:49 2015 UTC
# Line 80  bool SpeckleyDomain::isValidFunctionSpac Line 80  bool SpeckleyDomain::isValidFunctionSpac
80          case DegreesOfFreedom:          case DegreesOfFreedom:
81          case Nodes:          case Nodes:
82          case Elements:          case Elements:
83            case ReducedElements:
84          case Points:          case Points:
85              return true;              return true;
86          default:          default:
# Line 118  pair<int,dim_t> SpeckleyDomain::getDataS Line 119  pair<int,dim_t> SpeckleyDomain::getDataS
119              return pair<int,dim_t>(1, getNumDOF());              return pair<int,dim_t>(1, getNumDOF());
120          case Elements:          case Elements:
121              return pair<int,dim_t>(ptsPerSample, getNumElements());              return pair<int,dim_t>(ptsPerSample, getNumElements());
122            case ReducedElements:
123                return pair<int,dim_t>(1, getNumElements());
124          case Points:          case Points:
125              return pair<int,dim_t>(1, m_diracPoints.size());              return pair<int,dim_t>(1, m_diracPoints.size());
126          default:          default:
# Line 147  bool SpeckleyDomain::commonFunctionSpace Line 150  bool SpeckleyDomain::commonFunctionSpace
150      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
151      interpolated back and forth):      interpolated back and forth):
152      class 0: DOF <-> Nodes      class 0: DOF <-> Nodes
153        class 1: ReducedDOF <-> ReducedNodes
154      class 2: Points      class 2: Points
155      class 3: Elements      class 3: Elements
156        class 4: ReducedElements
157    
158      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
159      between lines.      between lines.
160      class 0 and 1 belong to all lines so aren't considered.      class 0 and 1 belong to all lines so aren't considered.
161      line 0: class 2      line 0: class 2
162      line 1: class 3      line 1: class 3,4
163    
164      For classes with multiple members (eg class 1) we have vars to record if      For classes with multiple members (eg class 1) we have vars to record if
165      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
# Line 165  bool SpeckleyDomain::commonFunctionSpace Line 170  bool SpeckleyDomain::commonFunctionSpace
170      vector<bool> hasclass(7, false);      vector<bool> hasclass(7, false);
171      vector<int> hasline(3, 0);      vector<int> hasline(3, 0);
172      bool hasnodes=false;      bool hasnodes=false;
173        bool hasrednodes=false;
174      for (size_t i=0; i<fs.size(); ++i) {      for (size_t i=0; i<fs.size(); ++i) {
175          switch (fs[i]) {          switch (fs[i]) {
176              case Nodes: hasnodes=true; // fall through              case Nodes: hasnodes=true; // fall through
177              case DegreesOfFreedom:              case DegreesOfFreedom:
178                  hasclass[0]=true;                  hasclass[0]=true;
179                  break;                  break;
180                case ReducedNodes: hasrednodes=true; // fall through
181                case ReducedDegreesOfFreedom:
182                    hasclass[1]=true;
183                    break;
184              case Points:              case Points:
185                  hasline[0]=1;                  hasline[0]=1;
186                  hasclass[2]=true;                  hasclass[2]=true;
# Line 179  bool SpeckleyDomain::commonFunctionSpace Line 189  bool SpeckleyDomain::commonFunctionSpace
189                  hasclass[3]=true;                  hasclass[3]=true;
190                  hasline[1]=1;                  hasline[1]=1;
191                  break;                  break;
192                case ReducedElements:
193                    hasclass[4]=true;
194                    hasline[1]=1;
195                    break;
196              default:              default:
197                  return false;                  return false;
198          }          }
# Line 194  bool SpeckleyDomain::commonFunctionSpace Line 208  bool SpeckleyDomain::commonFunctionSpace
208          if (hasline[0]==1)          if (hasline[0]==1)
209              resultcode=Points;              resultcode=Points;
210          else if (hasline[1]==1) {          else if (hasline[1]==1) {
211                if (hasclass[4])
212                    resultcode=ReducedElements;
213                else
214              resultcode=Elements;              resultcode=Elements;
215          }          }
216      } else { // numLines==0      } else { // numLines==0
217            if (hasclass[1])
218                // something from class 1
219                resultcode=(hasrednodes ? ReducedNodes : ReducedDegreesOfFreedom);
220            else // something from class 0
221          resultcode=(hasnodes ? Nodes : DegreesOfFreedom);          resultcode=(hasnodes ? Nodes : DegreesOfFreedom);
222      }      }
223      return true;      return true;
# Line 216  bool SpeckleyDomain::probeInterpolationO Line 237  bool SpeckleyDomain::probeInterpolationO
237          case Nodes:          case Nodes:
238          case DegreesOfFreedom:          case DegreesOfFreedom:
239              return true;              return true;
240            case ReducedNodes:
241            case ReducedDegreesOfFreedom:
242                return (fsType_target != Nodes &&
243                        fsType_target != DegreesOfFreedom);
244          case Elements:          case Elements:
245              return (fsType_target==Elements || fsType_target==Nodes);              return (fsType_target==Elements || fsType_target==ReducedElements || fsType_target==Nodes);
246          case Points:          case Points:
247              return (fsType_target==fsType_source);              return (fsType_target==fsType_source);
248            case ReducedElements:
249                return (fsType_target==Elements || fsType_target==Nodes);
250    
251          default: {          default: {
252              stringstream msg;              stringstream msg;
# Line 255  signed char SpeckleyDomain::preferredInt Line 282  signed char SpeckleyDomain::preferredInt
282          case Nodes:          case Nodes:
283          case DegreesOfFreedom:          case DegreesOfFreedom:
284              return 1;              return 1;
285          case Elements:          case ReducedNodes:
286              return 0;          case ReducedDegreesOfFreedom:
287                return (fsType_target != Nodes &&
288                        fsType_target != DegreesOfFreedom) ? -1 : 0;
289            case Elements:
290                return (fsType_target==ReducedElements) ? -1 : 0;
291            case ReducedElements:
292                return (fsType_target==Elements) ? 1 : 0;
293          case Points:          case Points:
294              return 0;  // other case caught by the if above              return 0;  // other case caught by the if above
295          default: {          default: {
# Line 290  void SpeckleyDomain::interpolateOnDomain Line 323  void SpeckleyDomain::interpolateOnDomain
323      // simplest case: 1:1 copy      // simplest case: 1:1 copy
324      if (inFS==outFS) {      if (inFS==outFS) {
325          copyData(target, in);          copyData(target, in);
326      } else if (inFS==Elements && outFS==Nodes) {      } else if ((inFS==Elements || inFS==ReducedElements) && outFS==Nodes) {
327              interpolateElementsOnNodes(target, in);          interpolateElementsOnNodes(target, in);
328        } else if (inFS==ReducedElements && outFS==Elements) {
329            multiplyData(target, in);
330        } else if (inFS==Elements && outFS==ReducedElements) {
331            reduceElements(target, in);
332      } else {      } else {
333          switch (inFS) {          switch (inFS) {
334              case Nodes:              case Nodes:
# Line 301  void SpeckleyDomain::interpolateOnDomain Line 338  void SpeckleyDomain::interpolateOnDomain
338                          copyData(target, in);                          copyData(target, in);
339                          break;                          break;
340                      case Elements:                      case Elements:
341                          interpolateNodesOnElements(target, in);                          interpolateNodesOnElements(target, in, false);
342                            break;
343                        case ReducedElements:
344                            interpolateNodesOnElements(target, in, true);
345                          break;                          break;
   
346                      case Points:                      case Points:
347                          {                          {
348                              const dim_t numComp = in.getDataPointSize();                              const dim_t numComp = in.getDataPointSize();
# Line 329  void SpeckleyDomain::interpolateOnDomain Line 368  void SpeckleyDomain::interpolateOnDomain
368                          break;                          break;
369    
370                      case Elements:                      case Elements:
371                        case ReducedElements:
372                          if (getMPISize()==1) {                          if (getMPISize()==1) {
373                              interpolateNodesOnElements(target, in);                              interpolateNodesOnElements(target, in, outFS==ReducedElements);
374                          } else {                          } else {
375                              escript::Data contIn(in, escript::continuousFunction(*this));                              escript::Data contIn(in, escript::continuousFunction(*this));
376                              interpolateNodesOnElements(target, contIn);                              interpolateNodesOnElements(target, contIn, outFS==ReducedElements);
377                          }                          }
378                          break;                          break;
379    
# Line 411  void SpeckleyDomain::setToGradient(escri Line 451  void SpeckleyDomain::setToGradient(escri
451      switch (grad.getFunctionSpace().getTypeCode()) {      switch (grad.getFunctionSpace().getTypeCode()) {
452          case Nodes:          case Nodes:
453          case Elements:          case Elements:
454            case ReducedElements:
455              break;              break;
456          default: {          default: {
457              stringstream msg;              stringstream msg;
# Line 458  void SpeckleyDomain::setToIntegrals(vect Line 499  void SpeckleyDomain::setToIntegrals(vect
499              }              }
500              break;              break;
501          case Elements:          case Elements:
502            case ReducedElements:
503              assembleIntegrate(integrals, arg);              assembleIntegrate(integrals, arg);
504              break;              break;
505          default: {          default: {
# Line 477  bool SpeckleyDomain::isCellOriented(int Line 519  bool SpeckleyDomain::isCellOriented(int
519          case DegreesOfFreedom:          case DegreesOfFreedom:
520              return false;              return false;
521          case Elements:          case Elements:
522            case ReducedElements:
523          case Points:          case Points:
524              return true;              return true;
525          default:          default:
# Line 493  bool SpeckleyDomain::canTag(int fsType) Line 536  bool SpeckleyDomain::canTag(int fsType)
536      switch(fsType) {      switch(fsType) {
537          case Nodes:          case Nodes:
538          case Elements:          case Elements:
539            case ReducedElements:
540          case Points:          case Points:
541              return true;              return true;
542          case DegreesOfFreedom:          case DegreesOfFreedom:
# Line 547  int SpeckleyDomain::getTagFromSampleNo(i Line 591  int SpeckleyDomain::getTagFromSampleNo(i
591                  return m_nodeTags[sampleNo];                  return m_nodeTags[sampleNo];
592              break;              break;
593          case Elements:          case Elements:
594            case ReducedElements:
595              if (m_elementTags.size() > sampleNo)              if (m_elementTags.size() > sampleNo)
596                  return m_elementTags[sampleNo];                  return m_elementTags[sampleNo];
597              break;              break;
# Line 569  int SpeckleyDomain::getNumberOfTagsInUse Line 614  int SpeckleyDomain::getNumberOfTagsInUse
614          case Nodes:          case Nodes:
615              return m_nodeTagsInUse.size();              return m_nodeTagsInUse.size();
616          case Elements:          case Elements:
617            case ReducedElements:
618              return m_elementTagsInUse.size();              return m_elementTagsInUse.size();
619          default: {          default: {
620              stringstream msg;              stringstream msg;
# Line 585  const int* SpeckleyDomain::borrowListOfT Line 631  const int* SpeckleyDomain::borrowListOfT
631          case Nodes:          case Nodes:
632              return &m_nodeTagsInUse[0];              return &m_nodeTagsInUse[0];
633          case Elements:          case Elements:
634            case ReducedElements:
635              return &m_elementTagsInUse[0];              return &m_elementTagsInUse[0];
636          default: {          default: {
637              stringstream msg;              stringstream msg;
# Line 724  void SpeckleyDomain::copyData(escript::D Line 771  void SpeckleyDomain::copyData(escript::D
771  }  }
772    
773  //protected  //protected
774    void SpeckleyDomain::multiplyData(escript::Data& out, const escript::Data& in) const
775    {
776        const dim_t numComp = in.getDataPointSize();
777        const dim_t dpp = out.getNumDataPointsPerSample();
778        const dim_t numSamples = in.getNumSamples();
779        out.requireWrite();
780    #pragma omp parallel for
781        for (index_t i=0; i<numSamples; i++) {
782            const double* src = in.getSampleDataRO(i);
783            double* dest = out.getSampleDataRW(i);
784            for (index_t c=0; c<numComp; c++) {
785                for (index_t q=0; q<dpp; q++)
786                    dest[c+q*numComp] = src[c];
787            }
788        }
789    }
790    
791    //protected
792  void SpeckleyDomain::updateTagsInUse(int fsType) const  void SpeckleyDomain::updateTagsInUse(int fsType) const
793  {  {
794      std::vector<int>* tagsInUse=NULL;      std::vector<int>* tagsInUse=NULL;
# Line 748  void SpeckleyDomain::updateTagsInUse(int Line 813  void SpeckleyDomain::updateTagsInUse(int
813      tagsInUse->clear();      tagsInUse->clear();
814      int lastFoundValue = numeric_limits<int>::min();      int lastFoundValue = numeric_limits<int>::min();
815      int minFoundValue, local_minFoundValue;      int minFoundValue, local_minFoundValue;
816      const long numTags = tags->size();      const int numTags = tags->size();
817    
818      while (true) {      while (true) {
819          // find smallest value bigger than lastFoundValue          // find smallest value bigger than lastFoundValue
# Line 756  void SpeckleyDomain::updateTagsInUse(int Line 821  void SpeckleyDomain::updateTagsInUse(int
821  #pragma omp parallel private(local_minFoundValue)  #pragma omp parallel private(local_minFoundValue)
822          {          {
823              local_minFoundValue = minFoundValue;              local_minFoundValue = minFoundValue;
824              long i;     // should be size_t but omp mutter mutter  #pragma omp for schedule(static) nowait
825  #pragma omp for schedule(static) private(i) nowait              for (int i = 0; i < numTags; i++) {
             for (i = 0; i < numTags; i++) {  
826                  const int v = (*tags)[i];                  const int v = (*tags)[i];
827                  if ((v > lastFoundValue) && (v < local_minFoundValue))                  if ((v > lastFoundValue) && (v < local_minFoundValue))
828                      local_minFoundValue = v;                      local_minFoundValue = v;

Legend:
Removed from v.5477  
changed lines
  Added in v.5478

  ViewVC Help
Powered by ViewVC 1.1.26