From fb6bee2ac444226cf80c5cd70e2909a214b845c9 Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Tue, 1 Apr 2008 18:20:02 +0000
Subject: [PATCH] MED Gauss points
---
Post/PViewDataGModel.cpp | 34 +++++++++++++++++++++++---------
Post/PViewDataGModel.h | 8 ++++++++
Post/PViewDataGModelIO.cpp | 40 +++++++++++++++++++++++---------------
3 files changed, 57 insertions(+), 25 deletions(-)
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index edd7b4062b..de305b1625 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -1,4 +1,4 @@
-// $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();
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index 2187ae3fcb..d467a61739 100644
--- a/Post/PViewDataGModel.h
+++ b/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
diff --git a/Post/PViewDataGModelIO.cpp b/Post/PViewDataGModelIO.cpp
index 006f948c9c..969ea0cb58 100644
--- a/Post/PViewDataGModelIO.cpp
+++ b/Post/PViewDataGModelIO.cpp
@@ -1,4 +1,4 @@
-// $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)
--
GitLab