smoothing for ElementData fields

Post/PViewDataGModel.cpp

-// $Id: PViewDataGModel.cpp,v 1.39 2008-03-30 10:25:09 geuzaine Exp $
+// $Id: PViewDataGModel.cpp,v 1.40 2008-03-30 13:21:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -24,6 +24,7 @@
 
 #include "PViewDataGModel.h"
 #include "MElement.h"
+#include "Numeric.h"
 #include "Message.h"
 
 PViewDataGModel::PViewDataGModel() 
@@ -180,6 +181,48 @@ int PViewDataGModel::getNumEdges(int step, int ent, int ele)
   return _steps[step]->getEntity(ent)->getMeshElement(ele)->getNumEdges();
 }
 
+void PViewDataGModel::smooth()
+{
+  if(_type == NodeData) return;
+  std::vector<stepData<double>*> _steps2;
+  for(unsigned int step = 0; step < _steps.size(); step++){
+    GModel *m = _steps[step]->getModel();
+    int numComp = _steps[step]->getNumComponents();
+    _steps2.push_back(new stepData<double>(m, numComp, _steps[step]->getFileName(),
+					   _steps[step]->getFileIndex()));
+    std::map<int, int> nodeConnect;
+    for(int ent = 0; ent < getNumEntities(step); ent++){
+      for(int ele = 0; ele < getNumElements(step, ent); ele++){
+	MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
+	for(int nod = 0; nod < e->getNumVertices(); nod++){
+	  MVertex *v = e->getVertex(nod);
+	  if(nodeConnect.count(v->getNum()))
+	    nodeConnect[v->getNum()]++;
+	  else
+	    nodeConnect[v->getNum()] = 1;
+	  double *d = _steps2.back()->getData(v->getNum(), true), val;
+	  for(int j = 0; j < numComp; j++)
+	    if(getValue(step, e->getNum(), j, val)) d[j] += val;
+	}
+      }
+    }
+    for(int i = 0; i < _steps2.back()->getNumData(); i++){
+      double *d = _steps2.back()->getData(i);
+      if(d){
+	double f = nodeConnect[i];
+	if(f) for(int j = 0; j < numComp; j++) d[j] /= f;
+	double s = ComputeScalarRep(numComp, d);
+	_steps2[step]->setMin(std::min(_steps2[step]->getMin(), s));
+	_steps2[step]->setMax(std::max(_steps2[step]->getMax(), s));
+      }
+    }
+  }
+  for(unsigned int i = 0; i < _steps.size(); i++) delete _steps[i];
+  _steps = _steps2;
+  _type = NodeData;
+  finalize();
+}
+
 bool PViewDataGModel::skipEntity(int step, int ent)
 {
   return !_steps[step]->getEntity(ent)->getVisibility();

Post/PViewDataGModel.h

@@ -95,7 +95,10 @@ class stepData{
   {
     if(allocIfNeeded){
       if(index >= getNumData()) resizeData(index + 100); // optimize this
-      if(!(*_data)[index]) (*_data)[index] = new real[_numComp * mult];
+      if(!(*_data)[index]){
+	(*_data)[index] = new real[_numComp * mult];
+	for(int i = 0; i < _numComp * mult; i++) (*_data)[index][i] = 0.;
+      }
     }
     else{
       if(index >= getNumData()) return 0;
@@ -150,6 +153,7 @@ class PViewDataGModel : public PViewData {
   int getNumComponents(int step, int ent, int ele);
   void getValue(int step, int ent, int ele, int node, int comp, double &val);
   int getNumEdges(int step, int ent, int ele);
+  void smooth();
   bool skipEntity(int step, int ent);
   bool skipElement(int step, int ent, int ele);
   bool hasTimeStep(int step);

Post/PViewDataGModelIO.cpp

-// $Id: PViewDataGModelIO.cpp,v 1.23 2008-03-30 11:53:04 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.24 2008-03-30 13:21:04 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -320,8 +320,6 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	double *d = _steps[step]->getData(num, true, mult);
 	for(int j = 0; j < numComp * mult; j++)
 	  d[j] = val[numComp * i + j];
-	for(int j = numComp; j < _steps[step]->getNumComponents(); j++)
-	  d[j] = 0.;
 	double s = ComputeScalarRep(_steps[step]->getNumComponents(), d);
 	_steps[step]->setMin(std::min(_steps[step]->getMin(), s));
 	_steps[step]->setMax(std::max(_steps[step]->getMax(), s));