MED Gauss points

Post/PViewDataGModel.cpp

-// $Id: PViewDataGModel.cpp,v 1.44 2008-03-31 21:12:41 geuzaine Exp $
+// $Id: PViewDataGModel.cpp,v 1.45 2008-04-01 18:20:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -140,12 +140,13 @@ int PViewDataGModel::getDimension(int step, int ent, int ele)
 int PViewDataGModel::getNumNodes(int step, int ent, int ele)
 {
   // no sanity checks (assumed to be guarded by skipElement)
+  MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
   if(_type == GaussPointData){
-    return 1; // FIXME
+    return _steps[step]->getGaussPoints(e->getTypeForMSH()).size() / 3;
   }
   else{
-    return _steps[step]->getEntity(ent)->getMeshElement(ele)->getNumVertices();
-    //return _steps[step]->getEntity(ent)->getMeshElement(ele)->getNumPrimaryVertices();
+    //return e->getNumVertices();
+    return e->getNumPrimaryVertices();
   }
 }
 
@@ -153,14 +154,29 @@ void PViewDataGModel::getNode(int step, int ent, int ele, int nod,
                               double &x, double &y, double &z)
 {
   // no sanity checks (assumed to be guarded by skipElement)
+  MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
   if(_type == GaussPointData){ 
-    // FIXME
-    MElement *e = _steps[step]->getEntity(ent)->getMeshElement(ele);
-    SPoint3 bc = e->barycenter();
-    x = bc.x(); y = bc.y(); z = bc.z();
+    std::vector<double> &p(_steps[step]->getGaussPoints(e->getTypeForMSH()));
+    if(p[0] == 1.e22){
+      // hack: the points are the element vertices
+      x = e->getVertex(nod)->x();
+      y = e->getVertex(nod)->y();
+      z = e->getVertex(nod)->z();
+    }
+    else{
+      double vx[8], vy[8], vz[8];
+      for(int i = 0; i < e->getNumPrimaryVertices(); i++){
+	vx[i] = e->getVertex(i)->x();
+	vy[i] = e->getVertex(i)->y();
+	vz[i] = e->getVertex(i)->z();
+      }
+      x = e->interpolate(vx, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2]);
+      y = e->interpolate(vy, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2]);
+      z = e->interpolate(vz, p[3 * nod], p[3 * nod + 1], p[3 * nod + 2]);
+    }
   }
   else{
-    MVertex *v = _steps[step]->getEntity(ent)->getMeshElement(ele)->getVertex(nod);
+    MVertex *v = e->getVertex(nod);
     x = v->x();
     y = v->y();
     z = v->z();

Post/PViewDataGModel.h

@@ -50,6 +50,9 @@ class stepData{
   // the data and 2) not to store any additional info in MVertex or
   // MElement)
   std::vector<real*> *_data;
+  // a vector, indexed by MSH element type, of Gauss point locations
+  // in parametric space
+  std::vector<std::vector<double> > _gaussPoints;
  public:
   stepData(GModel *model, int numComp, std::string fileName="", int fileIndex=-1, 
 	   double time=0., double min=VAL_INF, double max=-VAL_INF)
@@ -114,6 +117,11 @@ class stepData{
       _data = 0;
     }
   }
+  std::vector<double> &getGaussPoints(int msh)
+  {
+    if(_gaussPoints.size() <= msh) _gaussPoints.resize(msh + 1);
+    return _gaussPoints[msh];
+  }
 };
 
 // data container using elements from a GModel

Post/PViewDataGModelIO.cpp

-// $Id: PViewDataGModelIO.cpp,v 1.33 2008-04-01 13:41:33 geuzaine Exp $
+// $Id: PViewDataGModelIO.cpp,v 1.34 2008-04-01 18:20:02 geuzaine Exp $
 //
 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
 //
@@ -275,22 +275,30 @@ bool PViewDataGModel::readMED(std::string fileName, int fileIndex)
 	return false;
       }
 
-      // read Gauss point data (if locname is MED_GAUSS_ELNO, the
-      // points are the element vertices)
-      if(_type == GaussPointData && std::string(locname) != MED_GAUSS_ELNO){
-	std::vector<med_float> refcoo((ele % 100) * (ele / 100));
-	std::vector<med_float> gscoo(ngauss * ele / 100);
-	std::vector<med_float> wg(ngauss);
-	if(MEDgaussLire(fid, &refcoo[0], &gscoo[0], &wg[0], MED_FULL_INTERLACE,
-			locname) < 0){
-	  Msg(GERROR, "Could not read Gauss points");
-	  return false;
+      // read Gauss point data
+      if(_type == GaussPointData){
+	std::vector<double> &p(_steps[step]->getGaussPoints(med2mshElementType(ele)));
+	if(std::string(locname) == MED_GAUSS_ELNO){
+	  // hack: the points are the vertices
+	  p.resize(ngauss * 3, 1.e22);
+	}
+	else{
+	  int dim = ele / 100;
+	  std::vector<med_float> refcoo((ele % 100) * dim);
+	  std::vector<med_float> gscoo(ngauss * dim);
+	  std::vector<med_float> wg(ngauss);
+	  if(MEDgaussLire(fid, &refcoo[0], &gscoo[0], &wg[0], MED_FULL_INTERLACE,
+			  locname) < 0){
+	    Msg(GERROR, "Could not read Gauss points");
+	    return false;
+	  }
+	  // we should check that refcoo corresponds to our internal
+	  // reference element
+	  for(unsigned int i = 0; i < gscoo.size(); i++){
+	    p.push_back(gscoo[i]);
+	    if(i % dim == dim - 1) for(int j = 0; j < 3 - dim; j++) p.push_back(0.); 
+	  }
 	}
-	// FIXME: store this in stepData, e.g. in a vector indexed by
-	// mshEleType std::vector<std::vector<u,v,w, u,v,w, u,v,w, ...> >
-	// (ele/100==geo dim, ele%100==num nodes)
-	// int msh = med2mshElementTupe(ele);
-	// gaussPointsCoordinates[msh] = gscoo; // zero pad
       }
 
       // compute profile (indices in full array of entities of given type)