diff --git a/Fltk/GUI.cpp b/Fltk/GUI.cpp
index 90e588a0921b5253150a75cd0f6f47fa45c6dd4c..02154eab05a8eda39828a6ebf02e754d45e19ab9 100644
--- a/Fltk/GUI.cpp
+++ b/Fltk/GUI.cpp
@@ -1,4 +1,4 @@
-// $Id: GUI.cpp,v 1.662 2008-03-18 08:41:21 remacle Exp $
+// $Id: GUI.cpp,v 1.663 2008-03-18 19:30:13 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -72,7 +72,6 @@ Fl_Menu_Item m_menubar_table[] = {
{"&Tools", 0, 0, 0, FL_SUBMENU},
{"&Options...", FL_CTRL+FL_SHIFT+'n', (Fl_Callback *)options_cb, 0},
{"Pl&ugins...", FL_CTRL+FL_SHIFT+'u', (Fl_Callback *)view_plugin_cb, (void*)(-1)},
- {"&Fields...", FL_CTRL+FL_SHIFT+'f', (Fl_Callback *)view_field_cb, (void*)(-1)},
{"&Visibility", FL_CTRL+FL_SHIFT+'v', (Fl_Callback *)visibility_cb, 0},
{"&Clipping Planes", FL_CTRL+FL_SHIFT+'c', (Fl_Callback *)clip_cb, 0},
{"&Manipulator", FL_CTRL+FL_SHIFT+'m', (Fl_Callback *)manip_cb, 0, FL_MENU_DIVIDER},
@@ -114,7 +113,6 @@ Fl_Menu_Item m_sys_menubar_table[] = {
{"Tools", 0, 0, 0, FL_SUBMENU},
{"Options...", FL_META+FL_SHIFT+'n', (Fl_Callback *)options_cb, 0},
{"Plugins...", FL_META+FL_SHIFT+'u', (Fl_Callback *)view_plugin_cb, (void*)(-1)},
- {"Fields...", FL_META+FL_SHIFT+'f', (Fl_Callback *)view_field_cb, (void*)(-1)},
{"Visibility", FL_META+FL_SHIFT+'v', (Fl_Callback *)visibility_cb, 0},
{"Clipping Planes", FL_META+FL_SHIFT+'c', (Fl_Callback *)clip_cb, 0},
{"Manipulator", FL_META+FL_SHIFT+'m', (Fl_Callback *)manip_cb, 0, FL_MENU_DIVIDER},
@@ -312,6 +310,7 @@ Context_Item menu_mesh[] = {
};
Context_Item menu_mesh_define[] = {
{"1Mesh>Define", NULL} ,
+ {"Fields", (Fl_Callback *)view_field_cb, (void*)(-1) },
{"Characteristic length", (Fl_Callback *)mesh_define_length_cb } ,
{"Recombine", (Fl_Callback *)mesh_define_recombine_cb } ,
{"Transfinite", (Fl_Callback *)mesh_define_transfinite_cb } ,
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 0ce655c6ab09d631abf3ea247ea3c0f79ccc6833..a3c97196422267afb7a69cbd8e8dbce8e36aed26 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1,4 +1,4 @@
-// $Id: GModel.cpp,v 1.71 2008-03-18 08:41:21 remacle Exp $
+// $Id: GModel.cpp,v 1.72 2008-03-18 19:30:13 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -99,8 +99,8 @@ void GModel::destroy()
MVertex::resetGlobalNumber();
MElement::resetGlobalNumber();
- fields.reset();
#if !defined(HAVE_GMSH_EMBEDDED)
+ fields.reset();
gmshSurface::reset();
BGMReset();
#endif
diff --git a/Geo/GModel.h b/Geo/GModel.h
index 022a4e36342b3fc9fb6ecdb7ea479dbb7ec5cb29..9958c592725913e7b0936823d75c9d2d336faf88 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -28,7 +28,9 @@
#include "GFace.h"
#include "GRegion.h"
#include "SBoundingBox3d.h"
+#if !defined(HAVE_GMSH_EMBEDDED)
#include "Field.h"
+#endif
class GEO_Internals;
class OCC_Internals;
@@ -38,7 +40,6 @@ class smooth_normals;
class GModel
{
private:
- std::set<Field*>mesh_size_fields;
// Vertex cache to speed-up direct access by vertex number (used for
// post-processing I/O)
std::vector<MVertex*> _vertexVectorCache;
@@ -248,8 +249,10 @@ class GModel
// Med interface ("Modele d'Echange de Donnees")
int writeMED(const std::string &name);
- // Characteristic Lengths fields
- FieldManager fields;
+#if !defined(HAVE_GMSH_EMBEDDED)
+ // Characteristic Lengths fields
+ FieldManager fields;
+#endif
};
#endif
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 5bedf5f37e173f16ded63f1548432edea5a79e5b..f6db22a9fbdae7f70d9a8c6babd05069be1a2503 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1,4 +1,4 @@
-// $Id: GModelIO_Mesh.cpp,v 1.42 2008-03-11 22:51:08 geuzaine Exp $
+// $Id: GModelIO_Mesh.cpp,v 1.43 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -261,7 +261,7 @@ int GModel::readMSH(const std::string &name)
double version = 1.0;
bool binary = false, swap = false;
- char str[256];
+ char str[256] = "XXX";
std::map<int, std::vector<MVertex*> > points;
std::map<int, std::vector<MElement*> > elements[7];
std::map<int, std::map<int, std::string> > physicals[4];
@@ -273,11 +273,11 @@ int GModel::readMSH(const std::string &name)
while(1) {
- do {
+ while(str[0] != '$'){
if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
- } while(str[0] != '$');
-
+ }
+
if(feof(fp))
break;
@@ -467,7 +467,6 @@ int GModel::readMSH(const std::string &name)
if(!fgets(str, sizeof(str), fp) || feof(fp))
Msg(GERROR, "Prematured end of mesh file");
} while(str[0] != '$');
-
}
// store the elements in their associated elementary entity. If the
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index ce7f8fecb952b8d5ca205a57ff97a4b05fdb47b7..248f98f1895baee7386309bc969a6dab7251a759 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -1,4 +1,4 @@
-// $Id: MElement.cpp,v 1.58 2008-02-22 07:19:07 geuzaine Exp $
+// $Id: MElement.cpp,v 1.59 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -132,19 +132,6 @@ double MTetrahedron::etaShapeMeasure()
#endif
}
-void MTetrahedron::getMat(double mat[3][3])
-{
- mat[0][0] = _v[1]->x() - _v[0]->x();
- mat[0][1] = _v[2]->x() - _v[0]->x();
- mat[0][2] = _v[3]->x() - _v[0]->x();
- mat[1][0] = _v[1]->y() - _v[0]->y();
- mat[1][1] = _v[2]->y() - _v[0]->y();
- mat[1][2] = _v[3]->y() - _v[0]->y();
- mat[2][0] = _v[1]->z() - _v[0]->z();
- mat[2][1] = _v[2]->z() - _v[0]->z();
- mat[2][2] = _v[3]->z() - _v[0]->z();
-}
-
double MTetrahedron::getVolume()
{
double mat[3][3];
@@ -169,31 +156,6 @@ bool MTetrahedron::invertmapping(double *p, double *uvw, double tol)
return false;
}
-void MTetrahedron::circumcenter(double X[4], double Y[4], double Z[4], double *res)
-{
- double mat[3][3], b[3], dum;
- b[0] = X[1] * X[1] - X[0] * X[0] +
- Y[1] * Y[1] - Y[0] * Y[0] + Z[1] * Z[1] - Z[0] * Z[0];
- b[1] = X[2] * X[2] - X[1] * X[1] +
- Y[2] * Y[2] - Y[1] * Y[1] + Z[2] * Z[2] - Z[1] * Z[1];
- b[2] = X[3] * X[3] - X[2] * X[2] +
- Y[3] * Y[3] - Y[2] * Y[2] + Z[3] * Z[3] - Z[2] * Z[2];
- for(int i = 0; i < 3; i++)
- b[i] *= 0.5;
- mat[0][0] = X[1] - X[0];
- mat[0][1] = Y[1] - Y[0];
- mat[0][2] = Z[1] - Z[0];
- mat[1][0] = X[2] - X[1];
- mat[1][1] = Y[2] - Y[1];
- mat[1][2] = Z[2] - Z[1];
- mat[2][0] = X[3] - X[2];
- mat[2][1] = Y[3] - Y[2];
- mat[2][2] = Z[3] - Z[2];
- if(!sys3x3(mat, b, res, &dum)) {
- res[0] = res[1] = res[2] = 10.0e10;
- }
-}
-
int MHexahedron::getVolumeSign()
{
double mat[3][3];
@@ -487,161 +449,6 @@ void MElement::writeBDF(FILE *fp, int format, int elementary)
}
}
-bool MTriangle::invertmappingXY(double *p, double *uv, double tol)
-{
- double mat[2][2];
- double b[2];
- getMat(mat);
- b[0] = p[0] - getVertex(0)->x();
- b[1] = p[1] - getVertex(0)->y();
- sys2x2(mat, b, uv);
-
- if(uv[0] >= -tol &&
- uv[1] >= -tol &&
- uv[0] <= 1. + tol &&
- uv[1] <= 1. + tol &&
- 1. - uv[0] - uv[1] > -tol) {
- return true;
- }
- return false;
-}
-
-bool MTriangle::invertmappingUV(GFace* gf, double *p, double *uv, double tol)
-{
- double mat[2][2];
- double b[2];
- double u0, v0, u1, v1, u2, v2;
-
- parametricCoordinates(getVertex(0), gf, u0, v0);
- parametricCoordinates(getVertex(1), gf, u1, v1);
- parametricCoordinates(getVertex(2), gf, u2, v2);
-
- mat[0][0] = u1 - u0;
- mat[0][1] = u2 - u0;
- mat[1][0] = v1 - v0;
- mat[1][1] = v2 - v0;
-
- b[0] = p[0] - u0;
- b[1] = p[1] - v0;
- sys2x2(mat, b, uv);
-
- if(uv[0] >= -tol &&
- uv[1] >= -tol &&
- uv[0] <= 1. + tol &&
- uv[1] <= 1. + tol &&
- 1. - uv[0] - uv[1] > -tol) {
- return true;
- }
- return false;
-}
-
-double MTriangle::getSurfaceUV(GFace *gf)
-{
- double u3, v3, u1, v1, u2, v2;
-
- parametricCoordinates(getVertex(0), gf, u1, v1);
- parametricCoordinates(getVertex(1), gf, u2, v2);
- parametricCoordinates(getVertex(2), gf, u3, v3);
-
- const double vv1 [2] = {u2 - u1, v2 - v1};
- const double vv2 [2] = {u3 - u1, v3 - v1};
-
- double s = vv1[0] * vv2[1] - vv1[1] * vv2[0];
- return s * 0.5;
-}
-
-double MTriangle::getSurfaceXY() const
-{
- const double x1 = _v[0]->x();
- const double x2 = _v[1]->x();
- const double x3 = _v[2]->x();
- const double y1 = _v[0]->y();
- const double y2 = _v[1]->y();
- const double y3 = _v[2]->y();
-
- const double v1 [2] = {x2 - x1, y2 - y1};
- const double v2 [2] = {x3 - x1, y3 - y1};
-
- double s = v1[0] * v2[1] - v1[1] * v2[0];
- return s * 0.5;
-}
-
-void MTriangle::circumcenterXYZ(double *p1, double *p2, double *p3,
- double *res, double *uv)
-{
- double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
- double v2[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
- double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
- double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
- double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
- norme(vx); norme(vy); norme(vz);
- double p1P[2] = {0.0, 0.0};
- double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
- double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
- double resP[2];
-
- circumcenterXY(p1P, p2P, p3P,resP);
-
- if(uv){
- double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
- {p2P[1] - p1P[1], p3P[1] - p1P[1]}};
- double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
- sys2x2(mat, rhs, uv);
- }
-
- res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
- res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
- res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
-}
-
-void MTriangle::circumcenterXY(double *p1, double *p2, double *p3, double *res)
-{
- double d, a1, a2, a3;
-
- const double x1 = p1[0];
- const double x2 = p2[0];
- const double x3 = p3[0];
- const double y1 = p1[1];
- const double y2 = p2[1];
- const double y3 = p3[1];
-
- d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
- if(d == 0.0) {
- Msg(WARNING, "Colinear points in circum circle computation");
- res[0] = res[1] = -99999.;
- return ;
- }
-
- a1 = x1 * x1 + y1 * y1;
- a2 = x2 * x2 + y2 * y2;
- a3 = x3 * x3 + y3 * y3;
- res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
- res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
-}
-
-void MTriangle::circumcenterUV(GFace *gf, double *res)
-{
- double u3, v3, u1, v1, u2, v2;
-
- parametricCoordinates(getVertex(0), gf, u1, v1);
- parametricCoordinates(getVertex(1), gf, u2, v2);
- parametricCoordinates(getVertex(2), gf, u3, v3);
-
- double p1[2] = {u1, v1};
- double p2[2] = {u2, v2};
- double p3[2] = {u3, v3};
-
- circumcenterXY(p1, p2, p3, res);
-}
-
-void MTriangle::circumcenterXY(double *res) const
-{
- double p1[2] = {_v[0]->x(), _v[0]->y()};
- double p2[2] = {_v[1]->x(), _v[1]->y()};
- double p3[2] = {_v[2]->x(), _v[2]->y()};
- circumcenterXY(p1, p2, p3, res);
-}
-
void MTriangle::jac(int ord, MVertex *vs[], double uu, double vv, double j[2][3])
{
#if defined(HAVE_GMSH_EMBEDDED)
@@ -732,7 +539,8 @@ void MTriangleN::jac(double uu, double vv , double j[2][3])
MTriangle::jac(_order, &(*(_vs.begin())), uu, vv, j);
}
-void MTriangleN::pnt(double uu, double vv, SPoint3 &p){
+void MTriangleN::pnt(double uu, double vv, SPoint3 &p)
+{
MTriangle::pnt(_order, &(*(_vs.begin())), uu, vv, p);
}
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 04a55d28bc5ce946960d4d71e3c2df8755430dfb..27de345dcf0c027b00902857d33d7a3505f2f354 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -318,23 +318,14 @@ class MTriangle : public MElement {
}
~MTriangle(){}
virtual int getDim(){ return 2; }
- virtual void getMat(double mat[2][2])
+ void getMat(double mat[2][2])
{
mat[0][0] = _v[1]->x() - _v[0]->x();
mat[0][1] = _v[2]->x() - _v[0]->x();
mat[1][0] = _v[1]->y() - _v[0]->y();
mat[1][1] = _v[2]->y() - _v[0]->y();
}
- void circumcenterXY(double *res) const;
virtual double gammaShapeMeasure();
- void circumcenterUV(GFace*,double *res);
- static void circumcenterXYZ(double *p1, double *p2, double *p3, double *res,
- double *uv = 0);
- static void circumcenterXY(double *p1, double *p2, double *p3, double *res);
- double getSurfaceXY() const;
- double getSurfaceUV(GFace*);
- bool invertmappingXY(double *p, double *uv, double tol = 1.e-8);
- bool invertmappingUV(GFace*,double *p, double *uv, double tol = 1.e-8);
virtual int getNumVertices(){ return 3; }
virtual MVertex *getVertex(int num){ return _v[num]; }
virtual MVertex *getOtherVertex(MVertex *v1, MVertex *v2)
@@ -753,21 +744,24 @@ class MTetrahedron : public MElement {
{
MVertex *tmp = _v[0]; _v[0] = _v[1]; _v[1] = tmp;
}
- virtual void getMat(double mat[3][3]);
+ void getMat(double mat[3][3])
+ {
+ mat[0][0] = _v[1]->x() - _v[0]->x();
+ mat[0][1] = _v[2]->x() - _v[0]->x();
+ mat[0][2] = _v[3]->x() - _v[0]->x();
+ mat[1][0] = _v[1]->y() - _v[0]->y();
+ mat[1][1] = _v[2]->y() - _v[0]->y();
+ mat[1][2] = _v[3]->y() - _v[0]->y();
+ mat[2][0] = _v[1]->z() - _v[0]->z();
+ mat[2][1] = _v[2]->z() - _v[0]->z();
+ mat[2][2] = _v[3]->z() - _v[0]->z();
+ }
virtual double getVolume();
virtual int getVolumeSign(){ return (getVolume() >= 0) ? 1 : -1; }
virtual double gammaShapeMeasure();
virtual double etaShapeMeasure();
// returns true if the point lies inside the tet
bool invertmapping(double *p, double *uvw, double tol = 1.e-8);
- static void circumcenter(double X[4],double Y[4],double Z[4],double *res);
- void circumcenter(double *res)
- {
- double X[4] = {_v[0]->x(), _v[1]->x(), _v[2]->x(), _v[3]->x()};
- double Y[4] = {_v[0]->y(), _v[1]->y(), _v[2]->y(), _v[3]->y()};
- double Z[4] = {_v[0]->z(), _v[1]->z(), _v[2]->z(), _v[3]->z()};
- MTetrahedron::circumcenter(X, Y, Z, res);
- }
};
class MTetrahedron10 : public MTetrahedron {
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 30bdabc60dacd32b7955e1e25e233eb5dd04e0e7..9edb7e115fcac014a5f1e30af9502ef7770ae446 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.125 2008-03-12 14:52:58 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.126 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -41,10 +41,10 @@
extern Context_T CTX;
-void fourthPoint (double *p1, double *p2, double *p3, double *p4)
+void fourthPoint(double *p1, double *p2, double *p3, double *p4)
{
double c[3];
- MTriangle::circumcenterXYZ(p1, p2, p3, c);
+ circumCenterXYZ(p1, p2, p3, c);
double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
double vz[3]; prodve(vx, vy, vz);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index ae7b68ebe57cc5e401215ea1d98edc4ae50c6863..a9297656cbfdb79bc241b7e3f3089ecc5a1528f6 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceDelaunayInsertion.cpp,v 1.13 2008-03-12 14:52:58 geuzaine Exp $
+// $Id: meshGFaceDelaunayInsertion.cpp,v 1.14 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -31,6 +31,58 @@
#include <map>
#include <algorithm>
+void circumCenterXY(double *p1, double *p2, double *p3, double *res)
+{
+ double d, a1, a2, a3;
+
+ const double x1 = p1[0];
+ const double x2 = p2[0];
+ const double x3 = p3[0];
+ const double y1 = p1[1];
+ const double y2 = p2[1];
+ const double y3 = p3[1];
+
+ d = 2. * (double)(y1 * (x2 - x3) + y2 * (x3 - x1) + y3 * (x1 - x2));
+ if(d == 0.0) {
+ Msg(WARNING, "Colinear points in circum circle computation");
+ res[0] = res[1] = -99999.;
+ return ;
+ }
+
+ a1 = x1 * x1 + y1 * y1;
+ a2 = x2 * x2 + y2 * y2;
+ a3 = x3 * x3 + y3 * y3;
+ res[0] = (double)((a1 * (y3 - y2) + a2 * (y1 - y3) + a3 * (y2 - y1)) / d);
+ res[1] = (double)((a1 * (x2 - x3) + a2 * (x3 - x1) + a3 * (x1 - x2)) / d);
+}
+
+void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv)
+{
+ double v1[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
+ double v2[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
+ double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
+ double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
+ double vz[3]; prodve(vx, vy, vz); prodve(vz, vx, vy);
+ norme(vx); norme(vy); norme(vz);
+ double p1P[2] = {0.0, 0.0};
+ double p2P[2]; prosca(v1, vx, &p2P[0]); prosca(v1, vy, &p2P[1]);
+ double p3P[2]; prosca(v2, vx, &p3P[0]); prosca(v2, vy, &p3P[1]);
+ double resP[2];
+
+ circumCenterXY(p1P, p2P, p3P,resP);
+
+ if(uv){
+ double mat[2][2] = {{p2P[0] - p1P[0], p3P[0] - p1P[0]},
+ {p2P[1] - p1P[1], p3P[1] - p1P[1]}};
+ double rhs[2] = {resP[0] - p1P[0], resP[1] - p1P[1]};
+ sys2x2(mat, rhs, uv);
+ }
+
+ res[0] = p1[0] + resP[0] * vx[0] + resP[1] * vy[0];
+ res[1] = p1[1] + resP[0] * vx[1] + resP[1] * vy[1];
+ res[2] = p1[2] + resP[0] * vx[2] + resP[1] * vy[2];
+}
+
bool circumCenterMetricInTriangle(MTriangle *base,
const double *metric,
const std::vector<double> &Us,
@@ -147,7 +199,7 @@ MTri3::MTri3(MTriangle *t, double lc) : deleted(false), base(t)
double pb[3] = {base->getVertex(1)->x(), base->getVertex(1)->y(), base->getVertex(1)->z()};
double pc[3] = {base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
double center[3];
- base->circumcenterXYZ(pa, pb, pc, center);
+ circumCenterXYZ(pa, pb, pc, center);
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
@@ -273,13 +325,13 @@ bool circUV(MTriangle *t, std::vector<double> & Us, std::vector<double> &Vs,
double u1 [3] = {Us[t->getVertex(0)->getNum()], Vs[t->getVertex(0)->getNum()], 0};
double u2 [3] = {Us[t->getVertex(1)->getNum()], Vs[t->getVertex(1)->getNum()], 0};
double u3 [3] = {Us[t->getVertex(2)->getNum()], Vs[t->getVertex(2)->getNum()], 0};
- t->circumcenterXY(u1, u2, u3, res);
+ circumCenterXY(u1, u2, u3, res);
return true;
double p1 [3] = {t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z()};
double p2 [3] = {t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z()};
double p3 [3] = {t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z()};
double resxy[3], uv[2];
- t->circumcenterXYZ(p1, p2, p3, resxy,uv);
+ circumCenterXYZ(p1, p2, p3, resxy,uv);
return true;
}
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 377e0a08e96f312f64f0ad3e4324a7c49ba49ad5..ab9160c99682f6b279ee4c11d9c6561329e17eb8 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -35,6 +35,7 @@ int inCircumCircleAniso(GFace *gf, double *p1, double *p2, double *p3, double *p
double *metric);
int inCircumCircleAniso(GFace *gf, MTriangle *base, const double *uv, const double *metric,
const std::vector<double> &Us, const std::vector<double> &Vs);
+void circumCenterXYZ(double *p1, double *p2, double *p3, double *res, double *uv=0);
void circumCenterMetric(MTriangle *base,
const double *metric,
const std::vector<double> &Us,
@@ -74,8 +75,6 @@ class MTri3
{
return inCircumCircle(v->x(), v->y());
}
- double getSurfaceXY () const { return base -> getSurfaceXY(); }
- double getSurfaceUV (GFace* gf) const { return base -> getSurfaceUV(gf); }
inline void setDeleted (bool d){ deleted = d; }
inline bool assertNeigh() const
{
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index a6f8af177a4c1f936b3adff5787abe632fdc7ba3..82df613ff9003101164d1df4a0d8925d4e53d89a 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.40 2008-03-12 08:36:49 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.41 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -809,7 +809,7 @@ void insertVerticesInRegion (GRegion *gr)
vSizes.size(), worst->getRadius());
if(worst->getRadius() < 1) break;
double center[3];
- worst->tet()->circumcenter(center);
+ worst->circumcenter(center);
double uvw[3];
bool inside = worst->tet()->invertmapping(center, uvw);
if(inside){
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index c4971ec0249ec8c69c5d27944586bf1114c701f4..ff09da3593f4c0c4eaa4ad31f8c99592fe518349 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -25,6 +25,7 @@
#include <map>
#include <stack>
#include "MElement.h"
+#include "Numeric.h"
#include "BackgroundMesh.h"
#include "qualityMeasures.h"
@@ -86,18 +87,47 @@ class MTet4
double vol;
circum_radius = qmTet(t, qm, &vol);
}
-
+ void circumcenter(double *res)
+ {
+ MVertex *v0 = base->getVertex(0);
+ MVertex *v1 = base->getVertex(1);
+ MVertex *v2 = base->getVertex(2);
+ MVertex *v3 = base->getVertex(3);
+ double X[4] = {v0->x(), v1->x(), v2->x(), v3->x()};
+ double Y[4] = {v0->y(), v1->y(), v2->y(), v3->y()};
+ double Z[4] = {v0->z(), v1->z(), v2->z(), v3->z()};
+ double b[3], mat[3][3], dum;
+ b[0] = X[1] * X[1] - X[0] * X[0] +
+ Y[1] * Y[1] - Y[0] * Y[0] + Z[1] * Z[1] - Z[0] * Z[0];
+ b[1] = X[2] * X[2] - X[1] * X[1] +
+ Y[2] * Y[2] - Y[1] * Y[1] + Z[2] * Z[2] - Z[1] * Z[1];
+ b[2] = X[3] * X[3] - X[2] * X[2] +
+ Y[3] * Y[3] - Y[2] * Y[2] + Z[3] * Z[3] - Z[2] * Z[2];
+ for(int i = 0; i < 3; i++)
+ b[i] *= 0.5;
+ mat[0][0] = X[1] - X[0];
+ mat[0][1] = Y[1] - Y[0];
+ mat[0][2] = Z[1] - Z[0];
+ mat[1][0] = X[2] - X[1];
+ mat[1][1] = Y[2] - Y[1];
+ mat[1][2] = Z[2] - Z[1];
+ mat[2][0] = X[3] - X[2];
+ mat[2][1] = Y[3] - Y[2];
+ mat[2][2] = Z[3] - Z[2];
+ if(!sys3x3(mat, b, res, &dum)) {
+ res[0] = res[1] = res[2] = 10.0e10;
+ }
+ }
void setup(MTetrahedron *t, std::vector<double> &sizes, std::vector<double> &sizesBGM)
{
base = t;
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double center[3];
- base->circumcenter(center);
+ circumcenter(center);
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
const double dz = base->getVertex(0)->z() - center[2];
- circum_radius = sqrt ( dx*dx + dy*dy + dz*dz);
-
+ circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
double lc1 = 0.25*(sizes[base->getVertex(0)->getNum()]+
sizes[base->getVertex(1)->getNum()]+
sizes[base->getVertex(2)->getNum()]+
@@ -107,17 +137,16 @@ class MTet4
sizesBGM[base->getVertex(2)->getNum()]+
sizesBGM[base->getVertex(3)->getNum()]);
double lc = Extend2dMeshIn3dVolumes() ? std::min(lc1, lcBGM) : lcBGM;
-
circum_radius /= lc;
deleted = false;
}
- inline GRegion *onWhat () const { return gr; }
- inline void setOnWhat (GRegion *g) { gr = g; }
- inline bool isDeleted () const { return deleted; }
- inline void forceRadius (double r){ circum_radius = r; }
- inline double getRadius () const { return circum_radius; }
- inline double getQuality () const { return circum_radius; }
- inline void setQuality (const double &q){ circum_radius = q; }
+ inline GRegion *onWhat() const { return gr; }
+ inline void setOnWhat(GRegion *g) { gr = g; }
+ inline bool isDeleted() const { return deleted; }
+ inline void forceRadius(double r){ circum_radius = r; }
+ inline double getRadius() const { return circum_radius; }
+ inline double getQuality() const { return circum_radius; }
+ inline void setQuality(const double &q){ circum_radius = q; }
inline MTetrahedron *tet() const { return base; }
inline MTetrahedron *&tet() { return base; }
inline void setNeigh(int iN, MTet4 *n) { neigh[iN] = n; }
@@ -203,14 +232,14 @@ class MTet4Factory
allSlots = new MTet4[s_alloc];
#endif
}
- ~MTet4Factory ()
+ ~MTet4Factory()
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
delete [] allSlots;
#endif
}
- MTet4 * Create (MTetrahedron * t, std::vector<double> &sizes,
- std::vector<double> &sizesBGM)
+ MTet4 *Create(MTetrahedron * t, std::vector<double> &sizes,
+ std::vector<double> &sizesBGM)
{
#ifdef _GMSH_PRE_ALLOCATE_STRATEGY_
MTet4 *t4 = getAnEmptySlot();
diff --git a/Plugin/Levelset.cpp b/Plugin/Levelset.cpp
index c4f9d726e5c6eb54d19a19cd6c2229acbf01bd6e..4640b95353d9e07fec74e47bf399dcbea4f0dba0 100644
--- a/Plugin/Levelset.cpp
+++ b/Plugin/Levelset.cpp
@@ -1,4 +1,4 @@
-// $Id: Levelset.cpp,v 1.39 2008-02-17 08:48:07 geuzaine Exp $
+// $Id: Levelset.cpp,v 1.40 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -24,40 +24,81 @@
#include "List.h"
#include "Numeric.h"
-GMSH_LevelsetPlugin::GMSH_LevelsetPlugin()
+static const int exn[13][12][2] = {
+ {{0,0}}, // point
+ {{0,1}}, // line
+ {{}}, // -
+ {{0,1}, {0,2}, {1,2}}, // triangle
+ {{0,1}, {0,3}, {1,2}, {2,3}}, // quad
+ {{}}, // -
+ {{0,1}, {0,2}, {0,3}, {1,2}, {1,3}, {2,3}}, // tetra
+ {{}}, // -
+ {{0,1}, {0,3}, {0,4}, {1,2}, {1,4}, {2,3}, {2,4}, {3,4}}, // pyramid
+ {{0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5}, {3,4}, {3,5}, {4,5}}, // prism
+ {{}}, {{}}, // -
+ {{0,1}, {0,3}, {0,4}, {1,2}, {1,5}, {2,3},
+ {2,6}, {3,7}, {4,5}, {4,7}, {5,6}, {6,7}} // hexa
+};
+
+static int numSimplexDec(int numEdges)
{
- _invert = 0.;
- _ref[0] = _ref[1] = _ref[2] = 0.;
- _valueIndependent = 0; // "moving" levelset
- _valueView = -1; // use same view for levelset and field data
- _valueTimeStep = -1; // use same time step in levelset and field data views
- _recurLevel = 4;
- _targetError = 0.;
- _extractVolume = 0; // to create isovolumes (keep all elements < or > levelset)
- _orientation = GMSH_LevelsetPlugin::NONE;
+ switch(numEdges){
+ case 4: return 2;
+ case 12: return 6;
+ case 9: return 3;
+ case 8: return 2;
+ default: return 1;
+ }
}
-static void affect(double *xpi, double *ypi, double *zpi, double valpi[12][9], int epi[12],
- int i,
- double *xp, double *yp, double *zp, double valp[12][9], int ep[12],
- int j, int nb)
+static void getSimplexDec(int numNodes, int numEdges, int i,
+ int &n0, int &n1, int &n2, int &n3,
+ int &nn, int &ne)
+{
+ static const int qua[2][3] = {{0,1,2}, {0,2,3}};
+ static const int hex[6][4] = {{0,1,3,7}, {0,4,1,7}, {1,4,5,7},
+ {1,2,3,7}, {1,6,2,7}, {1,5,6,7}};
+ static const int pri[3][4] = {{0,1,2,4}, {0,2,4,5}, {0,3,4,5}};
+ static const int pyr[2][4] = {{0,1,3,4}, {1,2,3,4}};
+ switch(numEdges){
+ case 4:
+ n0 = qua[i][0]; n1 = qua[i][1]; n2 = qua[i][2]; nn = 3; ne = 3;
+ break;
+ case 12:
+ n0 = hex[i][0]; n1 = hex[i][1]; n2 = hex[i][2]; n3 = hex[i][3]; nn = 4; ne = 6;
+ break;
+ case 9:
+ n0 = pri[i][0]; n1 = pri[i][1]; n2 = pri[i][2]; n3 = pri[i][3]; nn = 4; ne = 6;
+ break;
+ case 8:
+ n0 = pri[i][0]; n1 = pri[i][1]; n2 = pri[i][2]; n3 = pri[i][3]; nn = 4; ne = 6;
+ break;
+ default:
+ n0 = 0; n1 = 1; n2 = 2; n3 = 3; nn = numNodes; ne = numEdges;
+ break;
+ }
+}
+
+static void affect(double *xpi, double *ypi, double *zpi,
+ double valpi[12][9], int epi[12], int i,
+ double *xp, double *yp, double *zp,
+ double valp[12][9], int ep[12], int j, int nb)
{
xpi[i] = xp[j];
ypi[i] = yp[j];
zpi[i] = zp[j];
- for(int k = 0; k < nb; k++)
- valpi[i][k] = valp[j][k];
+ for(int k = 0; k < nb; k++) valpi[i][k] = valp[j][k];
epi[i] = ep[j];
}
-static void removeIdenticalNodes(int *np, int nbComp,
- double xp[12], double yp[12], double zp[12],
+static void removeIdenticalNodes(int *np, int numComp,
+ double xp[12], double yp[12], double zp[12],
double valp[12][9], int ep[12])
{
double xpi[12], ypi[12], zpi[12], valpi[12][9];
int epi[12];
- affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, numComp);
int npi = 1;
for(int j = 1; j < *np; j++) {
for(int i = 0; i < npi; i++) {
@@ -67,32 +108,29 @@ static void removeIdenticalNodes(int *np, int nbComp,
break;
}
if(i == npi-1) {
- affect(xpi, ypi, zpi, valpi, epi, npi, xp, yp, zp, valp, ep, j, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, npi, xp, yp, zp, valp, ep, j, numComp);
npi++;
break;
}
}
}
for(int i = 0; i < npi; i++)
- affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, nbComp);
+ affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, numComp);
*np = npi;
}
-static void reorderQuad(int nbComp,
- double xp[12], double yp[12], double zp[12],
+static void reorderQuad(int numComp, double xp[12], double yp[12], double zp[12],
double valp[12][9], int ep[12])
{
double xpi[1], ypi[1], zpi[1], valpi[1][9];
int epi[12];
- affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xp, yp, zp, valp, ep, 3, xp, yp, zp, valp, ep, 2, nbComp);
- affect(xp, yp, zp, valp, ep, 2, xpi, ypi, zpi, valpi, epi, 0, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xp, yp, zp, valp, ep, 3, xp, yp, zp, valp, ep, 2, numComp);
+ affect(xp, yp, zp, valp, ep, 2, xpi, ypi, zpi, valpi, epi, 0, numComp);
}
-static void reorderPrism(int nbComp,
- double xp[12], double yp[12], double zp[12],
- double valp[12][9], int ep[12],
- int nbCut, int exn[12][2])
+static void reorderPrism(int numComp, double xp[12], double yp[12], double zp[12],
+ double valp[12][9], int ep[12], int nbCut)
{
double xpi[6], ypi[6], zpi[6], valpi[6][9];
int epi[12];
@@ -100,135 +138,120 @@ static void reorderPrism(int nbComp,
if(nbCut == 3){
// 3 first nodes come from zero levelset intersection, next 3 are
// endpoints of relative edges
- affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 5, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 5, numComp);
for(int i = 0; i < 3; i++){
int edgecut = ep[i]-1;
for(int j = 0; j < 3; j++){
int p = -epi[j]-1;
- if(exn[edgecut][0] == p || exn[edgecut][1] == p)
- affect(xp, yp, zp, valp, ep, 3+i, xpi, ypi, zpi, valpi, epi, j, nbComp);
+ if(exn[9][edgecut][0] == p || exn[9][edgecut][1] == p)
+ affect(xp, yp, zp, valp, ep, 3+i, xpi, ypi, zpi, valpi, epi, j, numComp);
}
}
}
else if(nbCut == 4){
// 4 first nodes come from zero levelset intersection
- affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 0, xp, yp, zp, valp, ep, 0, numComp);
int edgecut = ep[0]-1;
int p0 = -ep[4]-1;
- if(exn[edgecut][0] == p0 || exn[edgecut][1] == p0){
- affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, nbComp);
- if(exn[ep[1]-1][0] == p0 || exn[ep[1]-1][1] == p0){
- affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+ if(exn[9][edgecut][0] == p0 || exn[9][edgecut][1] == p0){
+ affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 4, numComp);
+ if(exn[9][ep[1]-1][0] == p0 || exn[9][ep[1]-1][1] == p0){
+ affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
}
else{
- affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 5, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
}
}
else{
- affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 5, nbComp);
- if(exn[ep[1]-1][0] == p0 || exn[ep[1]-1][1] == p0){
- affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 1, xp, yp, zp, valp, ep, 5, numComp);
+ if(exn[9][ep[1]-1][0] == p0 || exn[9][ep[1]-1][1] == p0){
+ affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 1, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
}
else{
- affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, nbComp);
- affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, nbComp);
+ affect(xpi, ypi, zpi, valpi, epi, 2, xp, yp, zp, valp, ep, 3, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 3, xp, yp, zp, valp, ep, 1, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 4, xp, yp, zp, valp, ep, 4, numComp);
+ affect(xpi, ypi, zpi, valpi, epi, 5, xp, yp, zp, valp, ep, 2, numComp);
}
}
for(int i = 0; i < 6; i++)
- affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, nbComp);
+ affect(xp, yp, zp, valp, ep, i, xpi, ypi, zpi, valpi, epi, i, numComp);
}
}
-void GMSH_LevelsetPlugin::evalLevelset(int nbNod, int nbComp,
- double *x, double *y, double *z, double *val,
- double *levels, double *scalarVal)
+GMSH_LevelsetPlugin::GMSH_LevelsetPlugin()
{
- if(_valueIndependent) {
- for(int k = 0; k < nbNod; k++)
- levels[k] = levelset(x[k], y[k], z[k], 0.0);
- }
- else{
- for(int k = 0; k < nbNod; k++) {
- double *vals = &val[nbComp * k];
- switch(nbComp) {
- case 1: // scalar
- scalarVal[k] = vals[0];
- break;
- case 3 : // vector
- scalarVal[k] = sqrt(DSQR(vals[0]) + DSQR(vals[1]) + DSQR(vals[2]));
- break;
- case 9 : // tensor
- scalarVal[k] = ComputeVonMises(vals);
- break;
- }
- levels[k] = levelset(x[k], y[k], z[k], scalarVal[k]);
- }
- }
+ _invert = 0.;
+ _ref[0] = _ref[1] = _ref[2] = 0.;
+ _valueIndependent = 0; // "moving" levelset
+ _valueView = -1; // use same view for levelset and field data
+ _valueTimeStep = -1; // use same time step in levelset and field data views
+ _recurLevel = 4;
+ _targetError = 0.;
+ _extractVolume = 0; // to create isovolumes (keep all elements < or > levelset)
+ _orientation = GMSH_LevelsetPlugin::NONE;
}
-void GMSH_LevelsetPlugin::addElement(int timeStep, int np, int nbEdg, int dNbComp,
- double xp[12], double yp[12], double zp[12],
- double valp[12][9], std::vector<PViewDataList *> &out)
+void GMSH_LevelsetPlugin::_addElement(int step, int np, int numEdges, int numComp,
+ double xp[12], double yp[12], double zp[12],
+ double valp[12][9], std::vector<PViewDataList*> &out)
{
- // select the output data
- PViewDataList *data = _valueIndependent ? out[0] : out[timeStep];
+ PViewDataList *data = _valueIndependent ? out[0] : out[step];
List_T *list;
int *nbPtr;
switch(np){
case 1:
- if(dNbComp == 1) { list = data->SP; nbPtr = &data->NbSP; }
- else if(dNbComp == 3) { list = data->VP; nbPtr = &data->NbVP; }
+ if(numComp == 1) { list = data->SP; nbPtr = &data->NbSP; }
+ else if(numComp == 3) { list = data->VP; nbPtr = &data->NbVP; }
else { list = data->TP; nbPtr = &data->NbTP; }
break;
case 2:
- if(dNbComp == 1) { list = data->SL; nbPtr = &data->NbSL; }
- else if(dNbComp == 3) { list = data->VL; nbPtr = &data->NbVL; }
+ if(numComp == 1) { list = data->SL; nbPtr = &data->NbSL; }
+ else if(numComp == 3) { list = data->VL; nbPtr = &data->NbVL; }
else { list = data->TL; nbPtr = &data->NbTL; }
break;
case 3:
- if(dNbComp == 1) { list = data->ST; nbPtr = &data->NbST; }
- else if(dNbComp == 3) { list = data->VT; nbPtr = &data->NbVT; }
+ if(numComp == 1) { list = data->ST; nbPtr = &data->NbST; }
+ else if(numComp == 3) { list = data->VT; nbPtr = &data->NbVT; }
else { list = data->TT; nbPtr = &data->NbTT; }
break;
case 4:
- if(!_extractVolume || nbEdg <= 4){
- if(dNbComp == 1) { list = data->SQ; nbPtr = &data->NbSQ; }
- else if(dNbComp == 3) { list = data->VQ; nbPtr = &data->NbVQ; }
+ if(!_extractVolume || numEdges <= 4){
+ if(numComp == 1) { list = data->SQ; nbPtr = &data->NbSQ; }
+ else if(numComp == 3) { list = data->VQ; nbPtr = &data->NbVQ; }
else { list = data->TQ; nbPtr = &data->NbTQ; }
}
else{
- if(dNbComp == 1) { list = data->SS; nbPtr = &data->NbSS; }
- else if(dNbComp == 3) { list = data->VS; nbPtr = &data->NbVS; }
+ if(numComp == 1) { list = data->SS; nbPtr = &data->NbSS; }
+ else if(numComp == 3) { list = data->VS; nbPtr = &data->NbVS; }
else { list = data->TS; nbPtr = &data->NbTS; }
}
break;
case 5:
- if(dNbComp == 1) { list = data->SY; nbPtr = &data->NbSY; }
- else if(dNbComp == 3) { list = data->VY; nbPtr = &data->NbVY; }
+ if(numComp == 1) { list = data->SY; nbPtr = &data->NbSY; }
+ else if(numComp == 3) { list = data->VY; nbPtr = &data->NbVY; }
else { list = data->TY; nbPtr = &data->NbTY; }
break;
case 6:
- if(dNbComp == 1) { list = data->SI; nbPtr = &data->NbSI; }
- else if(dNbComp == 3) { list = data->VI; nbPtr = &data->NbVI; }
+ if(numComp == 1) { list = data->SI; nbPtr = &data->NbSI; }
+ else if(numComp == 3) { list = data->VI; nbPtr = &data->NbVI; }
else { list = data->TI; nbPtr = &data->NbTI; }
break;
- case 8: // should never happen
- if(dNbComp == 1) { list = data->SH; nbPtr = &data->NbSH; }
- else if(dNbComp == 3) { list = data->VH; nbPtr = &data->NbVH; }
+ case 8:
+ if(numComp == 1) { list = data->SH; nbPtr = &data->NbSH; }
+ else if(numComp == 3) { list = data->VH; nbPtr = &data->NbVH; }
else { list = data->TH; nbPtr = &data->NbTH; }
break;
default:
@@ -236,7 +259,7 @@ void GMSH_LevelsetPlugin::addElement(int timeStep, int np, int nbEdg, int dNbCom
}
// copy the elements in the output data
- if(!timeStep || !_valueIndependent) {
+ if(!step || !_valueIndependent) {
for(int k = 0; k < np; k++)
List_Add(list, &xp[k]);
for(int k = 0; k < np; k++)
@@ -246,525 +269,217 @@ void GMSH_LevelsetPlugin::addElement(int timeStep, int np, int nbEdg, int dNbCom
(*nbPtr)++;
}
for(int k = 0; k < np; k++)
- for(int l = 0; l < dNbComp; l++)
+ for(int l = 0; l < numComp; l++)
List_Add(list, &valp[k][l]);
}
-void GMSH_LevelsetPlugin::nonZeroLevelset(int timeStep,
- int nbNod, int nbEdg, int exn[12][2],
- double *x, double *y, double *z,
- double *iVal, int iNbComp,
- double *dVal, int dNbComp,
- std::vector<PViewDataList*> &out)
+void GMSH_LevelsetPlugin::_cutAndAddElements(PViewData *vdata, PViewData *wdata,
+ int ent, int ele, int step, int wstep,
+ double x[8], double y[8], double z[8],
+ double levels[8], double scalarValues[8],
+ std::vector<PViewDataList*> &out)
{
- double levels[8], scalarVal[8];
-
- evalLevelset(nbNod, iNbComp, x, y, z, iVal, levels, scalarVal);
-
- int add = 1;
- for(int k = 0; k < nbNod; k++){
- if((_extractVolume < 0. && levels[k] > 0.) ||
- (_extractVolume > 0. && levels[k] < 0.)){
- add = 0;
- break;
- }
- }
-
- if(add){
- double xp[12], yp[12], zp[12], valp[12][9];
- for(int k = 0; k < nbNod; k++){
- xp[k] = x[k];
- yp[k] = y[k];
- zp[k] = z[k];
- for(int l = 0; l < dNbComp; l++)
- valp[k][l] = dVal[dNbComp * k + l];
- }
- addElement(timeStep, nbNod, nbEdg, dNbComp, xp, yp, zp, valp, out);
- }
-}
+ int numNodes = vdata->getNumNodes(ent, ele);
+ int numEdges = vdata->getNumEdges(ent, ele);
+ int numComp = wdata->getNumComponents(ent, ele, wstep);
-int GMSH_LevelsetPlugin::zeroLevelset(int timeStep,
- int nbNod, int nbEdg, int exn[12][2],
- double *x, double *y, double *z,
- double *iVal, int iNbComp,
- double *dVal, int dNbComp,
- std::vector<PViewDataList*> &out)
-{
- double levels[8], scalarVal[8];
-
- evalLevelset(nbNod, iNbComp, x, y, z, iVal, levels, scalarVal);
-
- // interpolate the zero levelset and the field to plot on it
- double xp[12], yp[12], zp[12], valp[12][9];
- int np = 0, ep[12];
- for(int k = 0; k < nbEdg; k++) {
- if(levels[exn[k][0]] * levels[exn[k][1]] <= 0.0) {
- if(iVal && dVal) {
- double coef = InterpolateIso(x, y, z, levels, 0.0, exn[k][0], exn[k][1],
- &xp[np], &yp[np], &zp[np]);
- double *val1 = &dVal[dNbComp * exn[k][0]];
- double *val2 = &dVal[dNbComp * exn[k][1]];
- for(int l = 0; l < dNbComp; l++)
- valp[np][l] = coef * (val2[l] - val1[l]) + val1[l];
+ for(int simplex = 0; simplex < numSimplexDec(numEdges); simplex++){
+ double xp[12], yp[12], zp[12], valp[12][9];
+ int n[4], np = 0, ep[12], nsn, nse;
+ getSimplexDec(numNodes, numEdges, simplex, n[0], n[1], n[2], n[3], nsn, nse);
+
+ for(int i = 0; i < nse; i++){
+ int n0 = exn[nse][i][0], n1 = exn[nse][i][1];
+ if(levels[n[n0]] * levels[n[n1]] <= 0.) {
+ double c = InterpolateIso(x, y, z, levels, 0., n[n0], n[n1],
+ &xp[np], &yp[np], &zp[np]);
+ for(int comp = 0; comp < numComp; comp++){
+ double v0, v1;
+ wdata->getValue(ent, ele, n[n0], comp, wstep, v0);
+ wdata->getValue(ent, ele, n[n1], comp, wstep, v1);
+ valp[np][comp] = v0 + c * (v1 - v0);
+ }
+ ep[np++] = i + 1;
}
- ep[np] = k+1;
- np++;
}
- }
- if(!iVal || !dVal)
- return np;
-
- // Remove identical nodes (this can happen if an edge actually
- // belongs to the zero levelset, i.e., if levels[] * levels[] == 0)
- if(np > 1)
- removeIdenticalNodes(&np, dNbComp, xp, yp, zp, valp, ep);
-
- if(nbEdg > 4 && np < 3) // 3D input should only lead to 2D output
- return 0;
- else if(nbEdg > 1 && np < 2) // 2D input should only lead to 1D output
- return 0;
- else if(np < 1 || np > 4) // can't deal with this
- return 0;
-
- // avoid "butterflies"
- if(np == 4)
- reorderQuad(dNbComp, xp, yp, zp, valp, ep);
-
- // orient the triangles and the quads to get the normals right
- if(!_extractVolume && (np == 3 || np == 4)) {
- if(!timeStep || !_valueIndependent) {
- // test this only once for spatially-fixed views
- double v1[3] = { xp[2] - xp[0], yp[2] - yp[0], zp[2] - zp[0] };
- double v2[3] = { xp[1] - xp[0], yp[1] - yp[0], zp[1] - zp[0] };
- double gr[3], n[3];
- prodve(v1, v2, n);
- switch (_orientation) {
- case MAP:
- gradSimplex(x, y, z, scalarVal, gr);
- prosca(gr, n, &_invert);
- break;
- case PLANE:
- prosca(n, _ref, &_invert);
- break;
- case SPHERE:
- gr[0] = xp[0] - _ref[0];
- gr[1] = yp[0] - _ref[1];
- gr[2] = zp[0] - _ref[2];
- prosca(gr, n, &_invert);
- case NONE:
- default:
- break;
+ // Remove identical nodes (this can happen if an edge actually
+ // belongs to the zero levelset, i.e., if levels[] * levels[] == 0)
+ if(np > 1) removeIdenticalNodes(&np, numComp, xp, yp, zp, valp, ep);
+ if(numEdges > 4 && np < 3) // 3D input should only lead to 2D output
+ continue;
+ else if(numEdges > 1 && np < 2) // 2D input should only lead to 1D output
+ continue;
+ else if(np < 1 || np > 4) // can't deal with this
+ continue;
+
+ // avoid "butterflies"
+ if(np == 4) reorderQuad(numComp, xp, yp, zp, valp, ep);
+
+ // orient the triangles and the quads to get the normals right
+ if(!_extractVolume && (np == 3 || np == 4)) {
+ if(!step || !_valueIndependent) {
+ // test this only once for spatially-fixed views
+ double v1[3] = {xp[2] - xp[0], yp[2] - yp[0], zp[2] - zp[0]};
+ double v2[3] = {xp[1] - xp[0], yp[1] - yp[0], zp[1] - zp[0]};
+ double gr[3], normal[3];
+ prodve(v1, v2, normal);
+ switch (_orientation) {
+ case MAP:
+ gradSimplex(x, y, z, scalarValues, gr);
+ prosca(gr, normal, &_invert);
+ break;
+ case PLANE:
+ prosca(normal, _ref, &_invert);
+ break;
+ case SPHERE:
+ gr[0] = xp[0] - _ref[0];
+ gr[1] = yp[0] - _ref[1];
+ gr[2] = zp[0] - _ref[2];
+ prosca(gr, normal, &_invert);
+ case NONE:
+ default:
+ break;
+ }
+ }
+ if(_invert > 0.) {
+ double xpi[12], ypi[12], zpi[12], valpi[12][9];
+ int epi[12];
+ for(int k = 0; k < np; k++)
+ affect(xpi, ypi, zpi, valpi, epi, k, xp, yp, zp, valp, ep, k, numComp);
+ for(int k = 0; k < np; k++)
+ affect(xp, yp, zp, valp, ep, k, xpi, ypi, zpi, valpi, epi, np - k - 1, numComp);
}
}
- if(_invert > 0.) {
- double xpi[12], ypi[12], zpi[12], valpi[12][9];
- int epi[12];
- for(int k = 0; k < np; k++)
- affect(xpi, ypi, zpi, valpi, epi, k, xp, yp, zp, valp, ep, k, dNbComp);
- for(int k = 0; k < np; k++)
- affect(xp, yp, zp, valp, ep, k, xpi, ypi, zpi, valpi, epi, np-k-1, dNbComp);
- }
- }
- // if we compute isovolumes, add the nodes on the chosen side
- if(_extractVolume){
- // FIXME: when cutting 2D views, the elts we get can have the wrong
- // orientation
- int nbCut = np;
- for(int k = 0; k < nbNod; k++){
- if((_extractVolume < 0. && levels[k] < 0.0) ||
- (_extractVolume > 0. && levels[k] > 0.0)){
- xp[np] = x[k];
- yp[np] = y[k];
- zp[np] = z[k];
- for(int l = 0; l < dNbComp; l++)
- valp[np][l] = dVal[dNbComp * k + l];
- ep[np] = -(k+1); // node num!
- np++;
+ // if we compute isovolumes, add the nodes on the chosen side
+ // (FIXME: when cutting 2D views, the elts can have the wrong
+ // orientation)
+ if(_extractVolume){
+ int nbCut = np;
+ for(int nod = 0; nod < nsn; nod++){
+ if((_extractVolume < 0. && levels[n[nod]] < 0.) ||
+ (_extractVolume > 0. && levels[n[nod]] > 0.)){
+ xp[np] = x[n[nod]];
+ yp[np] = y[n[nod]];
+ zp[np] = z[n[nod]];
+ for(int comp = 0; comp < numComp; comp++)
+ wdata->getValue(ent, ele, n[nod], comp, wstep, valp[np][comp]);
+ ep[np] = -(nod + 1); // store node num!
+ np++;
+ }
}
+ removeIdenticalNodes(&np, numComp, xp, yp, zp, valp, ep);
+ if(np == 4 && numEdges <= 4)
+ reorderQuad(numComp, xp, yp, zp, valp, ep);
+ if(np == 6)
+ reorderPrism(numComp, xp, yp, zp, valp, ep, nbCut);
+ if(np > 8) // can't deal with this
+ continue;
}
- removeIdenticalNodes(&np, dNbComp, xp, yp, zp, valp, ep);
- if(np == 4 && nbEdg <= 4)
- reorderQuad(dNbComp, xp, yp, zp, valp, ep);
- if(np == 6)
- reorderPrism(dNbComp, xp, yp, zp, valp, ep, nbCut, exn);
- if(np > 8) // can't deal with this
- return 0;
- }
- addElement(timeStep, np, nbEdg, dNbComp, xp, yp, zp, valp, out);
-
- return 0;
-}
-
-void GMSH_LevelsetPlugin::executeList(PViewDataList *iData, List_T *iList,
- int iNbElm, int iNbComp,
- PViewDataList *dData, List_T *dList,
- int dNbElm, int dNbComp,
- int nbNod, int nbEdg, int exn[12][2],
- std::vector<PViewDataList*> &out)
-{
- if(!iNbElm || !dNbElm)
- return;
-
- if(iNbElm != dNbElm) {
- Msg(GERROR, "Views have a different number of elements (%d != %d)", iNbElm, dNbElm);
- return;
+ _addElement(step, np, numEdges, numComp, xp, yp, zp, valp, out);
}
- int dTimeStep = _valueTimeStep;
- if(dTimeStep >= dData->getNumTimeSteps()) {
- Msg(GERROR, "Wrong time step %d in view", dTimeStep);
- return;
- }
-
- int iNb = List_Nbr(iList) / iNbElm;
- int dNb = List_Nbr(dList) / dNbElm;
- for(int i = 0, j = 0; i < List_Nbr(iList); i += iNb, j += dNb) {
- double *x = (double *)List_Pointer_Fast(iList, i);
- double *y = (double *)List_Pointer_Fast(iList, i + nbNod);
- double *z = (double *)List_Pointer_Fast(iList, i + 2 * nbNod);
-
- if(nbNod == 1 || nbNod == 2 || nbNod == 3 || (nbNod == 4 && nbEdg == 6)) {
- // easy for simplices: at most one element is created per time step
- for(int iTS = 0; iTS < iData->getNumTimeSteps(); iTS++) {
- int dTS = (dTimeStep < 0) ? iTS : dTimeStep;
- // don't compute the zero levelset of the value view
- if(dTimeStep < 0 || iData != dData || dTS != iTS) {
- double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod +
- iNbComp * nbNod * iTS);
- double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
- dNbComp * nbNod * dTS);
- zeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z,
- iVal, iNbComp, dVal, dNbComp, out);
- if(_extractVolume)
- nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z,
- iVal, iNbComp, dVal, dNbComp, out);
- }
+ if(_extractVolume){
+ // if we compute isovolumes, add full elements that are completely
+ // on one side
+ bool add = true;
+ for(int nod = 0; nod < numNodes; nod++){
+ if((_extractVolume < 0. && levels[nod] > 0.) ||
+ (_extractVolume > 0. && levels[nod] < 0.)){
+ add = false;
+ break;
}
}
- else{
- // we decompose the element into simplices
- MakeSimplex iDec(nbNod, iNbComp);
- MakeSimplex dDec(nbNod, dNbComp);
-
- int nbNodNew = iDec.numSimplexNodes();
- int nbEdgNew = (nbNodNew == 4) ? 6 : 3;
- int exnTri[12][2] = {{0,1},{0,2},{1,2}};
- int exnTet[12][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
- double xNew[4], yNew[4], zNew[4];
- double *iValNew = new double[4 * iNbComp];
- double *dValNew = new double[4 * dNbComp];
-
- if(_valueIndependent) {
- // since the intersection is value-independent, we can loop on
- // the time steps so that only one element gets created per
- // time step. This allows us to still easily generate a
- // *single* multi-timestep view.
- if(zeroLevelset(0, nbNod, nbEdg, exn, x, y, z, NULL, 0,
- NULL, 0, out)) {
- for(int k = 0; k < iDec.numSimplices(); k++) {
- for(int iTS = 0; iTS < iData->getNumTimeSteps(); iTS++) {
- int dTS = (dTimeStep < 0) ? iTS : dTimeStep;
- // don't compute the zero levelset of the value view
- if(dTimeStep < 0 || iData != dData || dTS != iTS) {
- double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod +
- iNbComp * nbNod * iTS);
- double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
- dNbComp * nbNod * dTS);
- iDec.decompose(k, x, y, z, iVal, xNew, yNew, zNew, iValNew);
- dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew);
- zeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri,
- xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
- if(_extractVolume)
- nonZeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri,
- xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
- }
- }
- }
- }
- else if(_extractVolume){
- for(int iTS = 0; iTS < iData->getNumTimeSteps(); iTS++) {
- int dTS = (dTimeStep < 0) ? iTS : dTimeStep;
- // don't compute the zero levelset of the value view
- if(dTimeStep < 0 || iData != dData || dTS != iTS) {
- double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod +
- iNbComp * nbNod * iTS);
- double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
- dNbComp * nbNod * dTS);
- nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z,
- iVal, iNbComp, dVal, dNbComp, out);
- }
- }
- }
+ if(add){
+ double xp[12], yp[12], zp[12], valp[12][9];
+ for(int nod = 0; nod < numNodes; nod++){
+ xp[nod] = x[nod];
+ yp[nod] = y[nod];
+ zp[nod] = z[nod];
+ for(int comp = 0; comp < numComp; comp++)
+ wdata->getValue(ent, ele, nod, comp, wstep, valp[nod][comp]);
}
- else{
- // since we generate one view for each time step, we can
- // generate multiple elements per time step without problem.
- for(int iTS = 0; iTS < iData->getNumTimeSteps(); iTS++) {
- int dTS = (dTimeStep < 0) ? iTS : dTimeStep;
- // don't compute the zero levelset of the value view
- if(dTimeStep < 0 || iData != dData || dTS != iTS) {
- double *iVal = (double *)List_Pointer_Fast(iList, i + 3 * nbNod +
- iNbComp * nbNod * iTS);
- double *dVal = (double *)List_Pointer_Fast(dList, j + 3 * nbNod +
- dNbComp * nbNod * dTS);
- if(zeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z, iVal, iNbComp,
- NULL, 0, out)) {
- for(int k = 0; k < iDec.numSimplices(); k++) {
- iDec.decompose(k, x, y, z, iVal, xNew, yNew, zNew, iValNew);
- dDec.decompose(k, x, y, z, dVal, xNew, yNew, zNew, dValNew);
- zeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri,
- xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
- if(_extractVolume)
- nonZeroLevelset(iTS, nbNodNew, nbEdgNew, (nbNodNew == 4) ? exnTet : exnTri,
- xNew, yNew, zNew, iValNew, iNbComp, dValNew, dNbComp, out);
- }
- }
- else if(_extractVolume)
- nonZeroLevelset(iTS, nbNod, nbEdg, exn, x, y, z,
- iVal, iNbComp, dVal, dNbComp, out);
- }
- }
- }
-
- delete [] iValNew;
- delete [] dValNew;
+ _addElement(step, numNodes, numEdges, numComp, xp, yp, zp, valp, out);
}
-
}
}
-
+
PView *GMSH_LevelsetPlugin::execute(PView *v)
{
- PView *w;
- std::vector<PViewDataList*> out;
-
- PViewDataList *dv = getDataList(v);
- if(!dv) return v;
-
- if(dv->isAdaptive()){
- dv->adaptive->setTolerance(_targetError);
- if(dv->NbST || dv->NbSS || dv->NbSQ || dv->NbSH){
- dv->adaptive->setAdaptiveResolutionLevel(dv, _recurLevel, this);
- v->setChanged(true);
+ if(v->getData()->isAdaptive()){
+ PViewDataList *dv = getDataList(v);
+ if(dv){
+ dv->adaptive->setTolerance(_targetError);
+ if(dv->NbST || dv->NbSS || dv->NbSQ || dv->NbSH){
+ dv->adaptive->setAdaptiveResolutionLevel(dv, _recurLevel, this);
+ v->setChanged(true);
+ }
}
}
+ PViewData *vdata = v->getData(), *wdata;
if(_valueView < 0) {
- w = v;
+ wdata = vdata;
}
else if(_valueView > (int)PView::list.size() - 1){
Msg(GERROR, "View[%d] does not exist: reverting to View[%d]",
_valueView, v->getIndex());
- w = v;
+ wdata = vdata;
}
else{
- w = PView::list[_valueView];
+ wdata = PView::list[_valueView]->getData();
}
- PViewDataList *dw = getDataList(w);
- if(!dw) return v;
+ // sanity checks
+ if(vdata->getNumEntities() != wdata->getNumEntities() ||
+ vdata->getNumElements() != wdata->getNumElements()) return v;
+ if(_valueTimeStep >= wdata->getNumTimeSteps()) {
+ Msg(GERROR, "Wrong time step %d in view", _valueTimeStep);
+ return v;
+ }
+ // create output dataset(s), in the list format
+ std::vector<PViewDataList*> out;
if(_valueIndependent) {
out.push_back(getDataList(new PView(true)));
}
else{
- for(int ts = 0; ts < dv->getNumTimeSteps(); ts++)
+ for(int ts = 0; ts < vdata->getNumTimeSteps(); ts++)
out.push_back(getDataList(new PView(true)));
}
-
- // We should definitely recode the View interface in C++ (and define
- // iterators on the different kinds of elements...)
-
- // To avoid surprising results, we don't try to plot values of
- // different types on the levelset if the value view (w) is the same
- // as the levelset view (v).
-
- // points
- executeList(dv, dv->SP, dv->NbSP, 1, dw, dw->SP, dw->NbSP, 1, 1, 0, 0, out);
- if(dv != dw) {
- executeList(dv, dv->SP, dv->NbSP, 1, dw, dw->VP, dw->NbVP, 3, 1, 0, 0, out);
- executeList(dv, dv->SP, dv->NbSP, 1, dw, dw->TP, dw->NbTP, 9, 1, 0, 0, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VP, dv->NbVP, 3, dw, dw->SP, dw->NbSP, 1, 1, 0, 0, out);
- }
- executeList(dv, dv->VP, dv->NbVP, 3, dw, dw->VP, dw->NbVP, 3, 1, 0, 0, out);
- if(dv != dw) {
- executeList(dv, dv->VP, dv->NbVP, 3, dw, dw->TP, dw->NbTP, 9, 1, 0, 0, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TP, dv->NbTP, 9, dw, dw->SP, dw->NbSP, 1, 1, 0, 0, out);
- executeList(dv, dv->TP, dv->NbTP, 9, dw, dw->VP, dw->NbVP, 3, 1, 0, 0, out);
- }
- executeList(dv, dv->TP, dv->NbTP, 9, dw, dw->TP, dw->NbTP, 9, 1, 0, 0, out);
-
- // lines
- int exnLin[12][2] = {{0,1}};
- executeList(dv, dv->SL, dv->NbSL, 1, dw, dw->SL, dw->NbSL, 1, 2, 1, exnLin, out);
- if(dv != dw) {
- executeList(dv, dv->SL, dv->NbSL, 1, dw, dw->VL, dw->NbVL, 3, 2, 1, exnLin, out);
- executeList(dv, dv->SL, dv->NbSL, 1, dw, dw->TL, dw->NbTL, 9, 2, 1, exnLin, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VL, dv->NbVL, 3, dw, dw->SL, dw->NbSL, 1, 2, 1, exnLin, out);
- }
- executeList(dv, dv->VL, dv->NbVL, 3, dw, dw->VL, dw->NbVL, 3, 2, 1, exnLin, out);
- if(dv != dw) {
- executeList(dv, dv->VL, dv->NbVL, 3, dw, dw->TL, dw->NbTL, 9, 2, 1, exnLin, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TL, dv->NbTL, 9, dw, dw->SL, dw->NbSL, 1, 2, 1, exnLin, out);
- executeList(dv, dv->TL, dv->NbTL, 9, dw, dw->VL, dw->NbVL, 3, 2, 1, exnLin, out);
- }
- executeList(dv, dv->TL, dv->NbTL, 9, dw, dw->TL, dw->NbTL, 9, 2, 1, exnLin, out);
-
- // triangles
- int exnTri[12][2] = {{0,1}, {0,2}, {1,2}};
- executeList(dv, dv->ST, dv->NbST, 1, dw, dw->ST, dw->NbST, 1, 3, 3, exnTri, out);
- if(dv != dw) {
- executeList(dv, dv->ST, dv->NbST, 1, dw, dw->VT, dw->NbVT, 3, 3, 3, exnTri, out);
- executeList(dv, dv->ST, dv->NbST, 1, dw, dw->TT, dw->NbTT, 9, 3, 3, exnTri, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VT, dv->NbVT, 3, dw, dw->ST, dw->NbST, 1, 3, 3, exnTri, out);
- }
- executeList(dv, dv->VT, dv->NbVT, 3, dw, dw->VT, dw->NbVT, 3, 3, 3, exnTri, out);
- if(dv != dw) {
- executeList(dv, dv->VT, dv->NbVT, 3, dw, dw->TT, dw->NbTT, 9, 3, 3, exnTri, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TT, dv->NbTT, 9, dw, dw->ST, dw->NbST, 1, 3, 3, exnTri, out);
- executeList(dv, dv->TT, dv->NbTT, 9, dw, dw->VT, dw->NbVT, 3, 3, 3, exnTri, out);
- }
- executeList(dv, dv->TT, dv->NbTT, 9, dw, dw->TT, dw->NbTT, 9, 3, 3, exnTri, out);
-
- // tets
- int exnTet[12][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
- executeList(dv, dv->SS, dv->NbSS, 1, dw, dw->SS, dw->NbSS, 1, 4, 6, exnTet, out);
- if(dv != dw) {
- executeList(dv, dv->SS, dv->NbSS, 1, dw, dw->VS, dw->NbVS, 3, 4, 6, exnTet, out);
- executeList(dv, dv->SS, dv->NbSS, 1, dw, dw->TS, dw->NbTS, 9, 4, 6, exnTet, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VS, dv->NbVS, 3, dw, dw->SS, dw->NbSS, 1, 4, 6, exnTet, out);
- }
- executeList(dv, dv->VS, dv->NbVS, 3, dw, dw->VS, dw->NbVS, 3, 4, 6, exnTet, out);
- if(dv != dw) {
- executeList(dv, dv->VS, dv->NbVS, 3, dw, dw->TS, dw->NbTS, 9, 4, 6, exnTet, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TS, dv->NbTS, 9, dw, dw->SS, dw->NbSS, 1, 4, 6, exnTet, out);
- executeList(dv, dv->TS, dv->NbTS, 9, dw, dw->VS, dw->NbVS, 3, 4, 6, exnTet, out);
- }
- executeList(dv, dv->TS, dv->NbTS, 9, dw, dw->TS, dw->NbTS, 9, 4, 6, exnTet, out);
-
- // quads
- int exnQua[12][2] = {{0,1}, {0,3}, {1,2}, {2,3}};
- executeList(dv, dv->SQ, dv->NbSQ, 1, dw, dw->SQ, dw->NbSQ, 1, 4, 4, exnQua, out);
- if(dv != dw) {
- executeList(dv, dv->SQ, dv->NbSQ, 1, dw, dw->VQ, dw->NbVQ, 3, 4, 4, exnQua, out);
- executeList(dv, dv->SQ, dv->NbSQ, 1, dw, dw->TQ, dw->NbTQ, 9, 4, 4, exnQua, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VQ, dv->NbVQ, 3, dw, dw->SQ, dw->NbSQ, 1, 4, 4, exnQua, out);
- }
- executeList(dv, dv->VQ, dv->NbVQ, 3, dw, dw->VQ, dw->NbVQ, 3, 4, 4, exnQua, out);
- if(dv != dw) {
- executeList(dv, dv->VQ, dv->NbVQ, 3, dw, dw->TQ, dw->NbTQ, 9, 4, 4, exnQua, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TQ, dv->NbTQ, 9, dw, dw->SQ, dw->NbSQ, 1, 4, 4, exnQua, out);
- executeList(dv, dv->TQ, dv->NbTQ, 9, dw, dw->VQ, dw->NbVQ, 3, 4, 4, exnQua, out);
- }
- executeList(dv, dv->TQ, dv->NbTQ, 9, dw, dw->TQ, dw->NbTQ, 9, 4, 4, exnQua, out);
-
- // hexes
- int exnHex[12][2] = {{0,1}, {0,3}, {0,4}, {1,2}, {1,5}, {2,3},
- {2,6}, {3,7}, {4,5}, {4,7}, {5,6}, {6,7}};
- executeList(dv, dv->SH, dv->NbSH, 1, dw, dw->SH, dw->NbSH, 1, 8, 12, exnHex, out);
- if(dv != dw) {
- executeList(dv, dv->SH, dv->NbSH, 1, dw, dw->VH, dw->NbVH, 3, 8, 12, exnHex, out);
- executeList(dv, dv->SH, dv->NbSH, 1, dw, dw->TH, dw->NbTH, 9, 8, 12, exnHex, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VH, dv->NbVH, 3, dw, dw->SH, dw->NbSH, 1, 8, 12, exnHex, out);
- }
- executeList(dv, dv->VH, dv->NbVH, 3, dw, dw->VH, dw->NbVH, 3, 8, 12, exnHex, out);
- if(dv != dw) {
- executeList(dv, dv->VH, dv->NbVH, 3, dw, dw->TH, dw->NbTH, 9, 8, 12, exnHex, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TH, dv->NbTH, 9, dw, dw->SH, dw->NbSH, 1, 8, 12, exnHex, out);
- executeList(dv, dv->TH, dv->NbTH, 9, dw, dw->VH, dw->NbVH, 3, 8, 12, exnHex, out);
- }
- executeList(dv, dv->TH, dv->NbTH, 9, dw, dw->TH, dw->NbTH, 9, 8, 12, exnHex, out);
-
- // prisms
- int exnPri[12][2] = {{0,1}, {0,2}, {0,3}, {1,2}, {1,4}, {2,5},
- {3,4}, {3,5}, {4,5}};
- executeList(dv, dv->SI, dv->NbSI, 1, dw, dw->SI, dw->NbSI, 1, 6, 9, exnPri, out);
- if(dv != dw) {
- executeList(dv, dv->SI, dv->NbSI, 1, dw, dw->VI, dw->NbVI, 3, 6, 9, exnPri, out);
- executeList(dv, dv->SI, dv->NbSI, 1, dw, dw->TI, dw->NbTI, 9, 6, 9, exnPri, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VI, dv->NbVI, 3, dw, dw->SI, dw->NbSI, 1, 6, 9, exnPri, out);
- }
- executeList(dv, dv->VI, dv->NbVI, 3, dw, dw->VI, dw->NbVI, 3, 6, 9, exnPri, out);
- if(dv != dw) {
- executeList(dv, dv->VI, dv->NbVI, 3, dw, dw->TI, dw->NbTI, 9, 6, 9, exnPri, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TI, dv->NbTI, 9, dw, dw->SI, dw->NbSI, 1, 6, 9, exnPri, out);
- executeList(dv, dv->TI, dv->NbTI, 9, dw, dw->VI, dw->NbVI, 3, 6, 9, exnPri, out);
- }
- executeList(dv, dv->TI, dv->NbTI, 9, dw, dw->TI, dw->NbTI, 9, 6, 9, exnPri, out);
-
- // pyramids
- int exnPyr[12][2] = {{0,1}, {0,3}, {0,4}, {1,2}, {1,4}, {2,3}, {2,4}, {3,4}};
- executeList(dv, dv->SY, dv->NbSY, 1, dw, dw->SY, dw->NbSY, 1, 5, 8, exnPyr, out);
- if(dv != dw) {
- executeList(dv, dv->SY, dv->NbSY, 1, dw, dw->VY, dw->NbVY, 3, 5, 8, exnPyr, out);
- executeList(dv, dv->SY, dv->NbSY, 1, dw, dw->TY, dw->NbTY, 9, 5, 8, exnPyr, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->VY, dv->NbVY, 3, dw, dw->SY, dw->NbSY, 1, 5, 8, exnPyr, out);
- }
- executeList(dv, dv->VY, dv->NbVY, 3, dw, dw->VY, dw->NbVY, 3, 5, 8, exnPyr, out);
- if(dv != dw) {
- executeList(dv, dv->VY, dv->NbVY, 3, dw, dw->TY, dw->NbTY, 9, 5, 8, exnPyr, out);
- }
-
- if(dv != dw) {
- executeList(dv, dv->TY, dv->NbTY, 9, dw, dw->SY, dw->NbSY, 1, 5, 8, exnPyr, out);
- executeList(dv, dv->TY, dv->NbTY, 9, dw, dw->VY, dw->NbVY, 3, 5, 8, exnPyr, out);
+
+ for(int ent = 0; ent < vdata->getNumEntities(); ent++){
+ for(int ele = 0; ele < vdata->getNumElements(ent); ele++){
+ int numNodes = vdata->getNumNodes(ent, ele);
+ int numEdges = vdata->getNumEdges(ent, ele);
+ double x[8], y[8], z[8], levels[8];
+ double scalarValues[8] = {0., 0., 0., 0., 0., 0., 0., 0.};
+ for(int nod = 0; nod < numNodes; nod++){
+ vdata->getNode(ent, ele, nod, x[nod], y[nod], z[nod]);
+ if(_valueIndependent)
+ levels[nod] = levelset(x[nod], y[nod], z[nod], 0.);
+ }
+ for(int step = 0; step < vdata->getNumTimeSteps(); step++){
+ if(!_valueIndependent){
+ for(int nod = 0; nod < numNodes; nod++){
+ vdata->getScalarValue(ent, ele, nod, step, scalarValues[nod]);
+ levels[nod] = levelset(x[nod], y[nod], z[nod], scalarValues[nod]);
+ }
+ }
+ int wstep = (_valueTimeStep < 0) ? step : _valueTimeStep;
+ _cutAndAddElements(vdata, wdata, ent, ele, step, wstep, x, y, z,
+ levels, scalarValues, out);
+ }
+ }
}
- executeList(dv, dv->TY, dv->NbTY, 9, dw, dw->TY, dw->NbTY, 9, 5, 8, exnPyr, out);
-
+
for(unsigned int i = 0; i < out.size(); i++) {
- // FIXME: create time data
char tmp[246];
sprintf(tmp, "_Levelset_%d", i);
- out[i]->setName(dv->getName() + tmp);
- out[i]->setFileName(dv->getFileName() + tmp + ".pos");
+ out[i]->setName(vdata->getName() + tmp);
+ out[i]->setFileName(vdata->getFileName() + tmp + ".pos");
out[i]->finalize();
}
@@ -846,7 +561,7 @@ static bool recur_sign_change(adapt_tet *t, double val,
static bool recur_sign_change(adapt_hex *t, double val,
const GMSH_LevelsetPlugin *plug)
{
- if (!t->e[0]|| t->visible){
+ if (!t->e[0] || t->visible){
double v1 = plug->levelset(t->p[0]->X, t->p[0]->Y, t->p[0]->Z, t->p[0]->val);
double v2 = plug->levelset(t->p[1]->X, t->p[1]->Y, t->p[1]->Z, t->p[1]->val);
double v3 = plug->levelset(t->p[2]->X, t->p[2]->Y, t->p[2]->Z, t->p[2]->val);
@@ -887,8 +602,8 @@ static bool recur_sign_change(adapt_hex *t, double val,
}
}
-static bool recur_sign_change (adapt_quad *q, double val,
- const GMSH_LevelsetPlugin *plug)
+static bool recur_sign_change(adapt_quad *q, double val,
+ const GMSH_LevelsetPlugin *plug)
{
if(!q->e[0]|| q->visible){
double v1 = plug->levelset(q->p[0]->X, q->p[0]->Y, q->p[0]->Z, q->p[0]->val);
diff --git a/Plugin/Levelset.h b/Plugin/Levelset.h
index 7fa6257e14f39b185772f19db5762d05c07737e1..1575e9b149e71c6ae79fca7858308d5caababb3c 100644
--- a/Plugin/Levelset.h
+++ b/Plugin/Levelset.h
@@ -26,39 +26,26 @@
class GMSH_LevelsetPlugin : public GMSH_Post_Plugin
{
-public:
- typedef enum {NONE, PLANE, SPHERE, MAP} ORIENTATION ;
- virtual double levelset(double x, double y, double z, double val) const = 0;
-protected:
+ private:
+ double _invert;
+ void _addElement(int step, int np, int numEdges, int numComp,
+ double xp[12], double yp[12], double zp[12],
+ double valp[12][9], std::vector<PViewDataList*> &out);
+ void _cutAndAddElements(PViewData *vdata, PViewData *wdata,
+ int ent, int ele, int step, int wstep,
+ double x[8], double y[8], double z[8],
+ double levels[8], double scalarValues[8],
+ std::vector<PViewDataList*> &out);
+ protected:
double _ref[3], _targetError;
int _valueTimeStep, _valueView, _valueIndependent, _recurLevel, _extractVolume;
+ typedef enum {NONE, PLANE, SPHERE, MAP} ORIENTATION;
ORIENTATION _orientation;
-private:
- double _invert;
- void addElement(int timeStep, int np, int nbEdg, int dNbComp,
- double xp[12], double yp[12], double zp[12],
- double valp[12][9], std::vector<PViewDataList*> &out);
- void evalLevelset(int nbNod, int nbComp,
- double *x, double *y, double *z, double *val,
- double *levels, double *scalarVal);
- void nonZeroLevelset(int timeStep, int nbVert, int nbEdg, int exn[12][2],
- double *x, double *y, double *z,
- double *iVal, int iNbComp, double *dVal, int dNbComp,
- std::vector<PViewDataList*> &out);
- int zeroLevelset(int timeStep, int nbVert, int nbEdg, int exn[12][2],
- double *x, double *y, double *z,
- double *iVal, int iNbComp, double *dVal, int dNbComp,
- std::vector<PViewDataList*> &out);
- void executeList(PViewDataList *iData, List_T *iList,
- int iNbElm, int iNbComp,
- PViewDataList *dData, List_T *dList,
- int dNbElm, int dNbComp,
- int nbVert, int nbEdg, int exn[12][2],
- std::vector<PViewDataList*> &out);
- virtual void assignSpecificVisibility () const;
public:
GMSH_LevelsetPlugin();
+ virtual double levelset(double x, double y, double z, double val) const = 0;
virtual PView *execute(PView *);
+ void assignSpecificVisibility() const;
};
#endif
diff --git a/Plugin/MakeSimplex.cpp b/Plugin/MakeSimplex.cpp
index a00c4576ecb2efdf79c79078bb08dc15142c0994..7831612b97e214ac4d2fa83014acff063e50afde 100644
--- a/Plugin/MakeSimplex.cpp
+++ b/Plugin/MakeSimplex.cpp
@@ -1,4 +1,4 @@
-// $Id: MakeSimplex.cpp,v 1.6 2008-02-17 08:48:07 geuzaine Exp $
+// $Id: MakeSimplex.cpp,v 1.7 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -184,26 +184,6 @@ void MakeSimplex::reorder(int map[4], int n,
}
int map2[4] = {map[0], map[1], map[2], map[3]};
-#if 0
- // make sure to create tets with positive volume?
- if(n == 4){ // tets
- double mat[3][3];
- mat[0][0] = xn[1] - xn[0]; mat[0][1] = xn[2] - xn[0]; mat[0][2] = xn[3] - xn[0];
- mat[1][0] = yn[1] - yn[0]; mat[1][1] = yn[2] - yn[0]; mat[1][2] = yn[3] - yn[0];
- mat[2][0] = zn[1] - zn[0]; mat[2][1] = zn[2] - zn[0]; mat[2][2] = zn[3] - zn[0];
- if(det3x3(mat) < 0.){
- map2[0] = map[1];
- map2[1] = map[0];
- xn[0] = x[map2[0]];
- yn[0] = y[map2[0]];
- zn[0] = z[map2[0]];
- xn[1] = x[map2[1]];
- yn[1] = y[map2[1]];
- zn[1] = z[map2[1]];
- }
- }
-#endif
-
for(int ts = 0; ts < _numTimeSteps; ts++)
for(int i = 0; i < n; i++) {
for(int j = 0; j < _numComponents; j++)
diff --git a/Post/OctreePost.cpp b/Post/OctreePost.cpp
index 97f913d9129bf319727161c61082ddd584ce18e0..1b8e745a64acb87efd3a3d301216cb4a2d9a4371 100644
--- a/Post/OctreePost.cpp
+++ b/Post/OctreePost.cpp
@@ -1,4 +1,4 @@
-// $Id: OctreePost.cpp,v 1.3 2008-02-17 08:48:08 geuzaine Exp $
+// $Id: OctreePost.cpp,v 1.4 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -24,6 +24,7 @@
#include "List.h"
#include "PView.h"
#include "PViewDataList.h"
+#include "PViewDataGModel.h"
#include "Numeric.h"
#include "Message.h"
#include "ShapeFunctions.h"
@@ -216,107 +217,120 @@ static void addListOfStuff(Octree *o, List_T *l, int nbelm)
OctreePost::~OctreePost()
{
- Octree_Delete(ST); Octree_Delete(VT); Octree_Delete(TT);
- Octree_Delete(SQ); Octree_Delete(VQ); Octree_Delete(TQ);
- Octree_Delete(SS); Octree_Delete(VS); Octree_Delete(TS);
- Octree_Delete(SH); Octree_Delete(VH); Octree_Delete(TH);
- Octree_Delete(SI); Octree_Delete(VI); Octree_Delete(TI);
- Octree_Delete(SY); Octree_Delete(VY); Octree_Delete(TY);
+ if(_ST) Octree_Delete(_ST); if(_VT) Octree_Delete(_VT); if(_TT) Octree_Delete(_TT);
+ if(_SQ) Octree_Delete(_SQ); if(_VQ) Octree_Delete(_VQ); if(_TQ) Octree_Delete(_TQ);
+ if(_SS) Octree_Delete(_SS); if(_VS) Octree_Delete(_VS); if(_TS) Octree_Delete(_TS);
+ if(_SH) Octree_Delete(_SH); if(_VH) Octree_Delete(_VH); if(_TH) Octree_Delete(_TH);
+ if(_SI) Octree_Delete(_SI); if(_VI) Octree_Delete(_VI); if(_TI) Octree_Delete(_TI);
+ if(_SY) Octree_Delete(_SY); if(_VY) Octree_Delete(_VY); if(_TY) Octree_Delete(_TY);
+ if(_GModel) Octree_Delete(_GModel);
}
OctreePost::OctreePost(PView *v)
- : theView(v)
+ : _ST(0), _VT(0), _TT(0), _SQ(0), _VQ(0), _TQ(0), _SS(0), _VS(0), _TS(0),
+ _SH(0), _VH(0), _TH(0), _SI(0), _VI(0), _TI(0), _SY(0), _VY(0), _TY(0),
+ _GModel(0), _theView(v), _viewType(-1)
{
- PViewDataList *l = dynamic_cast<PViewDataList*>(theView->getData());
-
- if(!l){
- Msg(GERROR, "OctreePost only accepts list-based view data for now");
- return;
- }
-
- SBoundingBox3d bb = l->getBoundingBox();
+ SBoundingBox3d bb = v->getData()->getBoundingBox();
double min[3] = {bb.min().x(), bb.min().y(), bb.min().z()};
double size[3] = {bb.max().x() - bb.min().x(),
bb.max().y() - bb.min().y(),
bb.max().z() - bb.min().z()};
-
+
const int maxElePerBucket = 100; // memory vs. speed trade-off
- SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
- addListOfStuff(SL, l->SL, 6 + 2 * l->getNumTimeSteps());
- Octree_Arrange(SL);
- VL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
- addListOfStuff(VL, l->VL, 6 + 6 * l->getNumTimeSteps());
- Octree_Arrange(VL);
- TL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
- addListOfStuff(TL, l->TL, 6 + 18 * l->getNumTimeSteps());
- Octree_Arrange(TL);
-
- ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
- addListOfStuff(ST, l->ST, 9 + 3 * l->getNumTimeSteps());
- Octree_Arrange(ST);
- VT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
- addListOfStuff(VT, l->VT, 9 + 9 * l->getNumTimeSteps());
- Octree_Arrange(VT);
- TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
- addListOfStuff(TT, l->TT, 9 + 27 * l->getNumTimeSteps());
- Octree_Arrange(TT);
-
- SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
- addListOfStuff(SQ, l->SQ, 12 + 4 * l->getNumTimeSteps());
- Octree_Arrange(SQ);
- VQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
- addListOfStuff(VQ, l->VQ, 12 + 12 * l->getNumTimeSteps());
- Octree_Arrange(VQ);
- TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
- addListOfStuff(TQ, l->TQ, 12 + 36 * l->getNumTimeSteps());
- Octree_Arrange(TQ);
-
- SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
- addListOfStuff(SS, l->SS, 12 + 4 * l->getNumTimeSteps());
- Octree_Arrange(SS);
- VS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
- addListOfStuff(VS, l->VS, 12 + 12 * l->getNumTimeSteps());
- Octree_Arrange(VS);
- TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
- addListOfStuff(TS, l->TS, 12 + 36 * l->getNumTimeSteps());
- Octree_Arrange(TS);
-
- SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
- addListOfStuff(SH, l->SH, 24 + 8 * l->getNumTimeSteps());
- Octree_Arrange(SH);
- VH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
- addListOfStuff(VH, l->VH, 24 + 24 * l->getNumTimeSteps());
- Octree_Arrange(VH);
- TH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
- addListOfStuff(TH, l->TH, 24 + 72 * l->getNumTimeSteps());
- Octree_Arrange(TH);
-
- SI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
- addListOfStuff(SI, l->SI, 18 + 6 * l->getNumTimeSteps());
- Octree_Arrange(SI);
- VI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
- addListOfStuff(VI, l->VI, 18 + 18 * l->getNumTimeSteps());
- Octree_Arrange(VI);
- TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
- addListOfStuff(TI, l->TI, 18 + 54 * l->getNumTimeSteps());
- Octree_Arrange(TI);
-
- SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
- addListOfStuff(SY, l->SY, 15 + 5 * l->getNumTimeSteps());
- Octree_Arrange(SY);
- VY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
- addListOfStuff(VY, l->VY, 15 + 15 * l->getNumTimeSteps());
- Octree_Arrange(VY);
- TY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
- addListOfStuff(TY, l->TY, 15 + 45 * l->getNumTimeSteps());
- Octree_Arrange(TY);
+ PViewDataList *l = dynamic_cast<PViewDataList*>(_theView->getData());
+ if(l){
+ _viewType = 0;
+
+ _SL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(_SL, l->SL, 6 + 2 * l->getNumTimeSteps());
+ Octree_Arrange(_SL);
+ _VL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(_VL, l->VL, 6 + 6 * l->getNumTimeSteps());
+ Octree_Arrange(_VL);
+ _TL = Octree_Create(maxElePerBucket, min, size, linBB, linCentroid, linInEle);
+ addListOfStuff(_TL, l->TL, 6 + 18 * l->getNumTimeSteps());
+ Octree_Arrange(_TL);
+
+ _ST = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(_ST, l->ST, 9 + 3 * l->getNumTimeSteps());
+ Octree_Arrange(_ST);
+ _VT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(_VT, l->VT, 9 + 9 * l->getNumTimeSteps());
+ Octree_Arrange(_VT);
+ _TT = Octree_Create(maxElePerBucket, min, size, triBB, triCentroid, triInEle);
+ addListOfStuff(_TT, l->TT, 9 + 27 * l->getNumTimeSteps());
+ Octree_Arrange(_TT);
+
+ _SQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(_SQ, l->SQ, 12 + 4 * l->getNumTimeSteps());
+ Octree_Arrange(_SQ);
+ _VQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(_VQ, l->VQ, 12 + 12 * l->getNumTimeSteps());
+ Octree_Arrange(_VQ);
+ _TQ = Octree_Create(maxElePerBucket, min, size, quaBB, quaCentroid, quaInEle);
+ addListOfStuff(_TQ, l->TQ, 12 + 36 * l->getNumTimeSteps());
+ Octree_Arrange(_TQ);
+
+ _SS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(_SS, l->SS, 12 + 4 * l->getNumTimeSteps());
+ Octree_Arrange(_SS);
+ _VS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(_VS, l->VS, 12 + 12 * l->getNumTimeSteps());
+ Octree_Arrange(_VS);
+ _TS = Octree_Create(maxElePerBucket, min, size, tetBB, tetCentroid, tetInEle);
+ addListOfStuff(_TS, l->TS, 12 + 36 * l->getNumTimeSteps());
+ Octree_Arrange(_TS);
+
+ _SH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(_SH, l->SH, 24 + 8 * l->getNumTimeSteps());
+ Octree_Arrange(_SH);
+ _VH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(_VH, l->VH, 24 + 24 * l->getNumTimeSteps());
+ Octree_Arrange(_VH);
+ _TH = Octree_Create(maxElePerBucket, min, size, hexBB, hexCentroid, hexInEle);
+ addListOfStuff(_TH, l->TH, 24 + 72 * l->getNumTimeSteps());
+ Octree_Arrange(_TH);
+
+ _SI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(_SI, l->SI, 18 + 6 * l->getNumTimeSteps());
+ Octree_Arrange(_SI);
+ _VI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(_VI, l->VI, 18 + 18 * l->getNumTimeSteps());
+ Octree_Arrange(_VI);
+ _TI = Octree_Create(maxElePerBucket, min, size, priBB, priCentroid, priInEle);
+ addListOfStuff(_TI, l->TI, 18 + 54 * l->getNumTimeSteps());
+ Octree_Arrange(_TI);
+
+ _SY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(_SY, l->SY, 15 + 5 * l->getNumTimeSteps());
+ Octree_Arrange(_SY);
+ _VY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(_VY, l->VY, 15 + 15 * l->getNumTimeSteps());
+ Octree_Arrange(_VY);
+ _TY = Octree_Create(maxElePerBucket, min, size, pyrBB, pyrCentroid, pyrInEle);
+ addListOfStuff(_TY, l->TY, 15 + 45 * l->getNumTimeSteps());
+ Octree_Arrange(_TY);
+ }
+
+ PViewDataGModel *g = dynamic_cast<PViewDataGModel*>(_theView->getData());
+ if(g){
+ _viewType = 1;
+
+ //_GModel = Octree_Create(maxElePerBucket, min, size,
+ // MElementBB, MElementCentroid, MElementInEle);
+ // add MElement pointers in
+ //Octree_Arrange(_GModel);
+ }
+
}
-bool OctreePost::getValue(void *in, int dim, int nbNod, int nbComp,
- double P[3], int timestep, double *values, double *size_elem)
+bool OctreePost::_getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int step, double *values,
+ double *elementSize)
{
if(!in) return false;
@@ -327,89 +341,106 @@ bool OctreePost::getValue(void *in, int dim, int nbNod, int nbComp,
if(!e) return false;
e->xyz2uvw(P, U);
- if(timestep < 0){
- for(int i = 0; i < theView->getData()->getNumTimeSteps(); i++)
+ if(step < 0){
+ for(int i = 0; i < _theView->getData()->getNumTimeSteps(); i++)
for(int j = 0; j < nbComp; j++)
values[nbComp * i + j] = e->interpolate(&V[nbNod * nbComp * i + j],
U[0], U[1], U[2], nbComp);
}
else{
for(int j = 0; j < nbComp; j++)
- values[j] = e->interpolate(&V[nbNod * nbComp * timestep + j],
+ values[j] = e->interpolate(&V[nbNod * nbComp * step + j],
U[0], U[1], U[2], nbComp);
}
- if(size_elem)
- *size_elem = e->maxEdgeLength();
+ if(elementSize) *elementSize = e->maxEdgeLength();
delete e;
return true;
}
bool OctreePost::searchScalar(double x, double y, double z, double *values,
- int timestep, double *size_elem)
+ int step, double *size)
{
double P[3] = {x, y, z};
- if(timestep < 0)
- for(int i = 0; i < theView->getData()->getNumTimeSteps(); i++)
+ if(step < 0)
+ for(int i = 0; i < _theView->getData()->getNumTimeSteps(); i++)
values[i] = 0.0;
else
values[0] = 0.0;
- if(getValue(Octree_Search(P, SS), 3, 4, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, SH), 3, 8, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, SI), 3, 6, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, SY), 3, 5, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, ST), 2, 3, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, SQ), 2, 4, 1, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, SL), 1, 2, 1, P, timestep, values, size_elem)) return true;
+ if(_viewType == 0){
+ if(_getValue(Octree_Search(P, _SS), 3, 4, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _SH), 3, 8, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _SI), 3, 6, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _SY), 3, 5, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _ST), 2, 3, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _SQ), 2, 4, 1, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _SL), 1, 2, 1, P, step, values, size)) return true;
+ }
+
+ if(_viewType == 1){
+
+ }
return false;
}
bool OctreePost::searchVector(double x, double y, double z, double *values,
- int timestep, double *size_elem)
+ int step, double *size)
{
double P[3] = {x, y, z};
- if(timestep < 0)
- for(int i = 0; i < 3 * theView->getData()->getNumTimeSteps(); i++)
+ if(step < 0)
+ for(int i = 0; i < 3 * _theView->getData()->getNumTimeSteps(); i++)
values[i] = 0.0;
else
for(int i = 0; i < 3; i++)
values[i] = 0.0;
- if(getValue(Octree_Search(P, VS), 3, 4, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VH), 3, 8, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VI), 3, 6, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VY), 3, 5, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VT), 2, 3, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VQ), 2, 4, 3, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, VL), 1, 2, 3, P, timestep, values, size_elem)) return true;
+ if(_viewType == 0){
+ if(_getValue(Octree_Search(P, _VS), 3, 4, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VH), 3, 8, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VI), 3, 6, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VY), 3, 5, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VT), 2, 3, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VQ), 2, 4, 3, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _VL), 1, 2, 3, P, step, values, size)) return true;
+ }
+
+ if(_viewType == 1){
+
+ }
return false;
}
bool OctreePost::searchTensor(double x, double y, double z, double *values,
- int timestep, double *size_elem)
+ int step, double *size)
{
double P[3] = {x, y, z};
- if(timestep < 0)
- for(int i = 0; i < 9 * theView->getData()->getNumTimeSteps(); i++)
+ if(step < 0)
+ for(int i = 0; i < 9 * _theView->getData()->getNumTimeSteps(); i++)
values[i] = 0.0;
else
for(int i = 0; i < 9; i++)
values[i] = 0.0;
- if(getValue(Octree_Search(P, TS), 3, 4, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TH), 3, 8, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TI), 3, 6, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TY), 3, 5, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TT), 2, 3, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TQ), 2, 4, 9, P, timestep, values, size_elem)) return true;
- if(getValue(Octree_Search(P, TL), 1, 2, 9, P, timestep, values, size_elem)) return true;
+ if(_viewType == 0){
+ if(_getValue(Octree_Search(P, _TS), 3, 4, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TH), 3, 8, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TI), 3, 6, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TY), 3, 5, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TT), 2, 3, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TQ), 2, 4, 9, P, step, values, size)) return true;
+ if(_getValue(Octree_Search(P, _TL), 1, 2, 9, P, step, values, size)) return true;
+ }
+
+ if(_viewType == 1){
+
+ }
return false;
}
diff --git a/Post/OctreePost.h b/Post/OctreePost.h
index ec7301751f6afd83c351b6f63821acc947fccca5..b1e0911d245097f863c22fa41af0f44a898bccab 100644
--- a/Post/OctreePost.h
+++ b/Post/OctreePost.h
@@ -26,16 +26,20 @@ class PView;
class OctreePost
{
- Octree *SL, *VL, *TL;
- Octree *ST, *VT, *TT;
- Octree *SQ, *VQ, *TQ;
- Octree *SS, *VS, *TS;
- Octree *SH, *VH, *TH;
- Octree *SI, *VI, *TI;
- Octree *SY, *VY, *TY;
- PView *theView;
- bool getValue(void *in, int dim, int nbNod, int nbComp,
- double P[3], int timestep, double *values, double *size_elem);
+ private:
+ Octree *_SL, *_VL, *_TL;
+ Octree *_ST, *_VT, *_TT;
+ Octree *_SQ, *_VQ, *_TQ;
+ Octree *_SS, *_VS, *_TS;
+ Octree *_SH, *_VH, *_TH;
+ Octree *_SI, *_VI, *_TI;
+ Octree *_SY, *_VY, *_TY;
+ Octree *_GModel;
+ PView *_theView;
+ int _viewType; // internal view type (0=list, 1=GModel)
+ bool _getValue(void *in, int dim, int nbNod, int nbComp,
+ double P[3], int step, double *values,
+ double *elementSize);
public :
OctreePost(PView *);
~OctreePost();
@@ -45,11 +49,11 @@ class OctreePost
// interpolated unless time step is set to a different value than
// -1.
bool searchScalar(double x, double y, double z, double *values,
- int timestep = -1, double *size_elem = 0);
+ int step = -1, double *size = 0);
bool searchVector(double x, double y, double z, double *values,
- int timestep = -1, double *size_elem = 0);
+ int step = -1, double *size = 0);
bool searchTensor(double x, double y, double z, double *values,
- int timestep = -1, double *size_elem = 0);
+ int step = -1, double *size = 0);
};
#endif
diff --git a/Post/PView.cpp b/Post/PView.cpp
index 626cb03631c8305f7a95c1fa979869ab4d403faa..5100c4b73c2054b9766d1acfc20338d9ba5a9a5d 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -1,4 +1,4 @@
-// $Id: PView.cpp,v 1.22 2008-03-12 21:28:53 geuzaine Exp $
+// $Id: PView.cpp,v 1.23 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -231,16 +231,16 @@ bool PView::readPOS(std::string fileName, int fileIndex)
return false;
}
- char str[256];
+ char str[256] = "XXX";
double version;
int format, size, index = -1;
while(1) {
- do {
- if(!fgets(str, 256, fp) || feof(fp))
+ while(str[0] != '$'){
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
- } while(str[0] != '$');
+ }
if(feof(fp))
break;
@@ -291,7 +291,7 @@ bool PView::readPOS(std::string fileName, int fileIndex)
}
do {
- if(!fgets(str, 256, fp) || feof(fp))
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
} while(str[0] != '$');
@@ -310,16 +310,16 @@ bool PView::readMSH(std::string fileName, int fileIndex)
return false;
}
- char str[256];
+ char str[256] = "XXX";
int index = -1;
bool binary = false, swap = false;
while(1) {
- do {
- if(!fgets(str, 256, fp) || feof(fp))
- break;
- } while(str[0] != '$');
+ while(str[0] != '$'){
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
+ break;
+ }
if(feof(fp))
break;
@@ -345,7 +345,7 @@ bool PView::readMSH(std::string fileName, int fileIndex)
if(fileIndex < 0 || fileIndex == index){
// read data info
if(!fgets(str, sizeof(str), fp)) return false;
- std::string name = extractDoubleQuotedString(str, 256);
+ std::string name = extractDoubleQuotedString(str, sizeof(str));
int timeStep, partition, interpolationScheme, numComp, numNodes;
double time;
if(!fgets(str, sizeof(str), fp)) return false;
@@ -373,10 +373,9 @@ bool PView::readMSH(std::string fileName, int fileIndex)
}
do {
- if(!fgets(str, 256, fp) || feof(fp))
+ if(!fgets(str, sizeof(str), fp) || feof(fp))
break;
} while(str[0] != '$');
-
}
fclose(fp);
diff --git a/Post/PViewData.cpp b/Post/PViewData.cpp
index c153ca4c240e6f9966e9559b67d1346b5a0ecd4a..0956aeabec3c50daf7115717e1f344712f256f86 100644
--- a/Post/PViewData.cpp
+++ b/Post/PViewData.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewData.cpp,v 1.12 2008-03-12 21:28:53 geuzaine Exp $
+// $Id: PViewData.cpp,v 1.13 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -23,6 +23,7 @@
//
#include "PViewData.h"
+#include "Numeric.h"
PViewData::PViewData()
: _dirty(true), _fileIndex(0)
@@ -31,7 +32,15 @@ PViewData::PViewData()
bool PViewData::empty()
{
- return (!getNumElements() &&
- !getNumStrings2D() &&
- !getNumStrings3D());
+ return (!getNumElements() && !getNumStrings2D() && !getNumStrings3D());
}
+
+void PViewData::getScalarValue(int ent, int ele, int nod, int step, double &val)
+{
+ int numComp = getNumComponents(ent, ele, step);
+ double d[9];
+ for(int comp = 0; comp < numComp; comp++)
+ getValue(ent, ele, nod, comp, step, d[comp]);
+ val = ComputeScalarRep(numComp, d);
+}
+
diff --git a/Post/PViewData.h b/Post/PViewData.h
index 53ba037c6d626b1f93d239bd737d29a57d4facd7..9172f4e0d0ae6af333aad111114c854c87193fd8 100644
--- a/Post/PViewData.h
+++ b/Post/PViewData.h
@@ -88,6 +88,10 @@ class PViewData {
// associated with the node-th node from the ele-th element in the
// ent-th entity
virtual void getValue(int ent, int ele, int nod, int comp, int step, double &val) = 0;
+ // Returns a scalar value (same as value for scalars, norm for
+ // vectors, etc.) associated with the node-th node from the ele-th
+ // element in the ent-th entity
+ virtual void getScalarValue(int ent, int ele, int nod, int step, double &val);
// Returns the number of edges of the ele-th element in the ent-th
// entity
virtual int getNumEdges(int ent, int ele) = 0;
@@ -122,12 +126,4 @@ class nameData{
std::vector<PViewData*> data;
};
-class nameDataLessThan{
- public:
- bool operator()(const nameData &n1, const nameData &n2) const
- {
- return n1.name < n2.name;
- }
-};
-
#endif
diff --git a/Post/PViewDataGModel.cpp b/Post/PViewDataGModel.cpp
index f75eb1c4712b8a264d15c9e5013993e7b222e04a..e947442c777b1f629c3ab2e135e00d34f3ef5dfe 100644
--- a/Post/PViewDataGModel.cpp
+++ b/Post/PViewDataGModel.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataGModel.cpp,v 1.27 2008-03-13 22:02:08 geuzaine Exp $
+// $Id: PViewDataGModel.cpp,v 1.28 2008-03-18 19:30:14 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -27,7 +27,6 @@
PViewDataGModel::PViewDataGModel(GModel *model)
: _model(model), _min(VAL_INF), _max(-VAL_INF)
-
{
// store vector of GEntities so we can index them efficiently
for(GModel::eiter it = _model->firstEdge(); it != _model->lastEdge(); ++it)
@@ -53,6 +52,10 @@ bool PViewDataGModel::finalize()
_min = std::min(_min, _nodeData[i]->getMin());
_max = std::max(_max, _nodeData[i]->getMax());
}
+ for(unsigned int i = 0; i < _elementData.size(); i++){
+ _min = std::min(_min, _elementData[i]->getMin());
+ _max = std::max(_max, _elementData[i]->getMax());
+ }
setDirty(false);
return true;
}
@@ -66,20 +69,28 @@ double PViewDataGModel::getTime(int step)
{
if(step < (int)_nodeData.size())
return _nodeData[step]->getTime();
+ else if(step < (int)_elementData.size())
+ return _elementData[step]->getTime();
return 0.;
}
double PViewDataGModel::getMin(int step)
{
if(step < 0) return _min;
- if(step < (int)_nodeData.size()) return _nodeData[step]->getMin();
+ if(step < (int)_nodeData.size())
+ return _nodeData[step]->getMin();
+ else if(step < (int)_elementData.size())
+ return _elementData[step]->getMin();
return 0.;
}
double PViewDataGModel::getMax(int step)
{
if(step < 0) return _max;
- if(step < (int)_nodeData.size()) return _nodeData[step]->getMax();
+ if(step < (int)_nodeData.size())
+ return _nodeData[step]->getMax();
+ else if(step < (int)_elementData.size())
+ return _elementData[step]->getMax();
return 0.;
}
@@ -114,16 +125,21 @@ void PViewDataGModel::getNode(int ent, int ele, int nod, double &x, double &y, d
int PViewDataGModel::getNumComponents(int ent, int ele, int step)
{
- // no range check here: we assume this is guarded by skipElement()
- return _nodeData[step]->getNumComp();
+ if(step < (int)_nodeData.size())
+ return _nodeData[step]->getNumComp();
+ else if(step < (int)_elementData.size())
+ return _elementData[step]->getNumComp();
+ return 1;
}
void PViewDataGModel::getValue(int ent, int ele, int nod, int comp, int step, double &val)
{
MVertex *v = _entities[ent]->getMeshElement(ele)->getVertex(nod);
int index = v->getDataIndex();
- // no range check here: we assume this is guarded by skipElement()
- val = _nodeData[step]->getData(index)[comp];
+ if(step < (int)_nodeData.size())
+ val = _nodeData[step]->getData(index)[comp];
+ //else if(step < (int)_elementData.size())
+ // val = _elementData[step]->getData(index)[nod * numComp + comp];
}
int PViewDataGModel::getNumEdges(int ent, int ele)
diff --git a/Post/PViewDataGModel.h b/Post/PViewDataGModel.h
index af4d2067719d676b26f6745b7b8106718aa8a302..e81fb2732ce30dbd53154626d2aea9be1126f28f 100644
--- a/Post/PViewDataGModel.h
+++ b/Post/PViewDataGModel.h
@@ -28,9 +28,9 @@
template<class real>
class stepData{
private:
- // the file the data was read from
+ // the file the data was read from (if empty, refer to PViewData)
std::string _fileName;
- // the index in the file
+ // the index in the file (if negative, refer to PViewData)
int _fileIndex;
// the value of the time step and value min/max
double _time, _min, _max;
@@ -80,8 +80,8 @@ class stepData{
for(unsigned int i = 0; i < _data->size(); i++)
if((*_data)[i]) delete [] (*_data)[i];
delete _data;
+ _data = 0;
}
- _data = 0;
}
};
diff --git a/benchmarks/3d/hexaprism.geo b/benchmarks/3d/hexaprism.geo
index eaf12fe1e2d153812ccb24c3658a2a624f92c6aa..807bbd5b1d0a391582e652615760a02723f5d743 100644
--- a/benchmarks/3d/hexaprism.geo
+++ b/benchmarks/3d/hexaprism.geo
@@ -8,11 +8,12 @@ Line(3) = {3,4};
Line(4) = {4,1};
Line Loop(5) = {2,3,4,1};
Plane Surface(6) = {5};
-Recombine Surface {6};
-Extrude Surface {6, {0.0,0.0,1.0}}{Layers{1,1,1};Recombine;};
+Transfinite Surface{6} = {1,2,3,4};
+Recombine Surface{6};
+Extrude Surface {6, {0.0,0.0,1.0}}{Layers{1};Recombine;};
Point(15) = {2,0.0,0.0,1.0};
Line(29) = {15,2};
Line(30) = {3,15};
Line Loop(31) = {30,29,2};
Plane Surface(32) = {31};
-Extrude Surface {32, {0.0,0.0,1.0}}{Layers{1,1,1};Recombine;};
+Extrude Surface {32, {0.0,0.0,1.0}}{Layers{1};Recombine;};
diff --git a/utils/misc/driver.cpp b/utils/misc/driver.cpp
index 5e6b766210db434b9fac5fa6919fb91fe17af2f0..bd8b2b1c1aa8ceb41eae8fd02e99f3f35270f298 100644
--- a/utils/misc/driver.cpp
+++ b/utils/misc/driver.cpp
@@ -17,9 +17,10 @@ int main(int argc, char **argv)
m->readGEO("../../tutorial/t5.geo");
m->mesh(3);
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it){
- printf("volume %d contains %d elements:\n", (*it)->tag(), (*it)->getNumMeshElements());
- for(unsigned int i = 0; i < (*it)->getNumMeshElements(); i++)
- printf(" %d", (*it)->getMeshElement(i)->getNum());
+ GRegion *r = *it;
+ printf("volume %d contains %d elements:\n", r->tag(), r->getNumMeshElements());
+ for(unsigned int i = 0; i < r->getNumMeshElements(); i++)
+ printf(" %d", r->getMeshElement(i)->getNum());
printf("\n");
}
m->writeMSH("test.msh");