Skip to content
Snippets Groups Projects
Commit 7a3ca986 authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

- fixed isSphere test when all circle arcs ar in the same plane (should never 

do this, but some people actually do...)
- added some attractor examples
parent dbfc224e
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/GFace.cpp
+ 14
− 9
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -522,6 +522,8 @@ void GFace::XYZtoUV(const double X, const double Y, const double 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[0][1] * (X - P.x()) + jac[1][1] * (Y - P.y()) +
jac[2][1] * (Z - P.z())); jac[2][1] * (Z - P.z()));
//if(Unew > umax || Vnew > vmax || Unew < umin || Vnew < vmin) break;
err = DSQR(Unew - U) + DSQR(Vnew - V); err = DSQR(Unew - U) + DSQR(Vnew - V);
err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z())); err2 = sqrt(DSQR(X - P.x()) + DSQR(Y - P.y()) + DSQR(Z - P.z()));
iter++; iter++;
...@@ -529,6 +531,9 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z, ...@@ -529,6 +531,9 @@ void GFace::XYZtoUV(const double X, const double Y, const double Z,
V = Vnew; 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 && if(iter < MaxIter && err <= Precision &&
Unew <= umax && Vnew <= vmax && Unew <= umax && Vnew <= vmax &&
Unew >= umin && Vnew >= vmin){ Unew >= umin && Vnew >= vmin){
......
Geo/GVertex.cpp
+ 2
− 2
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
......
Geo/GeoInterpolation.cpp
+ 65
− 58
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -55,7 +55,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4], ...@@ -55,7 +55,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[i] = 0.0; vec[i] = 0.0;
} }
/* X */ // X
for(i = 0; i < 4; i++) { for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) { for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.X; vec[i] += mat[i][j] * v[j]->Pos.X;
...@@ -67,7 +67,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4], ...@@ -67,7 +67,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[j] = 0.0; vec[j] = 0.0;
} }
/* Y */ // Y
for(i = 0; i < 4; i++) { for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) { for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.Y; vec[i] += mat[i][j] * v[j]->Pos.Y;
...@@ -79,7 +79,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4], ...@@ -79,7 +79,7 @@ Vertex InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
vec[j] = 0.0; vec[j] = 0.0;
} }
/* Z */ // Z
for(i = 0; i < 4; i++) { for(i = 0; i < 4; i++) {
for(j = 0; j < 4; j++) { for(j = 0; j < 4; j++) {
vec[i] += mat[i][j] * v[j]->Pos.Z; vec[i] += mat[i][j] * v[j]->Pos.Z;
...@@ -126,7 +126,6 @@ SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4], ...@@ -126,7 +126,6 @@ SPoint2 InterpolateCubicSpline(Vertex * v[4], double t, double mat[4][4],
p += coord[j] * T[j] ; p += coord[j] * T[j] ;
} }
return p; return p;
} }
// Uniform BSplines // Uniform BSplines
...@@ -238,8 +237,8 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee) ...@@ -238,8 +237,8 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
V.u = u; V.u = u;
if(derivee) { if(derivee) {
double eps1 = u==0?0:1.e-5; double eps1 = (u == 0) ? 0 : 1.e-5;
double eps2 = u==1?0:1.e-5; double eps2 = (u == 1) ? 0 : 1.e-5;
Vertex D[2]; Vertex D[2];
D[0] = InterpolateCurve(c, u - eps1, 0); D[0] = InterpolateCurve(c, u - eps1, 0);
D[1] = InterpolateCurve(c, u + eps2, 0); D[1] = InterpolateCurve(c, u + eps2, 0);
...@@ -394,9 +393,9 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee) ...@@ -394,9 +393,9 @@ Vertex InterpolateCurve(Curve * c, double u, int derivee)
// Surfaces // Surfaces
/* Interpolation transfinie sur un quadrangle : // Interpolation transfinie sur un quadrangle :
f(u,v) = (1-u)c4(v) + u c2(v) + (1-v)c1(u) + v c3(u) // 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 ] */ // - [ (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) \ #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) (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, ...@@ -419,8 +418,8 @@ Vertex TransfiniteQua(Vertex c1, Vertex c2, Vertex c3, Vertex c4,
return (V); return (V);
} }
/* Interpolation transfinie sur un triangle : TRAN_QUA avec s1=s4=c4 // 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 */ // 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); #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, ...@@ -428,7 +427,6 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
Vertex s1, Vertex s2, Vertex s3, double u, double v) Vertex s1, Vertex s2, Vertex s3, double u, double v)
{ {
Vertex V; Vertex V;
V.lc = TRAN_TRI(c1.lc, c2.lc, c3.lc, s1.lc, s2.lc, s3.lc, u, 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.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, V.Pos.X = TRAN_TRI(c1.Pos.X, c2.Pos.X, c3.Pos.X,
...@@ -437,22 +435,20 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3, ...@@ -437,22 +435,20 @@ Vertex TransfiniteTri(Vertex c1, Vertex c2, Vertex c3,
s1.Pos.Y, s2.Pos.Y, s3.Pos.Y, u, v); 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, 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); 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; double r = sqrt(DSQR(S.Pos.X - center.Pos.X) + DSQR(S.Pos.Y - center.Pos.Y)
r = sqrt(DSQR(S.Pos.X - center.Pos.X) + DSQR(S.Pos.Y - center.Pos.Y)
+ DSQR(S.Pos.Z - center.Pos.Z)); + DSQR(S.Pos.Z - center.Pos.Z));
s = sqrt(DSQR(T->Pos.X - center.Pos.X) + DSQR(T->Pos.Y - center.Pos.Y) double s = sqrt(DSQR(T->Pos.X - center.Pos.X) + DSQR(T->Pos.Y - center.Pos.Y)
+ DSQR(T->Pos.Z - center.Pos.Z)); + DSQR(T->Pos.Z - center.Pos.Z));
dirx = (T->Pos.X - center.Pos.X) / s; double dirx = (T->Pos.X - center.Pos.X) / s;
diry = (T->Pos.Y - center.Pos.Y) / s; double diry = (T->Pos.Y - center.Pos.Y) / s;
dirz = (T->Pos.Z - center.Pos.Z) / s; double dirz = (T->Pos.Z - center.Pos.Z) / s;
T->Pos.X = center.Pos.X + r * dirx; T->Pos.X = center.Pos.X + r * dirx;
T->Pos.Y = center.Pos.Y + r * diry; T->Pos.Y = center.Pos.Y + r * diry;
...@@ -461,29 +457,41 @@ void TransfiniteSph(Vertex S, Vertex center, Vertex * T) ...@@ -461,29 +457,41 @@ void TransfiniteSph(Vertex S, Vertex center, Vertex * T)
Vertex InterpolateRuledSurface(Surface *s, double u, double v) Vertex InterpolateRuledSurface(Surface *s, double u, double v)
{ {
Curve *C[4]; Curve *C[4] = {0, 0, 0, 0};
Vertex *c1, *c2; Vertex *O = 0;
int issphere = 1; bool isSphere = true;
for(int i = 0; i < List_Nbr(s->Generatrices); i++) { for(int i = 0; i < std::min(List_Nbr(s->Generatrices), 4); i++) {
if(i == 4) break; // we accept holes in ruled surfaces!
List_Read(s->Generatrices, i, &C[i]); List_Read(s->Generatrices, i, &C[i]);
if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV){ if(C[i]->Typ != MSH_SEGM_CIRC && C[i]->Typ != MSH_SEGM_CIRC_INV){
issphere = 0; isSphere = false;
} }
else if(issphere) { else if(isSphere){
if(!i){ if(!i){
List_Read(C[i]->Control_Points, 1, &c1); List_Read(C[i]->Control_Points, 1, &O);
} }
else{ else{
List_Read(C[i]->Control_Points, 1, &c2); Vertex *tmp;
if(compareVertex(&c1, &c2)) List_Read(C[i]->Control_Points, 1, &tmp);
issphere = 0; if(compareVertex(&O, &tmp))
isSphere = false;
} }
} }
} }
Vertex *S[4], V[4], T; 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){ if(s->Typ == MSH_SURF_REGL && List_Nbr(s->Generatrices) >= 4){
S[0] = C[0]->beg; S[0] = C[0]->beg;
S[1] = C[1]->beg; S[1] = C[1]->beg;
...@@ -494,7 +502,7 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v) ...@@ -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[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); 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); 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){ else if(List_Nbr(s->Generatrices) >= 3){
S[0] = C[0]->beg; S[0] = C[0]->beg;
...@@ -504,7 +512,7 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v) ...@@ -504,7 +512,7 @@ Vertex InterpolateRuledSurface(Surface * s, double u, double v)
V[1] = InterpolateCurve(C[1], C[1]->ubeg + (C[1]->uend - C[1]->ubeg) * v, 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); 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; return T;
...@@ -550,7 +558,6 @@ Vertex InterpolateExtrudedSurface(Surface * s, double u, double v) ...@@ -550,7 +558,6 @@ 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) { if(derivee) {
double eps = 1.e-6; double eps = 1.e-6;
Vertex D[4], T; Vertex D[4], T;
......
Geo/gmshFace.cpp
+ 3
− 3
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
......
Post/PView.cpp
+ 8
− 1
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -211,6 +211,13 @@ void PView::combine(bool time, int how, bool remove) ...@@ -211,6 +211,13 @@ void PView::combine(bool time, int how, bool remove)
delete *it; 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) bool PView::readPOS(std::string filename, int fileIndex)
{ {
FILE *fp = fopen(filename.c_str(), "rb"); FILE *fp = fopen(filename.c_str(), "rb");
......
Post/PView.h
+ 3
− 0
View file @ 7a3ca986
...@@ -80,6 +80,9 @@ class PView{ ...@@ -80,6 +80,9 @@ class PView{
// combine view // combine view
static void combine(bool time, int how, bool remove); 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 // read view(s) in list format from a file
static bool readPOS(std::string filename, int fileIndex=-1); static bool readPOS(std::string filename, int fileIndex=-1);
// read view data from MSH file // read view data from MSH file
......
Post/PViewDataListIO.cpp
+ 32
− 25
View file @ 7a3ca986
// $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 // Copyright (C) 1997-2008 C. Geuzaine, J.-F. Remacle
// //
...@@ -469,33 +469,33 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8], ...@@ -469,33 +469,33 @@ static void writeElementMSH(FILE *fp, int num, int nbnod, pVertex nod[8],
switch(dim){ switch(dim){
case 0: 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; break;
case 1: 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; break;
case 2: case 2:
if(nbnod == 3) 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); nod[0].Num, nod[1].Num, nod[2].Num);
else 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); nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
break; break;
case 3: case 3:
default: default:
if(nbnod == 4) 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); nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num);
else if(nbnod == 5) 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); nod[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num);
else if(nbnod == 6) 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[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
nod[5].Num); nod[5].Num);
else 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[0].Num, nod[1].Num, nod[2].Num, nod[3].Num, nod[4].Num,
nod[5].Num, nod[6].Num, nod[7].Num); nod[5].Num, nod[6].Num, nod[7].Num);
break; break;
...@@ -566,15 +566,16 @@ bool PViewDataList::writeMSH(std::string name) ...@@ -566,15 +566,16 @@ bool PViewDataList::writeMSH(std::string name)
getNodeMSH(NbVY, VY, 5, 3, &nodes, &numelm); getNodeMSH(NbVY, VY, 5, 3, &nodes, &numelm);
getNodeMSH(NbTY, TY, 5, 9, &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()); fprintf(fp, "%d\n", (int)nodes.size());
for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin(); for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin();
it != nodes.end(); ++it){ it != nodes.end(); ++it){
fprintf(fp, "%d %.16g %.16g %.16g\n", it->Num, it->X, it->Y, it->Z); 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); fprintf(fp, "%d\n", numelm);
numelm = 0; numelm = 0;
writeElementsMSH(fp, NbSP, SP, 1, 1, 0, &nodes, &numelm); writeElementsMSH(fp, NbSP, SP, 1, 1, 0, &nodes, &numelm);
...@@ -601,19 +602,25 @@ bool PViewDataList::writeMSH(std::string name) ...@@ -601,19 +602,25 @@ bool PViewDataList::writeMSH(std::string name)
writeElementsMSH(fp, NbSY, SY, 5, 1, 3, &nodes, &numelm); writeElementsMSH(fp, NbSY, SY, 5, 1, 3, &nodes, &numelm);
writeElementsMSH(fp, NbVY, VY, 5, 3, 3, &nodes, &numelm); writeElementsMSH(fp, NbVY, VY, 5, 3, 3, &nodes, &numelm);
writeElementsMSH(fp, NbTY, TY, 5, 9, 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 #if 1 // test new postpro node-based storage
int numNodes = nodes.size();
if(numNodes){
fprintf(fp, "$NodeData\n"); fprintf(fp, "$NodeData\n");
fprintf(fp, "\"%s\"\n", getName().c_str()); fprintf(fp, "\"%s\"\n", getName().c_str());
fprintf(fp, "1 1 %d\n", nodes.size()); 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(); for(std::set<pVertex, pVertexLessThan>::iterator it = nodes.begin();
it != nodes.end(); ++it){ it != nodes.end(); ++it){
fprintf(fp, "%d", it->Num); fprintf(fp, "%d", it->Num);
for(int i = 0; i < it->Val.size(); i++) fprintf(fp, " %d", it->Val[i]); for(int i = 0; i < it->Val.size(); i++)
fprintf(fp, " %.16g", it->Val[i]);
fprintf(fp, "\n"); fprintf(fp, "\n");
} }
fprintf(fp, "$EndNodeData\n"); fprintf(fp, "$EndNodeData\n");
}
#endif #endif
fclose(fp); fclose(fp);
......
benchmarks/2d/Square-Attr1.geo
+ 1
− 1
View file @ 7a3ca986
...@@ -17,5 +17,5 @@ Point(11) = {0 ,0.,0,lc}; ...@@ -17,5 +17,5 @@ Point(11) = {0 ,0.,0,lc};
Point(22) = {1.,1.,0,lc}; Point(22) = {1.,1.,0,lc};
Line(5) = {11,22}; Line(5) = {11,22};
//Attractor Point{1,2,3,4,11} = {.0001,.0001,17}; //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};
benchmarks/2d/Square-Attr4.geo
+ 19
− 17
View file @ 7a3ca986
/****************************** lc = .15;
Square uniformly meshed
******************************/
lc = .149999;
Point(1) = {0.0,0.0,0,lc}; Point(1) = {0.0,0.0,0,lc};
Point(2) = {1,0.0,0,lc}; Point(2) = {1,0.0,0,lc};
Point(3) = {1,1,0,lc}; Point(3) = {1,1,0,lc};
...@@ -13,6 +10,11 @@ Line(4) = {4,3}; ...@@ -13,6 +10,11 @@ Line(4) = {4,3};
Point(55) = {0.2,.5,0,lc}; Point(55) = {0.2,.5,0,lc};
Line Loop(5) = {1,2,3,4}; Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5}; Plane Surface(6) = {5};
Attractor Point {55} = {.01,.1,3.0};
Mesh.Algorithm = 2;
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};
demos/attractors.geo 0 → 100644
+ 28
− 0
View file @ 7a3ca986
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};
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment