/[escript]/trunk/weipa/src/DataVar.cpp
ViewVC logotype

Diff of /trunk/weipa/src/DataVar.cpp

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

revision 2834 by caltinay, Thu Jan 7 06:06:56 2010 UTC revision 2880 by caltinay, Thu Jan 28 01:21:30 2010 UTC
# Line 47  DataVar::DataVar(const string& name) : Line 47  DataVar::DataVar(const string& name) :
47  // Copy constructor  // Copy constructor
48  //  //
49  DataVar::DataVar(const DataVar& d) :  DataVar::DataVar(const DataVar& d) :
50        initialized(d.initialized), finleyMesh(d.finleyMesh),
51      varName(d.varName), numSamples(d.numSamples),      varName(d.varName), numSamples(d.numSamples),
52      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),      rank(d.rank), ptsPerSample(d.ptsPerSample), centering(d.centering),
53      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)      funcSpace(d.funcSpace), shape(d.shape), sampleID(d.sampleID)
# Line 59  DataVar::DataVar(const DataVar& d) : Line 60  DataVar::DataVar(const DataVar& d) :
60              dataArray.push_back(c);              dataArray.push_back(c);
61          }          }
62      }      }
     initialized = d.initialized;  
63  }  }
64    
65  //  //
# Line 97  bool DataVar::initFromEscript(escript::D Line 97  bool DataVar::initFromEscript(escript::D
97          return false;          return false;
98      }      }
99    
100        finleyMesh = mesh;
101      rank = escriptData.getDataPointRank();      rank = escriptData.getDataPointRank();
102      ptsPerSample = escriptData.getNumDataPointsPerSample();      ptsPerSample = escriptData.getNumDataPointsPerSample();
103      shape = escriptData.getDataPointShape();      shape = escriptData.getDataPointShape();
# Line 114  bool DataVar::initFromEscript(escript::D Line 115  bool DataVar::initFromEscript(escript::D
115          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
116  #endif  #endif
117    
118        NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
119        if (nodes == NULL)
120            return false;
121    
122        meshName = nodes->getName();
123        siloMeshName = nodes->getFullSiloName();
124      initialized = true;      initialized = true;
125    
126        // no samples? Nothing more to do.
127      if (numSamples == 0)      if (numSamples == 0)
128          return true;          return true;
129    
# Line 151  bool DataVar::initFromEscript(escript::D Line 159  bool DataVar::initFromEscript(escript::D
159          }          }
160          delete[] tempData;          delete[] tempData;
161    
162          initialized = filterSamples(mesh);          initialized = reorderSamples();
163      }      }
164    
165      return initialized;      return initialized;
# Line 164  bool DataVar::initFromMesh(FinleyMesh_pt Line 172  bool DataVar::initFromMesh(FinleyMesh_pt
172  {  {
173      cleanup();      cleanup();
174            
175      const IntVec& data = mesh->getVarDataByName(varName);      finleyMesh = mesh;
176      rank = 0;      rank = 0;
177      ptsPerSample = 1;      ptsPerSample = 1;
178        NodeData_ptr nodes;
179    
180        if (varName.find("ContactElements_") != varName.npos) {
181            funcSpace = FINLEY_CONTACT_ELEMENTS_1;
182            centering = ZONE_CENTERED;
183            string elementName = varName.substr(0, varName.find('_'));
184            ElementData_ptr elements = mesh->getElementsByName(elementName);
185            nodes = elements->getNodeMesh();
186            sampleID = elements->getIDs();
187        } else if (varName.find("FaceElements_") != varName.npos) {
188            funcSpace = FINLEY_FACE_ELEMENTS;
189            centering = ZONE_CENTERED;
190            string elementName = varName.substr(0, varName.find('_'));
191            ElementData_ptr elements = mesh->getElementsByName(elementName);
192            nodes = elements->getNodeMesh();
193            sampleID = elements->getIDs();
194        } else if (varName.find("Elements_") != varName.npos) {
195            funcSpace = FINLEY_ELEMENTS;
196            centering = ZONE_CENTERED;
197            string elementName = varName.substr(0, varName.find('_'));
198            ElementData_ptr elements = mesh->getElementsByName(elementName);
199            nodes = elements->getNodeMesh();
200            sampleID = elements->getIDs();
201        } else if (varName.find("Nodes_") != varName.npos) {
202            funcSpace = FINLEY_NODES;
203            centering = NODE_CENTERED;
204            nodes = mesh->getNodes();
205            sampleID = nodes->getNodeIDs();
206        } else {
207            cerr << "WARNING: Unrecognized mesh variable '" << varName << "'\n";
208            return false;
209        }
210    
211        meshName = nodes->getName();
212        siloMeshName = nodes->getFullSiloName();
213    
214        const IntVec& data = mesh->getVarDataByName(varName);
215      numSamples = data.size();      numSamples = data.size();
216    
217      if (numSamples > 0) {      if (numSamples > 0) {
# Line 175  bool DataVar::initFromMesh(FinleyMesh_pt Line 220  bool DataVar::initFromMesh(FinleyMesh_pt
220          IntVec::const_iterator it;          IntVec::const_iterator it;
221          for (it=data.begin(); it != data.end(); it++)          for (it=data.begin(); it != data.end(); it++)
222              *c++ = static_cast<float>(*it);              *c++ = static_cast<float>(*it);
   
         if (varName.find("ContactElements_") != varName.npos) {  
             funcSpace = FINLEY_CONTACT_ELEMENTS_1;  
             centering = ZONE_CENTERED;  
             sampleID.insert(sampleID.end(),  
                     mesh->getContactElements()->getIDs().begin(),  
                     mesh->getContactElements()->getIDs().end());  
         } else if (varName.find("FaceElements_") != varName.npos) {  
             funcSpace = FINLEY_FACE_ELEMENTS;  
             centering = ZONE_CENTERED;  
             sampleID.insert(sampleID.end(),  
                     mesh->getFaceElements()->getIDs().begin(),  
                     mesh->getFaceElements()->getIDs().end());  
         } else if (varName.find("Elements_") != varName.npos) {  
             funcSpace = FINLEY_ELEMENTS;  
             centering = ZONE_CENTERED;  
             sampleID.insert(sampleID.end(),  
                     mesh->getElements()->getIDs().begin(),  
                     mesh->getElements()->getIDs().end());  
         } else if (varName.find("Nodes_") != varName.npos) {  
             funcSpace = FINLEY_NODES;  
             centering = NODE_CENTERED;  
             sampleID.insert(sampleID.end(),  
                     mesh->getNodes()->getNodeIDs().begin(),  
                     mesh->getNodes()->getNodeIDs().end());  
         } else {  
             return false;  
         }  
   
         NodeData_ptr nodes = mesh->getMeshForFinleyFS(funcSpace);  
         meshName = nodes->getName();  
         siloMeshName = nodes->getFullSiloName();  
223      }      }
224      initialized = true;      initialized = true;
225    
# Line 263  bool DataVar::initFromNetCDF(const strin Line 276  bool DataVar::initFromNetCDF(const strin
276          << ptsPerSample << " pts/s,  rank: " << rank << endl;          << ptsPerSample << " pts/s,  rank: " << rank << endl;
277  #endif  #endif
278    
279      initialized = true;      finleyMesh = mesh;
280        NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
281      // if there are no data samples we're done      if (nodes == NULL) {
     if (numSamples == 0) {  
282          delete input;          delete input;
283          return true;          return false;
284      }      }
285    
286      sampleID.insert(sampleID.end(), numSamples, 0);      meshName = nodes->getName();
287      NcVar* var = input->get_var("id");      siloMeshName = nodes->getFullSiloName();
288      var->get(&sampleID[0], numSamples);      initialized = true;
289    
290      size_t dimSize = 1;      size_t dimSize = 1;
291      vector<long> counts;      vector<long> counts;
# Line 297  bool DataVar::initFromNetCDF(const strin Line 309  bool DataVar::initFromNetCDF(const strin
309          initialized = false;          initialized = false;
310      }      }
311    
312      if (initialized) {      if (initialized && numSamples > 0) {
313            sampleID.insert(sampleID.end(), numSamples, 0);
314            NcVar* var = input->get_var("id");
315            var->get(&sampleID[0], numSamples);
316    
317          size_t dataSize = dimSize*numSamples*ptsPerSample;          size_t dataSize = dimSize*numSamples*ptsPerSample;
318          counts.push_back(ptsPerSample);          counts.push_back(ptsPerSample);
319          counts.push_back(numSamples);          counts.push_back(numSamples);
320          float* tempData = new float[dataSize];          float* tempData = new float[dataSize];
321          NcVar* var = input->get_var("data");          var = input->get_var("data");
322          var->get(tempData, &counts[0]);          var->get(tempData, &counts[0]);
323    
324          const float* srcPtr = tempData;          const float* srcPtr = tempData;
# Line 312  bool DataVar::initFromNetCDF(const strin Line 328  bool DataVar::initFromNetCDF(const strin
328          }          }
329          delete[] tempData;          delete[] tempData;
330    
331          initialized = filterSamples(mesh);          initialized = reorderSamples();
332      }      }
333    
334      delete input;      delete input;
# Line 337  bool DataVar::isNodeCentered() const Line 353  bool DataVar::isNodeCentered() const
353  //  //
354  float* DataVar::averageData(const float* src, size_t stride)  float* DataVar::averageData(const float* src, size_t stride)
355  {  {
356      float* res = new float[numSamples];      float* res;
357    
358      if (ptsPerSample == 1) {      if (ptsPerSample == 1) {
359            res = new float[numSamples];
360          float* dest = res;          float* dest = res;
361          for (int i=0; i<numSamples; i++, src+=stride)          for (int i=0; i<numSamples; i++, src+=stride)
362              *dest++ = *src;              *dest++ = *src;
363      } else {      } else {
364            ElementData_ptr cells = finleyMesh->getElementsForFinleyFS(funcSpace);
365            int cellFactor = cells->getElementFactor();
366            res = new float[cellFactor * numSamples];
367          float* dest = res;          float* dest = res;
368          for (int i=0; i<numSamples; i++) {          QuadMaskInfo qmi = cells->getQuadMask(funcSpace);
369              double tmpVal = 0.0;          if (qmi.mask.size() > 0) {
370              for (int j=0; j<ptsPerSample; j++, src+=stride)              const float* tmpSrc = src;
371                  tmpVal += *src;              for (int i=0; i<numSamples; i++, tmpSrc+=stride*ptsPerSample) {
372              *dest++ = (float)(tmpVal / ptsPerSample);                  for (int l=0; l<cellFactor; l++) {
373                        double tmpVal = 0.0;
374                        for (int j=0; j<ptsPerSample; j++) {
375                            if (qmi.mask[l][j] != 0) {
376                                tmpVal += *(tmpSrc+stride*j);
377                            }
378                        }
379                        *dest++ = (float)(tmpVal / qmi.factor[l]);
380                    }
381                }
382            } else {
383                for (int i=0; i<numSamples; i++) {
384                    double tmpVal = 0.0;
385                    for (int j=0; j<ptsPerSample; j++, src+=stride) {
386                        tmpVal += *src;
387                    }
388                    tmpVal /= ptsPerSample;
389                    for (int l=0; l<cellFactor; l++) {
390                        *dest++ = static_cast<float>(tmpVal);
391                    }
392                }
393          }          }
394      }      }
395      return res;      return res;
396  }  }
397    
398  //  //
399  // Filters and reorders the raw sample values according to the IDs provided  // Filters and reorders the raw sample values according to the node/element
400  // in 'requiredIDs'. This is used to have data arrays ordered according to  // IDs. This is used to have data arrays ordered according to the underlying
401  // the underlying mesh (i.e. DataID[i]==MeshNodeID[i])  // mesh (i.e. DataID[i]==MeshNodeID[i])
402  //  //
403  bool DataVar::filterSamples(FinleyMesh_ptr finleyMesh)  bool DataVar::reorderSamples()
404  {  {
405      if (numSamples == 0)      if (numSamples == 0)
406          return true;          return true;
407    
     IndexMap id2idxMap;  
408      const IntVec* requiredIDs = NULL;      const IntVec* requiredIDs = NULL;
   
     NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);  
     if (nodes == NULL)  
         return false;  
   
409      int requiredNumSamples = 0;      int requiredNumSamples = 0;
410        int cellFactor = 1;
411    
412      if (centering == NODE_CENTERED) {      if (centering == NODE_CENTERED) {
413          id2idxMap = nodes->getIndexMap();          NodeData_ptr nodes = finleyMesh->getMeshForFinleyFS(funcSpace);
414          requiredIDs = &nodes->getNodeIDs();          requiredIDs = &nodes->getNodeIDs();
415          requiredNumSamples = nodes->getNumNodes();          requiredNumSamples = nodes->getNumNodes();
416      } else {      } else {
# Line 383  bool DataVar::filterSamples(FinleyMesh_p Line 418  bool DataVar::filterSamples(FinleyMesh_p
418          if (cells == NULL)          if (cells == NULL)
419              return false;              return false;
420    
         id2idxMap = cells->getIndexMap();  
421          requiredIDs = &cells->getIDs();          requiredIDs = &cells->getIDs();
422          requiredNumSamples = cells->getNumElements();          requiredNumSamples = cells->getNumElements();
423            cellFactor = cells->getElementFactor();
424            if (cellFactor > 1) {
425                numSamples *= cellFactor;
426                // update sample IDs
427                IntVec newSampleID(numSamples);
428                IntVec::const_iterator idIt = sampleID.begin();
429                IntVec::iterator newIDit = newSampleID.begin();
430                for (; idIt != sampleID.end(); idIt++, newIDit+=cellFactor) {
431                    fill(newIDit, newIDit+cellFactor, *idIt);
432                }
433                sampleID.swap(newSampleID);
434            }
435      }      }
436    
437      if (requiredNumSamples > numSamples) {      if (requiredNumSamples > numSamples) {
# Line 394  bool DataVar::filterSamples(FinleyMesh_p Line 440  bool DataVar::filterSamples(FinleyMesh_p
440          return false;          return false;
441      }      }
442    
     meshName = nodes->getName();  
     siloMeshName = nodes->getFullSiloName();  
   
443      IndexMap sampleID2idx = buildIndexMap();      IndexMap sampleID2idx = buildIndexMap();
444      numSamples = requiredNumSamples;      numSamples = requiredNumSamples;
445    
# Line 405  bool DataVar::filterSamples(FinleyMesh_p Line 448  bool DataVar::filterSamples(FinleyMesh_p
448          float* c = new float[numSamples];          float* c = new float[numSamples];
449          const float* src = dataArray[i];          const float* src = dataArray[i];
450          IntVec::const_iterator idIt = requiredIDs->begin();          IntVec::const_iterator idIt = requiredIDs->begin();
451          for (; idIt != requiredIDs->end(); idIt++) {          size_t destIdx = 0;
452            for (; idIt != requiredIDs->end(); idIt+=cellFactor, destIdx+=cellFactor) {
453              size_t srcIdx = sampleID2idx.find(*idIt)->second;              size_t srcIdx = sampleID2idx.find(*idIt)->second;
454              size_t destIdx = id2idxMap.find(*idIt)->second;              copy(&src[srcIdx], &src[srcIdx+cellFactor], &c[destIdx]);
             c[destIdx] = src[srcIdx];  
455          }          }
456          delete[] dataArray[i];          delete[] dataArray[i];
457          dataArray[i] = c;          dataArray[i] = c;
458      }      }
459    
460        // sample IDs now = mesh node/element IDs
461        sampleID = *requiredIDs;
462    
463      return true;      return true;
464  }  }
465    

Legend:
Removed from v.2834  
changed lines
  Added in v.2880

  ViewVC Help
Powered by ViewVC 1.1.26