diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp
index da8d096db907f58fcd11fa97c6e3a96d85256a5f..d53b3e0e251bbbc9f8fb650ab1ce1533824073a0 100644
--- a/Geo/GFace.cpp
+++ b/Geo/GFace.cpp
@@ -1,4 +1,4 @@
-// $Id: GFace.cpp,v 1.58 2008-02-23 15:40:29 geuzaine Exp $
+// $Id: GFace.cpp,v 1.59 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -497,13 +497,13 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
err = 1.0;
iter = 1;
- GPoint P = point(U,V);
+ GPoint P = point(U, V);
err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
if (err2 < 1.e-8 * CTX.lc) return;
while(err > Precision && iter < MaxIter) {
- P = point(U,V);
- Pair<SVector3,SVector3> der = firstDer(SPoint2(U,V));
+ P = point(U, V);
+ Pair<SVector3, SVector3> der = firstDer(SPoint2(U, V));
mat[0][0] = der.left().x();
mat[0][1] = der.left().y();
mat[0][2] = der.left().z();
@@ -522,6 +522,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
(jac[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
jac[2][1] * (Z - P.z()));
+ //if(Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin) break;
+
err = DSQR(Unew - U) + DSQR(Vnew - V);
err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
iter++;
@@ -529,11 +531,14 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
V = Vnew;
}
+ //printf("i=%d j=%d err=%g iter=%d err2=%g u=%.16g v=%.16g x=%g y=%g z=%g\n",
+ // i, j, err, iter, err2, U, V, X, Y, Z);
+
if(iter < MaxIter && err <= Precision &&
Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){
if (onSurface && err2 > 1.e-4 * CTX.lc)
- Msg(WARNING,"Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g",
+ Msg(WARNING, "Converged for i=%d j=%d (err=%g iter=%d) BUT xyz error = %g",
i, j, err, iter, err2);
return;
}
@@ -543,7 +548,7 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
if(relax < 1.e-6)
Msg(GERROR, "Could not converge: surface mesh will be wrong");
else {
- Msg(INFO, "point %g %g %g : Relaxation factor = %g",X,Y,Z, 0.75 * relax);
+ Msg(INFO, "point %g %g %g : Relaxation factor = %g", X, Y, Z, 0.75 * relax);
XYZtoUV(X, Y, Z, U, V, 0.75 * relax);
}
#endif
@@ -551,9 +556,9 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
SPoint2 GFace::parFromPoint(const SPoint3 &p) const
{
- double U,V;
- XYZtoUV(p.x(),p.y(),p.z(),U,V,1.0);
- return SPoint2(U,V);
+ double U, V;
+ XYZtoUV(p.x(), p.y(), p.z(), U, V, 1.0);
+ return SPoint2(U, V);
}
int GFace::containsParam(const SPoint2 &pt) const
diff --git a/Geo/GVertex.cpp b/Geo/GVertex.cpp
index 9b701cfb95ad8205d8fb179cc5360533d4ea3243..bab1210652406d81ff8ff95df9965c9538500576 100644
--- a/Geo/GVertex.cpp
+++ b/Geo/GVertex.cpp
@@ -1,4 +1,4 @@
-// $Id: GVertex.cpp,v 1.20 2008-02-22 07:59:00 geuzaine Exp $
+// $Id: GVertex.cpp,v 1.21 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -56,7 +56,7 @@ void GVertex::delEdge(GEdge *e)
l_edges.erase(std::find(l_edges.begin(), l_edges.end(), e));
}
-SPoint2 GVertex::reparamOnFace ( GFace *gf , int) const
+SPoint2 GVertex::reparamOnFace(GFace *gf, int) const
{
return gf->parFromPoint(SPoint3(x(), y(), z()));
}
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 89469dd539d9a3aee3ced30e06bd4ca0343333af..11cc7ddc97f13aec85200f90e6e71c728274e029 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -1,4 +1,4 @@
-// $Id: GeoInterpolation.cpp,v 1.32 2008-02-17 08:47:58 geuzaine Exp $
+// $Id: GeoInterpolation.cpp,v 1.33 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -27,7 +27,7 @@
// Curves
-Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
+Vertex InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
int derivee, double t1, double t2)
{
Vertex V;
@@ -55,7 +55,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[i] = 0.0;
}
- /* X */
+ // X
for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.X;
@@ -67,7 +67,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[j] = 0.0;
}
- /* Y */
+ // Y
for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.Y;
@@ -79,7 +79,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[j] = 0.0;
}
- /* Z */
+ // Z
for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.Z;
@@ -100,7 +100,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
}
// interpolation in the parametric space !
-SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
+SPoint2 InterpolateCubicSpline(Vertex *v[4], double t, double mat[4][4],
int derivee, double t1, double t2, gmshSurface *s)
{
Vertex V;
@@ -126,11 +126,10 @@ SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
p += coord[j] * T[j] ;
}
return p;
-
}
// Uniform BSplines
-Vertex InterpolateUBS(Curve * Curve, double u, int derivee)
+Vertex InterpolateUBS(Curve *Curve, double u, int derivee)
{
bool periodic = (Curve->end == Curve->beg);
int NbControlPoints = List_Nbr(Curve->Control_Points);
@@ -149,7 +148,7 @@ Vertex InterpolateUBS(Curve * Curve, double u, int derivee)
}
if(Curve->geometry){
- SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, 0, t1, t2,Curve->geometry);
+ SPoint2 pp = InterpolateCubicSpline(v, t, Curve->mat, 0, t1, t2, Curve->geometry);
SPoint3 pt = Curve->geometry->point(pp);
Vertex V;
V.Pos.X = pt.x();
@@ -204,7 +203,7 @@ void basisFuns(double u, int i, int deg, float *U, double *N)
}
}
-Vertex InterpolateNurbs(Curve * Curve, double u, int derivee)
+Vertex InterpolateNurbs(Curve *Curve, double u, int derivee)
{
static double Nb[1000];
int span = findSpan(u, Curve->degre, List_Nbr(Curve->Control_Points), Curve->k);
@@ -223,7 +222,7 @@ Vertex InterpolateNurbs(Curve * Curve, double u, int derivee)
return p;
}
-Vertex InterpolateCurve(Curve * c, double u, int derivee)
+Vertex InterpolateCurve(Curve *c, double u, int derivee)
{
if(c->Num < 0) {
Curve *C0 = FindCurve(-c->Num);
@@ -231,21 +230,21 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
Msg(GERROR, "Unknown curve %d", -c->Num);
return Vertex(0., 0., 0.);
}
- return InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1.-u), derivee);
+ return InterpolateCurve(C0, C0->ubeg + (C0->uend - C0->ubeg) * (1. - u), derivee);
}
Vertex V;
V.u = u;
if(derivee) {
- double eps1 = u==0?0:1.e-5;
- double eps2 = u==1?0:1.e-5;
+ double eps1 = (u == 0) ? 0 : 1.e-5;
+ double eps2 = (u == 1) ? 0 : 1.e-5;
Vertex D[2];
- D[0] = InterpolateCurve(c, u-eps1, 0);
- D[1] = InterpolateCurve(c, u+eps2, 0);
- V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1+eps2);
- V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1+eps2);
- V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1+eps2);
+ D[0] = InterpolateCurve(c, u - eps1, 0);
+ D[1] = InterpolateCurve(c, u + eps2, 0);
+ V.Pos.X = (D[1].Pos.X - D[0].Pos.X) / (eps1 + eps2);
+ V.Pos.Y = (D[1].Pos.Y - D[0].Pos.Y) / (eps1 + eps2);
+ V.Pos.Z = (D[1].Pos.Z - D[0].Pos.Z) / (eps1 + eps2);
return V;
}
@@ -394,9 +393,9 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
// Surfaces
-/* Interpolation transfinie sur un quadrangle :
- f(u,v) = (1-u)c4(v) + u c2(v) + (1-v)c1(u) + v c3(u)
- - [ (1-u)(1-v)s1 + u(1-v)s2 + uv s3 + (1-u)v s4 ] */
+// Interpolation transfinie sur un quadrangle :
+// f(u,v) = (1-u)c4(v) + u c2(v) + (1-v)c1(u) + v c3(u)
+// - [ (1-u)(1-v)s1 + u(1-v)s2 + uv s3 + (1-u)v s4 ]
#define TRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
(1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
@@ -419,8 +418,8 @@ Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
return (V);
}
-/* Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4
- f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3 */
+// Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4
+// f(u,v) = u c2 (v) + (1-v) c1(u) + v c3(u) - u(1-v) s2 - uv s3
#define TRAN_TRI(c1,c2,c3,s1,s2,s3,u,v) u*c2+(1.-v)*c1+v*c3-(u*(1.-v)*s2+u*v*s3);
@@ -428,7 +427,6 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v)
{
Vertex V;
-
V.lc = TRAN_TRI(c1.lc, c2.lc, c3.lc, s1.lc, s2.lc, s3.lc, u, v);
V.w = TRAN_TRI(c1.w, c2.w, c3.w, s1.w, s2.w, s3.w, u, v);
V.Pos.X = TRAN_TRI(c1.Pos.X, c2.Pos.X, c3.Pos.X,
@@ -437,53 +435,63 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
s1.Pos.Y, s2.Pos.Y, s3.Pos.Y, u, v);
V.Pos.Z = TRAN_TRI(c1.Pos.Z, c2.Pos.Z, c3.Pos.Z,
s1.Pos.Z, s2.Pos.Z, s3.Pos.Z, u, v);
- return (V);
+ return V;
}
-void TransfiniteSph(Vertex S, Vertex center, Vertex * T)
+void TransfiniteSph(Vertex S, Vertex center, Vertex *T)
{
- double r, s, dirx, diry, dirz;
-
- r = sqrt(DSQR(S.Pos.X - center.Pos.X) + DSQR(S.Pos.Y - center.Pos.Y)
- + DSQR(S.Pos.Z - center.Pos.Z));
+ double r = sqrt(DSQR(S.Pos.X - center.Pos.X) + DSQR(S.Pos.Y - center.Pos.Y)
+ + DSQR(S.Pos.Z - center.Pos.Z));
- s = sqrt(DSQR(T->Pos.X - center.Pos.X) + DSQR(T->Pos.Y - center.Pos.Y)
- + DSQR(T->Pos.Z - center.Pos.Z));
+ double s = sqrt(DSQR(T->Pos.X - center.Pos.X) + DSQR(T->Pos.Y - center.Pos.Y)
+ + DSQR(T->Pos.Z - center.Pos.Z));
- dirx = (T->Pos.X - center.Pos.X) / s;
- diry = (T->Pos.Y - center.Pos.Y) / s;
- dirz = (T->Pos.Z - center.Pos.Z) / s;
+ double dirx = (T->Pos.X - center.Pos.X) / s;
+ double diry = (T->Pos.Y - center.Pos.Y) / s;
+ double dirz = (T->Pos.Z - center.Pos.Z) / s;
T->Pos.X = center.Pos.X + r * dirx;
T->Pos.Y = center.Pos.Y + r * diry;
T->Pos.Z = center.Pos.Z + r * dirz;
}
-Vertex InterpolateRuledSurface(Surface * s, double u, double v)
+Vertex InterpolateRuledSurface(Surface *s, double u, double v)
{
- Curve *C[4];
- Vertex *c1, *c2;
- int issphere = 1;
- for(int i = 0; i < List_Nbr(s->Generatrices); i++) {
- if(i == 4) break; // we accept holes in ruled surfaces!
+ Curve *C[4] = {0, 0, 0, 0};
+ Vertex *O = 0;
+ bool isSphere = true;
+ for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++) {
List_Read(s->Generatrices, i, &C[i]);
- if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV) {
- issphere = 0;
+ if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV){
+ isSphere = false;
}
- else if(issphere) {
- if(!i) {
- List_Read(C[i]->Control_Points, 1, &c1);
+ else if(isSphere){
+ if(!i){
+ List_Read(C[i]->Control_Points, 1, &O);
}
- else {
- List_Read(C[i]->Control_Points, 1, &c2);
- if(compareVertex(&c1, &c2))
- issphere = 0;
+ else{
+ Vertex *tmp;
+ List_Read(C[i]->Control_Points, 1, &tmp);
+ if(compareVertex(&O, &tmp))
+ isSphere = false;
}
}
}
+ if(isSphere){
+ double n[3] = {C[0]->Circle.invmat[0][2],
+ C[0]->Circle.invmat[1][2],
+ C[0]->Circle.invmat[2][2]};
+ bool isPlane = true;
+ for(int i = 1; i < std::min(List_Nbr(s->Generatrices), 4); i++)
+ isPlane &= (n[0] == C[i]->Circle.invmat[0][2] &&
+ n[1] == C[i]->Circle.invmat[1][2] &&
+ n[2] == C[i]->Circle.invmat[2][2]);
+ if(isPlane)
+ isSphere = false;
+ }
+
Vertex *S[4], V[4], T;
-
if(s->Typ == MSH_SURF_REGL && List_Nbr(s->Generatrices) >= 4){
S[0] = C[0]->beg;
S[1] = C[1]->beg;
@@ -494,7 +502,7 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v)
V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
V[3] = InterpolateCurve(C[3], C[3]->ubeg + (C[3]->uend - C[3]->ubeg) * (1. - v), 0);
T = TransfiniteQua(V[0], V[1], V[2], V[3], *S[0], *S[1], *S[2], *S[3], u, v);
- if(issphere) TransfiniteSph(*S[0], *c1, &T);
+ if(isSphere) TransfiniteSph(*S[0], *O, &T);
}
else if(List_Nbr(s->Generatrices) >= 3){
S[0] = C[0]->beg;
@@ -502,15 +510,15 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v)
S[2] = C[2]->beg;
V[0] = InterpolateCurve(C[0], C[0]->ubeg + (C[0]->uend - C[0]->ubeg) * u, 0);
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 0);
- V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1.-u), 0);
+ V[2] = InterpolateCurve(C[2], C[2]->ubeg + (C[2]->uend - C[2]->ubeg) * (1. - u), 0);
T = TransfiniteTri(V[0], V[1], V[2], *S[0], *S[1], *S[2], u, v);
- if(issphere) TransfiniteSph(*S[0], *c1, &T);
+ if(isSphere) TransfiniteSph(*S[0], *O, &T);
}
return T;
}
-Vertex InterpolateExtrudedSurface(Surface * s, double u, double v)
+Vertex InterpolateExtrudedSurface(Surface *s, double u, double v)
{
Curve *c = FindCurve(s->Extrude->geo.Source);
@@ -548,9 +556,8 @@ Vertex InterpolateExtrudedSurface(Surface * s, double u, double v)
}
}
-Vertex InterpolateSurface(Surface * s, double u, double v, int derivee, int u_v)
+Vertex InterpolateSurface(Surface *s, double u, double v, int derivee, int u_v)
{
-
if(derivee) {
double eps = 1.e-6;
Vertex D[4], T;
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index e25d5ff2dc7a260eb65979141cae0e6972b88826..d6c72bfd7392357cb8360394425d8907767f526e 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.53 2008-02-23 17:38:34 geuzaine Exp $
+// $Id: gmshFace.cpp,v 1.54 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -175,8 +175,8 @@ Pair<SVector3,SVector3> gmshFace::firstDer(const SPoint2 ¶m) const
{
Vertex vu = InterpolateSurface(s, param[0], param[1], 1, 1);
Vertex vv = InterpolateSurface(s, param[0], param[1], 1, 2);
- return Pair<SVector3,SVector3>(SVector3(vu.Pos.X, vu.Pos.Y, vu.Pos.Z),
- SVector3(vv.Pos.X, vv.Pos.Y, vv.Pos.Z));
+ return Pair<SVector3, SVector3>(SVector3(vu.Pos.X, vu.Pos.Y, vu.Pos.Z),
+ SVector3(vv.Pos.X, vv.Pos.Y, vv.Pos.Z));
}
GPoint gmshFace::point(double par1, double par2) const
diff --git a/Post/PView.cpp b/Post/PView.cpp
index 2392cc74b43614d0c797081f2155dbc18ba054b6..d017befc217cbd3a21f917c4c4be67bfb59f1738 100644
--- a/Post/PView.cpp
+++ b/Post/PView.cpp
@@ -1,4 +1,4 @@
-// $Id: PView.cpp,v 1.18 2008-03-01 01:32:03 geuzaine Exp $
+// $Id: PView.cpp,v 1.19 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -211,6 +211,13 @@ void PView::combine(bool time, int how, bool remove)
delete *it;
}
+PView *PView::getViewByName(std::string name)
+{
+ for(unsigned int i = 0; i < list.size(); i++)
+ if(list[i]->getData()->getName() == name) return list[i];
+ return 0;
+}
+
bool PView::readPOS(std::string filename, int fileIndex)
{
FILE *fp = fopen(filename.c_str(), "rb");
diff --git a/Post/PView.h b/Post/PView.h
index 4e9b00dcc06731bd75bf9bde07aaaa77655b2f48..4813d122cc74703915293da0fd92e50d6729d62e 100644
--- a/Post/PView.h
+++ b/Post/PView.h
@@ -80,6 +80,9 @@ class PView{
// combine view
static void combine(bool time, int how, bool remove);
+ // combine view
+ static PView *getViewByName(std::string name);
+
// read view(s) in list format from a file
static bool readPOS(std::string filename, int fileIndex=-1);
// read view data from MSH file
diff --git a/Post/PViewDataListIO.cpp b/Post/PViewDataListIO.cpp
index 1f457e32aee6003ccfeb3574fa358baecb5c2e17..e7c1d7d57ce11c27adf8f2ed846dbb6e1fd732da 100644
--- a/Post/PViewDataListIO.cpp
+++ b/Post/PViewDataListIO.cpp
@@ -1,4 +1,4 @@
-// $Id: PViewDataListIO.cpp,v 1.11 2008-02-24 21:37:46 geuzaine Exp $
+// $Id: PViewDataListIO.cpp,v 1.12 2008-03-03 22:04:22 geuzaine Exp $
//
// Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
//
@@ -425,7 +425,7 @@ static void getNodeMSH(int nbelm, List_T *list, int nbnod, int nbcomp,
{
if(!nbelm) return;
int nb = List_Nbr(list) / nbelm;
- for(int i = 0; i < List_Nbr(list); i+=nb){
+ for(int i = 0; i < List_Nbr(list); i += nb){
double *x = (double *)List_Pointer_Fast(list, i);
double *y = (double *)List_Pointer_Fast(list, i + nbnod);
double *z = (double *)List_Pointer_Fast(list, i + 2 * nbnod);
@@ -469,33 +469,33 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
switch(dim){
case 0:
- fprintf(fp, "%d 15 %d %d 1 %d\n", num, phys, ele, nod[0].Num);
+ fprintf(fp, "%d 15 2 %d %d %d\n", num, phys, ele, nod[0].Num);
break;
case 1:
- fprintf(fp, "%d 1 %d %d 2 %d %d\n", num, phys, ele, nod[0].Num, nod[1].Num);
+ fprintf(fp, "%d 1 2 %d %d %d %d\n", num, phys, ele, nod[0].Num, nod[1].Num);
break;
case 2:
if(nbnod == 3)
- fprintf(fp, "%d 2 %d %d 3 %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 2 2 %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num);
else
- fprintf(fp, "%d 3 %d %d 4 %d %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 3 2 %d %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
break;
case 3:
default:
if(nbnod == 4)
- fprintf(fp, "%d 4 %d %d 4 %d %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 4 2 %d %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
else if(nbnod == 5)
- fprintf(fp, "%d 7 %d %d 5 %d %d %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 7 2 %d %d %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num);
else if(nbnod == 6)
- fprintf(fp, "%d 6 %d %d 6 %d %d %d %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 6 2 %d %d %d %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
nod[5].Num);
else
- fprintf(fp, "%d 5 %d %d 8 %d %d %d %d %d %d %d %d\n", num, phys, ele,
+ fprintf(fp, "%d 5 2 %d %d %d %d %d %d %d %d %d %d\n", num, phys, ele,
nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
nod[5].Num, nod[6].Num, nod[7].Num);
break;
@@ -510,7 +510,7 @@ static void writeElementsMSH(FILE *fp, int nbelm, List_T *list,
if(!nbelm) return;
pVertex nod[8];
int nb = List_Nbr(list) / nbelm;
- for(int i = 0; i < List_Nbr(list); i+=nb){
+ for(int i = 0; i < List_Nbr(list); i += nb){
double *x = (double *)List_Pointer_Fast(list, i);
double *y = (double *)List_Pointer_Fast(list, i + nbnod);
double *z = (double *)List_Pointer_Fast(list, i + 2 * nbnod);
@@ -566,15 +566,16 @@ bool PViewDataList::writeMSH(std::string name)
getNodeMSH(NbVY, VY, 5, 3, &nodes, &numelm);
getNodeMSH(NbTY, TY, 5, 9, &nodes, &numelm);
- fprintf(fp, "$NOD\n");
+ fprintf(fp, "$MeshFormat\n2 0 8\n$EndMeshFormat\n");
+ fprintf(fp, "$Nodes\n");
fprintf(fp, "%d\n", (int)nodes.size());
for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin();
it != nodes.end(); ++it){
fprintf(fp, "%d %.16g %.16g %.16g\n", it->Num, it->X, it->Y, it->Z);
}
- fprintf(fp, "$ENDNOD\n");
+ fprintf(fp, "$EndNodes\n");
- fprintf(fp, "$ELM\n");
+ fprintf(fp, "$Elements\n");
fprintf(fp, "%d\n", numelm);
numelm = 0;
writeElementsMSH(fp, NbSP, SP, 1, 1, 0, &nodes, &numelm);
@@ -601,19 +602,25 @@ bool PViewDataList::writeMSH(std::string name)
writeElementsMSH(fp, NbSY, SY, 5, 1, 3, &nodes, &numelm);
writeElementsMSH(fp, NbVY, VY, 5, 3, 3, &nodes, &numelm);
writeElementsMSH(fp, NbTY, TY, 5, 9, 3, &nodes, &numelm);
- fprintf(fp, "$ENDELM\n");
+ fprintf(fp, "$EndElements\n");
-#if 0 // test new postpro node-based storage
- fprintf(fp, "$NodeData\n");
- fprintf(fp, "\"%s\"\n", getName().c_str());
- fprintf(fp, "1 1 %d\n", nodes.size());
- for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin();
- it != nodes.end(); ++it){
- fprintf(fp, "%d", it->Num);
- for(int i = 0; i < it->Val.size(); i++) fprintf(fp, " %d", it->Val[i]);
- fprintf(fp, "\n");
+#if 1 // test new postpro node-based storage
+ int numNodes = nodes.size();
+ if(numNodes){
+ fprintf(fp, "$NodeData\n");
+ fprintf(fp, "\"%s\"\n", getName().c_str());
+ int timeStep = 1, numComp = nodes.begin()->Val.size();
+ double time = 0.;
+ fprintf(fp, "%d %.16g %d %d\n", timeStep, time, numComp, numNodes);
+ for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin();
+ it != nodes.end(); ++it){
+ fprintf(fp, "%d", it->Num);
+ for(int i = 0; i < it->Val.size(); i++)
+ fprintf(fp, " %.16g", it->Val[i]);
+ fprintf(fp, "\n");
+ }
+ fprintf(fp, "$EndNodeData\n");
}
- fprintf(fp, "$EndNodeData\n");
#endif
fclose(fp);
diff --git a/benchmarks/2d/Square-Attr1.geo b/benchmarks/2d/Square-Attr1.geo
index 362b087f5d0170c98bebba5761bac65cfd263385..1beae2053366da2e396461c40c4cf96f6a0516d4 100644
--- a/benchmarks/2d/Square-Attr1.geo
+++ b/benchmarks/2d/Square-Attr1.geo
@@ -17,5 +17,5 @@ Point(11) = {0 ,0.,0,lc};
Point(22) = {1.,1.,0,lc};
Line(5) = {11,22};
//Attractor Point{1,2,3,4,11} = {.0001,.0001,17};
-Attractor Line{3} = {.1,lc2/30,lc2,1000,3};
+Attractor Line{3,5} = {.1,lc2/30,lc2,1000,3};
diff --git a/benchmarks/2d/Square-Attr4.geo b/benchmarks/2d/Square-Attr4.geo
index 8c5cf71d4c8eeb9bbd42b3248572668c5eb4546c..bb8fa75a6f40032225b49ff08306e2be18b2c476 100644
--- a/benchmarks/2d/Square-Attr4.geo
+++ b/benchmarks/2d/Square-Attr4.geo
@@ -1,18 +1,20 @@
-/******************************
-Square uniformly meshed
-******************************/
-lc = .149999;
-Point(1) = {0.0,0.0,0,lc};
-Point(2) = {1,0.0,0,lc};
-Point(3) = {1,1,0,lc};
-Point(4) = {0,1,0,lc};
-Line(1) = {3,2};
-Line(2) = {2,1};
-Line(3) = {1,4};
-Line(4) = {4,3};
-Point(55) = {0.2,.5,0,lc};
-Line Loop(5) = {1,2,3,4};
-Plane Surface(6) = {5};
-Attractor Point {55} = {.01,.1,3.0};
-Mesh.Algorithm = 2;
+lc = .15;
+Point(1) = {0.0,0.0,0,lc};
+Point(2) = {1,0.0,0,lc};
+Point(3) = {1,1,0,lc};
+Point(4) = {0,1,0,lc};
+Line(1) = {3,2};
+Line(2) = {2,1};
+Line(3) = {1,4};
+Line(4) = {4,3};
+Point(55) = {0.2,.5,0,lc};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+num_pts = 100; // number of points on the attractor, unused for points
+lc_min = lc/50; // lc inside r_min
+lc_max = lc; // lc outside r_max
+r_min = 0.05;
+r_max = 0.4;
+Attractor Point{1,55} = {r_min, lc_min, lc_max, num_pts, r_max / r_min};
+Attractor Line{1} = {r_min, lc_min, lc_max, num_pts, r_max / r_min};
diff --git a/demos/attractors.geo b/demos/attractors.geo
new file mode 100644
index 0000000000000000000000000000000000000000..78d148a19346d96ceeb4ce86588500c6660988f5
--- /dev/null
+++ b/demos/attractors.geo
@@ -0,0 +1,28 @@
+lc = .15;
+Point(1) = {0.0,0.0,0,lc};
+Point(2) = {1,0.0,0,lc};
+Point(3) = {1,1,0,lc};
+Point(4) = {0,1,0,lc};
+Line(1) = {3,2};
+Line(2) = {2,1};
+Line(3) = {1,4};
+Line(4) = {4,3};
+Point(55) = {0.2,.5,0,lc};
+Line Loop(5) = {1,2,3,4};
+Plane Surface(6) = {5};
+
+// Point and line attractors (shortcuts for Threshold fields,
+// automatically added to the list of char length fields)
+
+num_pts = 100; // number of points on the attractor, unused for points
+lc_min = lc/20; // lc inside r_min
+lc_max = lc; // lc outside r_max
+r_min = 0.15;
+r_max = 0.5;
+Attractor Point{1,55} = {r_min, lc_min, lc_max, num_pts, r_max / r_min};
+Attractor Line{1} = {r_min, lc_min, lc_max, num_pts, r_max / r_min};
+
+// Function field
+
+Function Field(1) = "Cos(4*3.14*x) * Sin(4*3.14*y) / 10 + 0.101";
+Characteristic Length Field{1};