Spencer,
Rather than try and pull this off my machine from the aeronautics dept, I
thought I'd just send you the patch.
This is the full, unadulterated patch to the source as I have it here. It
includes all the trial and error changes to the global variables as well as
the genuine leak fixes.
Angus
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Basis.C clean2/Hlib/src/Basis.C
--- dirty/Hlib/src/Basis.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Basis.C Wed May 1 17:51:41 2002
@@ -603,142 +603,205 @@ static double norm(int np, double *d, ch
/* interpolation matrix */
-static double ***BasisA_G, ***BasisA_GL, ***BasisA_GR;
-static double ****BasisB_G, ****BasisB_GR;
-static double ****BasisC_GR, ****BasisC_G;
-
-/* This function initialises the orthogonal basis statically defined
- matrices based on QGmax. These arrays may be accessed form [1-QGmax]*/
-
-void init_ortho_basis(void){
- BasisA_G = (double ***) calloc(QGmax, sizeof(double **));
- BasisA_G -= 1;
- BasisA_GL = (double ***) calloc(QGmax, sizeof(double **));
- BasisA_GL-= 1;
- BasisA_GR = (double ***) calloc(QGmax, sizeof(double **));
- BasisA_GR-= 1;
+// This class initialises the orthogonal basis matrices based on QGmax.
+// These arrays may be accessed from [1-QGmax]
+class OBMatrices {
+public:
+ // This is a singleton class. Get the instance.
+ static OBMatrices & get();
+
+ // Some compilers don't like the d-tor being private.
+ ~OBMatrices();
+
+ double *** A_G;
+ double *** A_GL;
+ double *** A_GR;
+ double **** B_G;
+ double **** B_GR;
+ double **** C_GR;
+ double **** C_G;
+
+private:
+ // Make the c-tor private so we can control how many objects
+ // are instantiated.
+ OBMatrices();
+};
+
+OBMatrices & OBMatrices::get()
+{
+ static OBMatrices singleton;
+ return singleton;
+}
+
+OBMatrices::OBMatrices()
+{
+ A_G = (double ***) calloc(QGmax, sizeof(double **));
+ A_G -= 1;
+ A_GL = (double ***) calloc(QGmax, sizeof(double **));
+ A_GL-= 1;
+ A_GR = (double ***) calloc(QGmax, sizeof(double **));
+ A_GR-= 1;
- BasisB_G = (double ****)calloc(QGmax, sizeof(double ***));
- BasisB_G -= 1;
- BasisB_GR = (double ****)calloc(QGmax, sizeof(double ***));
- BasisB_GR-= 1;
+ B_G = (double ****)calloc(QGmax, sizeof(double ***));
+ B_G -= 1;
+ B_GR = (double ****)calloc(QGmax, sizeof(double ***));
+ B_GR-= 1;
- BasisC_GR = (double ****)calloc(QGmax, sizeof(double ***));
- BasisC_GR-= 1;
+ C_GR = (double ****)calloc(QGmax, sizeof(double ***));
+ C_GR-= 1;
- BasisC_G = (double ****)calloc(QGmax, sizeof(double ***));
- BasisC_G-= 1;
+ C_G = (double ****)calloc(QGmax, sizeof(double ***));
+ C_G-= 1;
+}
+OBMatrices::~OBMatrices()
+{
+ // i starts at 1 because of the offset, above.
+ for (int i = 1; i <= QGmax; ++i) {
+ if (A_G[i]) {
+ free_dmatrix(A_G[i], 0, 0);
+ }
+ if (A_GL[i]) {
+ free_dmatrix(A_GL[i], 0, 0);
+ }
+ if (A_GR[i]) {
+ free_dmatrix(A_GR[i], 0, 0);
+ }
+ if (B_G[i]) {
+ free(B_G[i][0][0]); free(B_G[i][0]); free(B_G[i]);
+ }
+ if (B_GR[i]) {
+ free(B_GR[i][0][0]); free(B_GR[i][0]); free(B_GR[i]);
+ }
+ if (C_G[i]) {
+ free(C_G[i][0][0]); free(C_G[i][0]); free(C_G[i]);
+ }
+ if (C_GR[i]) {
+ free(C_GR[i][0][0]); free(C_GR[i][0]); free(C_GR[i]);
+ }
+ }
+ // The increment is needed to return to the address of the allocation.
+ free(++A_G);
+ free(++A_GL);
+ free(++A_GR);
+ free(++B_G);
+ free(++B_GR);
+ free(++C_GR);
+ free(++C_G);
}
+// Now that we use a singleton class to store the Matrices, this function is
+// redundant.
+void init_ortho_basis()
+{}
+
/* Gets the values of "mode" (0..LGmax) degree Leg. poly on "size"
- Gauss points, stores the result in BasisA_G[size][mode] array */
-void get_moda_G (int size, double ***mat){
+ Gauss points, stores the result in OBMatrices::get().A_G[size][mode] array */void get_moda_G (int size, double ***mat){
- if (!BasisA_G[size]){
+ if (!OBMatrices::get().A_G[size]){
double *w, *z;
int i;
- BasisA_G[size] = dmatrix(0,LGmax-1,0,size-1);
+ OBMatrices::get().A_G[size] = dmatrix(0,LGmax-1,0,size-1);
getzw(size,&z,&w,'g');
for (i=0; i < LGmax; ++i){
- jacobf(size, z, BasisA_G[size][i], i, 0.0, 0.0);
+ jacobf(size, z, OBMatrices::get().A_G[size][i], i, 0.0, 0.0);
/* normalise so that (g,g)=1*/
#ifndef NONORM
- dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_G[size][i],1);
+ dscal(size,sqrt(0.5*(2.0*i+1.0)),OBMatrices::get().A_G[size][i],1);
#endif
}
}
- *mat = BasisA_G[size];
+ *mat = OBMatrices::get().A_G[size];
return;
}
/* computes the values of a_mode[mode_a][j], mode_a=0..LGmax,
j=0..size-1, at size Gauss-Lob. points, stores the result in
- BasisA_GL[size][mode_a] array */
+ OBMatrices::get().A_GL[size][mode_a] array */
void get_moda_GL (int size, double ***mat){
- if (!BasisA_GL[size]){
+ if (!OBMatrices::get().A_GL[size]){
double *w, *z;
int i;
double *tmp = dvector(0, size-1);
- BasisA_GL[size] = dmatrix(0,LGmax-1,0,size-1);
+ OBMatrices::get().A_GL[size] = dmatrix(0,LGmax-1,0,size-1);
getzw(size,&z,&w,'a');
for (i=0; i < LGmax; ++i){
- jacobf(size, z, BasisA_GL[size][i], i, 0.0, 0.0);
+ jacobf(size, z, OBMatrices::get().A_GL[size][i], i, 0.0, 0.0);
#ifndef NONORM
/* normalise so that (g,g)=1*/
- dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_GL[size][i],1);
+ dscal(size,sqrt(0.5*(2.0*i+1.0)),OBMatrices::get().A_GL[size][i],1);
#endif
}
free(tmp);
}
- *mat = BasisA_GL[size];
+ *mat = OBMatrices::get().A_GL[size];
return;
}
/* computes the values of a_mode[mode_a][j], mode_a=0..LGmax,
j=0..size-1, at size Gauss-Lob. points, stores the result in
- BasisA_GL[size][mode_a] array */
+ OBMatrices::get().A_GL[size][mode_a] array */
void get_moda_GR (int size, double ***mat){
- if (!BasisA_GR[size]){
+ if (!OBMatrices::get().A_GR[size]){
double *w, *z;
int i;
double *tmp = dvector(0, size-1);
- BasisA_GR[size] = dmatrix(0,LGmax-1,0,size-1);
+ OBMatrices::get().A_GR[size] = dmatrix(0,LGmax-1,0,size-1);
getzw(size,&z,&w,'b');
for (i=0; i < LGmax; ++i){
- jacobf(size, z, BasisA_GR[size][i], i, 0.0, 0.0);
+ jacobf(size, z, OBMatrices::get().A_GR[size][i], i, 0.0, 0.0);
#ifndef NONORM
/* normalise so that (g,g)=1*/
- dscal(size,sqrt(0.5*(2.0*i+1.0)),BasisA_GR[size][i],1);
+ dscal(size,sqrt(0.5*(2.0*i+1.0)),OBMatrices::get().A_GR[size][i],1);
#endif
}
free(tmp);
}
- *mat = BasisA_GR[size];
+ *mat = OBMatrices::get().A_GR[size];
return;
}
/* computes the values of b dependent modes at either the gauss or
gauss radau points depending on the function and stores the result
- in BasisB_G array */
+ in OBMatrices::get().B_G array */
static double ***get_Bmodes(int size, double *z);
static double ***get_Cmodes(int size, double *z);
void get_modb_G (int size, double ****mat){
- if (!BasisB_G[size]){
+ if (!OBMatrices::get().B_G[size]){
double *w,*z;
getzw (size,&z,&w,'h');
- BasisB_G[size] = get_Bmodes(size,z);
+ OBMatrices::get().B_G[size] = get_Bmodes(size,z);
}
- *mat = BasisB_G[size];
+ *mat = OBMatrices::get().B_G[size];
return;
}
void get_modb_GR (int size, double ****mat){
- if (!BasisB_GR[size]){
+ if (!OBMatrices::get().B_GR[size]){
double *w,*z;
getzw (size,&z,&w,'b');
- BasisB_GR[size] = get_Bmodes(size,z);
+ OBMatrices::get().B_GR[size] = get_Bmodes(size,z);
}
- *mat = BasisB_GR[size];
+ *mat = OBMatrices::get().B_GR[size];
return;
}
static double ***get_Bmodes(int size, double *z){
@@ -792,27 +855,27 @@ static double ***get_Bmodes(int size, do
}
void get_modc_GR (int size, double ****mat){
- if (!BasisC_GR[size]){
+ if (!OBMatrices::get().C_GR[size]){
double *w,*z;
getzw (size,&z,&w,'c');
- BasisC_GR[size] = get_Cmodes(size,z);
+ OBMatrices::get().C_GR[size] = get_Cmodes(size,z);
}
- *mat = BasisC_GR[size];
+ *mat = OBMatrices::get().C_GR[size];
return;
}
void get_modc_G (int size, double ****mat){
- if (!BasisC_G[size]){
+ if (!OBMatrices::get().C_G[size]){
double *w,*z;
getzw (size,&z,&w,'h');
- BasisC_G[size] = get_Cmodes(size,z);
+ OBMatrices::get().C_G[size] = get_Cmodes(size,z);
}
- *mat = BasisC_G[size];
+ *mat = OBMatrices::get().C_G[size];
return;
}
@@ -884,31 +947,75 @@ void reset_bases(){
Hex_reset_basis();
}
-static double ***Mod_BasisA, ****Mod_BasisB, ****Mod_BasisC;
-
-void init_mod_basis(void){
- Mod_BasisA = (double ***) calloc(QGmax, sizeof(double **));
- Mod_BasisA -= 1;
+// This class initialises the modified basis matrices based on QGmax.
+// These arrays may be accessed from [1-QGmax]
+class MBMatrices {
+public:
+ // This is a singleton class. Get the instance.
+ static MBMatrices & get();
+
+ // Some compilers don't like the d-tor being private.
+ ~MBMatrices();
+
+ double ***A;
+ double ****B;
+
+private:
+ // Make the c-tor private so we can control how many objects
+ // are instantiated.
+ MBMatrices();
+};
+
+MBMatrices & MBMatrices::get()
+{
+ static MBMatrices singleton;
+ return singleton;
+}
+
+MBMatrices::MBMatrices()
+{
+ A = (double ***) calloc(QGmax, sizeof(double **));
+ A -= 1;
- Mod_BasisB = (double ****)calloc(QGmax, sizeof(double ***));
- Mod_BasisB -= 1;
+ B = (double ****)calloc(QGmax, sizeof(double ***));
+ B -= 1;
+}
- Mod_BasisC = (double ****)calloc(QGmax, sizeof(double ***));
- Mod_BasisC -= 1;
+MBMatrices::~MBMatrices()
+{
+ // i starts at 1 because of the offset, above.
+ for (int i = 1; i <= QGmax; ++i) {
+ if (A[i]) {
+ free_dmatrix(A[i], 0, 0);
+ }
+ if (B[i]) {
+ free(B[i][0][0]); free(B[i][0]); free(B[i]);
+ }
+ }
+ // The increment is needed to return to the address of the allocation.
+ free(++A);
+ free(++B);
}
+// Now that we use a singleton class to store the Matrices, this function is
+// redundant.
+void init_mod_basis()
+{}
+
/* Computes the values of modified a_mode[mode_a][j], mode_a=0..LGmax,
j=0..size-1, at size Gauss-Lob. points, stores the result in
- Mod_BasisA[size][mode_a] array
+ MBMatrics::get().A[size][mode_a] array
(1+a)/2
(1-a)/2
(1-a)/2 (1+a)/2
(1-a)/2 (1+a)/2 P_1^{1,1}(a)
(1-a)/2 (1+a)/2 P_{l-2}^{1,1}(a)*/
-void get_modA(int size, double ***mat){
+void get_modA(int size, double ***mat)
+{
+ double *** Mod_BasisA = MBMatrices::get().A;
if (!Mod_BasisA[size]){
double *w, *z;
int i;
@@ -940,9 +1047,12 @@ void get_modA(int size, double ***mat){
return;
}
static double ***get_Mod_Bmodes(int size, double *z);
-void get_modB(int size, double ****mat){
+
+void get_modB(int size, double ****mat)
+{
+ double **** Mod_BasisB = MBMatrices::get().B;
if (!Mod_BasisB[size]){
double *w,*z;
getzw (size,&z,&w,'b');
@@ -1010,9 +1120,33 @@ static double ***get_Mod_Bmodes(int size
return mat;
}
-static Basis ****Tri_bbase=NULL;
+struct Tri_bbase_singleton {
+ static Tri_bbase_singleton & get();
+ ~Tri_bbase_singleton();
+
+ Basis **** Tri_bbase;
+private:
+ Tri_bbase_singleton();
+};
+Tri_bbase_singleton & Tri_bbase_singleton::get()
+{
+ static Tri_bbase_singleton singleton;
+ return singleton;
+}
+
+Tri_bbase_singleton::Tri_bbase_singleton()
+ : Tri_bbase(0)
+{}
+
+Tri_bbase_singleton::~Tri_bbase_singleton()
+{
+
+ Tri_reset_basis();
+}
+
+
static Basis ****Quad_bbase=NULL;
static Basis *Tet_bbinf=NULL,*Tet_bbase=NULL;
static Basis *Pyr_bbinf=NULL,*Pyr_bbase=NULL;
static Basis *Prism_bbinf=NULL,*Prism_bbase=NULL;
@@ -1045,17 +1179,22 @@ Function Notes:
Basis *Tri::getbasis(){
- int i,j;
- if(!Tri_bbase){
- Tri_bbase = (Basis****) calloc(LGmax+1, sizeof(Basis***));
+ if(!Tri_bbase_singleton::get().Tri_bbase){
+ Tri_bbase_singleton::get().Tri_bbase =
+ (Basis****) calloc(LGmax+1, sizeof(Basis***));
+
+ Basis **** Tri_bbase = Tri_bbase_singleton::get().Tri_bbase;
+ int i,j;
+
for(i=0;i<LGmax+1;++i){
Tri_bbase[i] = (Basis***) calloc(QGmax+1, sizeof(Basis**));
for(j=0;j<QGmax+1;++j)
Tri_bbase[i][j] = (Basis**) calloc(QGmax+1, sizeof(Basis*));
}
}
+ Basis **** Tri_bbase = Tri_bbase_singleton::get().Tri_bbase;
if(!Tri_bbase[lmax][qa][qb])
Tri_bbase[lmax][qa][qb] = Tri_addbase(lmax,qa,qb,qc);
return Tri_bbase[lmax][qa][qb];
}
@@ -1322,32 +1461,31 @@ static Basis *Hex_addbase(int L, int qa,
return b;
}
-void Tri_reset_basis(void){
+void Tri_reset_basis(Basis *B);
+
+void Tri_reset_basis(void)
+{
+ Basis **** Tri_bbase = Tri_bbase_singleton::get().Tri_bbase;
+ if(!Tri_bbase)
+ return;
int i,j,k;
- if(Tri_bbase){
- for(i=0;i<LGmax+1;++i)
- for(j=0;j<QGmax+1;++j)
- for(k=0;k<QGmax+1;++k)
- if(Tri_bbase[i][j][k]){
- free(Tri_bbase[i][j][k]->vert[2].a);
- free(Tri_bbase[i][j][k]->vert[2].b);
-
- free(Tri_bbase[i][j][k]->vert);
- free(Tri_bbase[i][j][k]->edge);
- free(Tri_bbase[i][j][k]->face);
- free(Tri_bbase[i][j][k]);
- }
- for(i=0;i<LGmax+1;++i){
- for(j=0;j<QGmax+1;++j)
- free(Tri_bbase[i][j]);
- free(Tri_bbase[i]);
- }
- free(Tri_bbase);
- Tri_bbase = NULL;
+
+ for(i=0;i<LGmax+1;++i)
+ for(j=0;j<QGmax+1;++j)
+ for(k=0;k<QGmax+1;++k)
+ Tri_reset_basis(Tri_bbase[i][j][k]);
+
+ for(i=0;i<LGmax+1;++i){
+ for(j=0;j<QGmax+1;++j)
+ free(Tri_bbase[i][j]);
+ free(Tri_bbase[i]);
}
+ free(Tri_bbase);
+ Tri_bbase = NULL;
+
}
void Tri_reset_basis(Basis *B){
@@ -1356,8 +1494,12 @@ void Tri_reset_basis(Basis *B){
free(B->vert[2].b);
free(B->vert);
free(B->edge);
+ if (B->face) {
+ for(int i = 0; i < NTri_faces; ++i)
+ free(B->face[i]);
+ }
free(B->face);
free(B);
}
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Boundary.C clean2/Hlib/src/Boundary.C
--- dirty/Hlib/src/Boundary.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Boundary.C Wed May 1 17:51:41 2002
@@ -1962,10 +1962,10 @@ void Tri::Add_Surface_Contrib(Element *E
Element *E = this;
double *za,*zb,*wa,*wb;
double **ima, **imb;
double one = 1.0e0 , two = 2.0e0;
- double *f = Tri_wk;
- double *fi = Tri_wk+QGmax;
+ double *f = Tri_wk.get();
+ double *fi = f+QGmax;
getzw(qa,&za,&wa,'a'); // ok
getzw(qb,&zb,&wb,'b'); // ok
@@ -2796,60 +2796,61 @@ Function Notes:
void Tri::fill_edges(double *ux, double *uy, double *){
double **im;
// one field case
+ double * wk = Tri_wk.get();
if(!uy){
getim(qa,edge[0].qedg,&im,a2g);
- GetFace(ux,0,Tri_wk);
- Interp(*im,Tri_wk,qa,edge[0].h,edge[0].qedg);
+ GetFace(ux,0,wk);
+ Interp(*im,wk,qa,edge[0].h,edge[0].qedg);
getim(qa,edge[1].qedg,&im,b2g);
- GetFace(ux,1,Tri_wk);
- Interp(*im,Tri_wk,qb,edge[1].h,edge[1].qedg);
+ GetFace(ux,1,wk);
+ Interp(*im,wk,qb,edge[1].h,edge[1].qedg);
getim(qa,edge[2].qedg,&im,b2g);
- GetFace(ux,2,Tri_wk);
- Interp(*im,Tri_wk,qb,edge[2].h,edge[2].qedg);
+ GetFace(ux,2,wk);
+ Interp(*im,wk,qb,edge[2].h,edge[2].qedg);
}
else{
- double *wka = Tri_wk+QGmax;
- double *wkb = Tri_wk+2*QGmax;
+ double *wka = wk+QGmax;
+ double *wkb = wk+2*QGmax;
int qedg;
Edge *e;
e = edge;
qedg = e->qedg;
getim(qa,qedg,&im,a2g);
- GetFace(ux,0,Tri_wk);
- Interp(*im,Tri_wk,qa,wka,qedg);
+ GetFace(ux,0,wk);
+ Interp(*im,wk,qa,wka,qedg);
- GetFace(uy,0,Tri_wk);
- Interp(*im,Tri_wk,qa,wkb,qedg);
+ GetFace(uy,0,wk);
+ Interp(*im,wk,qa,wkb,qedg);
dvmul (qedg, e->norm->x, 1, wka, 1, e->h, 1);
dvvtvp(qedg, e->norm->y, 1, wkb, 1, e->h, 1, e->h, 1);
e = edge+1;
qedg = e->qedg;
getim(qb,qedg,&im,b2g);
- GetFace(ux,1,Tri_wk);
- Interp(*im,Tri_wk,qb,wka,qedg);
+ GetFace(ux,1,wk);
+ Interp(*im,wk,qb,wka,qedg);
- GetFace(uy,1,Tri_wk);
- Interp(*im,Tri_wk,qb,wkb,qedg);
+ GetFace(uy,1,wk);
+ Interp(*im,wk,qb,wkb,qedg);
dvmul (qedg, e->norm->x, 1, wka, 1, e->h, 1);
dvvtvp(qedg, e->norm->y, 1, wkb, 1, e->h, 1, e->h, 1);
e = edge+2;
qedg = e->qedg;
getim(qb,qedg,&im,b2g);
- GetFace(ux,2,Tri_wk);
- Interp(*im,Tri_wk,qb,wka,qedg);
+ GetFace(ux,2,wk);
+ Interp(*im,wk,qb,wka,qedg);
- GetFace(uy,2,Tri_wk);
- Interp(*im,Tri_wk,qb,wkb,qedg);
+ GetFace(uy,2,wk);
+ Interp(*im,wk,qb,wkb,qedg);
dvmul (qedg, e->norm->x, 1, wka, 1, e->h, 1);
dvvtvp(qedg, e->norm->y, 1, wkb, 1, e->h, 1, e->h, 1);
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/ChangeLog clean2/Hlib/src/ChangeLog
--- dirty/Hlib/src/ChangeLog Thu Jan 1 01:00:00 1970
+++ clean2/Hlib/src/ChangeLog Wed May 1 17:51:41 2002
@@ -0,0 +1,13 @@
+2002-10-04 Angus Leeming <[EMAIL PROTECTED]>
+
+ * Basis.C (OBMatrices): a new singleton class. Used to encapsulate the
+ orthogonal basis matrices BasisA_G etc. No other change of code at all;
+ E.g., matrix BasisA_G is now accessed as OBMatrices::get().A_G.
+ Why do this? Because the allocated memory is now cleaned-up when the
+ instance of OBMatrices goes out of scope at program end, removing
+ a heap of valgrind error messages.
+
+ * Tri.C (Tri_setup_iprodlap): free Tri * "T" and its member pointers.
+
+ * Memory.C (Tri::Mem_shift): check whether h and face[0].hj have been
+ allocated memory already and, if so, free it.
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Element.C clean2/Hlib/src/Element.C
--- dirty/Hlib/src/Element.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Element.C Wed May 1 17:51:41 2002
@@ -74,8 +74,32 @@ Element::Element(const Element &E){
face = (Face *)calloc(Nfaces,sizeof(Face));
memcpy (face,E.face,Nfaces*sizeof(Face));
}
+
+Element::~Element()
+{
+ free(vert);
+ free(edge);
+ free(h);
+
+ if (face && face[0].hj)
+ free(face[0].hj);
+
+ free(face);
+
+ // Looks to me like curve, curvX and geom should be stored by a
+ // boost::shared_ptr, as multiple elements reference the same block of memory.
+ // Angus
+// double **h;
+// double ***h_3d;
+// double ***hj_3d;
+// Curve *curve;
+// Cmodes *curvX;
+// Geom *geom;
+// Element *next;
+}
+
void Element::EdgeJbwd(double *d, int edg){
Basis *b = getbasis();
int va = vnum(edg,0);
int vb = vnum(edg,1), i;
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Element_List.C clean2/Hlib/src/Element_List.C
--- dirty/Hlib/src/Element_List.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Element_List.C Wed May 1 17:51:41 2002
@@ -63,8 +63,22 @@ Element_List::Element_List(Element **hea
nztot = 1;
flevels = (Element_List**)NULL;
}
+Element_List::~Element_List()
+{
+ for (int i = 0; i < nel; ++i) {
+ delete flist[i];
+ }
+ free(flist);
+
+ if (base_h)
+ free(base_h);
+ if (base_hj)
+ free(base_hj);
+}
+
+
Element *Element_List::operator()(int i){
#ifdef DEBUG
if(i>nel-1)
fprintf(stderr,"Element_List::operator() bounds error\n");
@@ -80,8 +94,22 @@ Element *Element_List::operator[](int i)
return flist[i];
}
void Element_List::Cat_mem(){
+ if (base_h) {
+ fprintf(stderr, "Element_List::Cat_mem(): base_h is set already: %d\n",
+ base_h);
+ free(base_h);
+ base_h = 0;
+ }
+
+ if (base_hj) {
+ fprintf(stderr, "Element_List::Cat_mem(): base_hj is set already: %d\n",
+ base_hj);
+ free(base_hj);
+ base_hj = 0;
+ }
+
htot = 0;
hjtot = 0;
Eloop(fhead){
htot += el->qtot;
@@ -90,20 +118,19 @@ void Element_List::Cat_mem(){
base_h = dvector(0, nz*htot-1);
base_hj = dvector(0, nz*hjtot-1);
+ dzero(nz*htot, base_h, 1);
+ dzero(nz*hjtot, base_hj, 1);
+
double *tmp_base_h = base_h;
double *tmp_base_hj = base_hj;
- dzero(nz*htot, base_h, 1);
- dzero(nz*hjtot, base_hj, 1);
-
Eloop(fhead){
el->Mem_shift(tmp_base_h, tmp_base_hj, 'n');
tmp_base_h += el->qtot;
tmp_base_hj += el->Nmodes;
}
-
}
void Element_List::Trans(Element_List *EL, Nek_Trans_Type ntt){
// Treat FFT as a global operation
@@ -1199,9 +1226,8 @@ static void add_pf(int elmt, int side);
Element_List *Part_elmt_list(Element_List *GEL, char trip){
int i,k;
int lnel = pllinfo.nloop, ig;
- Element_List *EL = (Element_List*) new Element_List();
Element** new_E = (Element**) malloc(lnel*sizeof(Element*));
Element *Eg, *E;
// fprintf(stderr, "Entering: Part_elmt_list\n");
@@ -1449,9 +1475,10 @@ Element_List *Part_elmt_list(Element_Lis
}
}
free(inc);
- EL = (Element_List*) new Element_List(new_E, lnel);
+ Element_List * EL = (Element_List*) new Element_List(new_E, lnel);
+
EL->Cat_mem();
for(E=EL->fhead;E;E=E->next){
E->geom = (Geom*)0;
E->set_geofac();
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Geofac.C clean2/Hlib/src/Geofac.C
--- dirty/Hlib/src/Geofac.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Geofac.C Wed May 1 17:51:41 2002
@@ -77,11 +77,11 @@ void Tri::set_edge_geofac(){
Edge *e;
Vert *v;
double **mat,*wk1,**im,nx,ny,fac;
- double *wk = Tri_wk;
+ double *wk = Tri_wk.get();
- wk1 = Tri_wk+QGmax;
+ wk1 = wk+QGmax;
v = vert;
for(i = 0; i < Nedges; ++i){
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Gradient.C clean2/Hlib/src/Gradient.C
--- dirty/Hlib/src/Gradient.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Gradient.C Wed May 1 17:51:41 2002
@@ -236,9 +236,9 @@ void Tri::Grad_h(double *H, double *ux,
Geom *G = geom;
Mode *v = getbasis()->vert;
double *u = H; // fix ??
// double *wk = dvector(0,QGmax*QGmax-1);
- double *ur = Tri_wk;
+ double *ur = Tri_wk.get();
uz=uz; // compiler fix
getD(&da,&dat,&db,&dbt,NULL,NULL);
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/InnerProd.C clean2/Hlib/src/InnerProd.C
--- dirty/Hlib/src/InnerProd.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/InnerProd.C Wed May 1 17:51:41 2002
@@ -103,8 +103,10 @@ void Tri::Iprod_d(Element *T, Basis *ba,
Edge *e = T->edge;
Face *f = T->face;
size_t d = sizeof(double);
+ double *wk = Tri_wk.get();
+
T->state = 't';
getzw(T->qa,&wa,&wa,'a');
getzw(T->qb,&wb,&wb,'b');
@@ -117,24 +119,24 @@ void Tri::Iprod_d(Element *T, Basis *ba,
for(i = 0; i < qb; ++i) dvmul(qa,wa,1, T->h[i],1, T->h[i],1);
for(i = 0; i < qb; ++i) dscal(qa,wb[i],T->h[i],1);
- dgemm('T','N',qb,lmax,qa,1.0,T->h[0],qa,ba->vert[1].a,qa,0.0,Tri_wk,qb);
+ dgemm('T','N',qb,lmax,qa,1.0,T->h[0],qa,ba->vert[1].a,qa,0.0,wk,qb);
#if 0
H = *T->h;
- dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,Tri_wk ,1,0.0,H ,1);
- dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,Tri_wk+qb,1,0.0,H+lmax,1);
+ dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,wk ,1,0.0,H ,1);
+ dgemv('T',qb,lmax,1.0,bb->vert[2].b,qb,wk+qb,1,0.0,H+lmax,1);
H += 2*lmax;
- Tri_wk_a = Tri_wk+2*qb;
- for(i = 0; i < lmax-2; Tri_wk_a+=qb,H+=lmax-2-i, ++i)
- dgemv('T',qb,lmax-2-i,1.0,bb->edge[0][i].b,qb,1,Tri_wk_a,1,0.0,H,1);
+ wk_a = wk+2*qb;
+ for(i = 0; i < lmax-2; wk_a+=qb,H+=lmax-2-i, ++i)
+ dgemv('T',qb,lmax-2-i,1.0,bb->edge[0][i].b,qb,1,wk_a,1,0.0,H,1);
#else
H = *T->h;
- mxva(bb->vert[2].b,qb,1,Tri_wk ,1,H ,1,lmax,qb);
- mxva(bb->vert[2].b,qb,1,Tri_wk+qb,1,H+lmax,1,lmax,qb);
+ mxva(bb->vert[2].b,qb,1,wk ,1,H ,1,lmax,qb);
+ mxva(bb->vert[2].b,qb,1,wk+qb,1,H+lmax,1,lmax,qb);
H += 2*lmax;
- wk_a = Tri_wk+2*qb;
+ wk_a = wk+2*qb;
for(i = 0; i < lmax-2; wk_a+=qb,H+=lmax-2-i, ++i)
mxva(bb->edge[0][i].b,qb,1,wk_a,1,H,1,lmax-2-i,qb);
#endif
@@ -974,11 +976,13 @@ void Tri::HelmHoltz(Metric *lambda){
Basis *b,*db;
double *store,*u1,*u2;
static int SVV = dparam("SVVfactor") ? 1:0;
- store = Tri_wk+QGmax*QGmax;
- u1 = Tri_wk+2*QGmax*QGmax;
- u2 = Tri_wk+3*QGmax*QGmax;
+ double *wk = Tri_wk.get();
+
+ store = wk+QGmax*QGmax;
+ u1 = wk+2*QGmax*QGmax;
+ u2 = wk+3*QGmax*QGmax;
Nm = Nmodes;
b = getbasis();
db = derbasis();
@@ -1612,9 +1616,9 @@ void Tri::form_diprod(double *u1, double
{
register int i;
const int qt = qa*qb;
Geom *G = geom;
- double *tmp = Tri_wk+4*QGmax*QGmax; // sort out later
+ double *tmp = Tri_wk.get()+4*QGmax*QGmax; // sort out later
double *Ur = tmp;
// Grad_d(u1,u2,(double*)NULL,'a');
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Interp_point.C clean2/Hlib/src/Interp_point.C
--- dirty/Hlib/src/Interp_point.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Interp_point.C Wed May 1 17:51:41 2002
@@ -13,8 +13,9 @@
#include <math.h>
#include <veclib.h>
#include "polylib.h"
#include "hotel.h"
+#include <smart_ptr.hpp>
static int Find_elmt_coords_2d(Element_List *U,
double xo, double yo,
int *eid, double *a, double *b, char trip);
@@ -207,9 +208,9 @@ double eval_field_at_pt_3d(int qa, int q
}
-class BoundBox{
+class BoundBox {
public:
int nel;
int dim;
double *lx;
@@ -220,31 +221,41 @@ public:
double *uy;
double *uz;
};
-BoundBox *Bbox = NULL;
+namespace {
+
+struct DestroyBoundBox {
+ void operator()(BoundBox * Bbox)
+ {
+ free(Bbox->lx);
+ free(Bbox->ly);
+ free(Bbox->ux);
+ free(Bbox->uy);
+ if (Bbox->dim == 3) {
+ free(Bbox->lz);
+ free(Bbox->uz);
+ }
+ free (Bbox);
+ Bbox = 0;
+ }
+};
+
+} // namespace anon
+
+nektar::scoped_ptr<BoundBox, DestroyBoundBox> Bbox;
#define TOLBOX 1.1 // scaling factor of box
+
void setup_elmt_boxes(Element_List *U){
Element *E;
double mean;
int i,k;
if(U->fhead->dim() == 3){
- if(Bbox){
- free(Bbox->lx);
- free(Bbox->ly);
- free(Bbox->lz);
-
- free(Bbox->ux);
- free(Bbox->uy);
- free(Bbox->uz);
-
- free(Bbox);
- }
+ Bbox.reset((BoundBox*) calloc(1, sizeof(BoundBox)));
- Bbox = (BoundBox*) calloc(1, sizeof(BoundBox));
Bbox->dim = U->fhead->dim();
Bbox->nel = U->nel;
Bbox->lx = dvector(0, U->nel-1);
Bbox->ly = dvector(0, U->nel-1);
@@ -257,9 +268,9 @@ void setup_elmt_boxes(Element_List *U){
Coord X;
X.x = dvector(0, QGmax*QGmax*QGmax-1);
X.y = dvector(0, QGmax*QGmax*QGmax-1);
X.z = dvector(0, QGmax*QGmax*QGmax-1);
-
+
for(k=0;k<U->nel;++k){
E = U->flist[k];
E->coord(&X);
Bbox->lx[k] = X.x[idmin(E->qtot, X.x, 1)];
@@ -297,19 +308,11 @@ void setup_elmt_boxes(Element_List *U){
free(X.y);
free(X.z);
}
else{
- if(Bbox){
- free(Bbox->lx);
- free(Bbox->ly);
-
- free(Bbox->ux);
- free(Bbox->uy);
- free(Bbox);
- }
+ Bbox.reset((BoundBox*) calloc(1, sizeof(BoundBox)));
- Bbox = (BoundBox*) calloc(1, sizeof(BoundBox));
Bbox->dim = U->fhead->dim();
Bbox->nel = U->nel;
Bbox->lx = dvector(0, U->nel-1);
Bbox->ly = dvector(0, U->nel-1);
@@ -349,19 +352,19 @@ void setup_elmt_boxes(Element_List *U){
free(X.y);
}
}
-int point_in_box_2d(double x, double y, int boxid){
-
+int point_in_box_2d(double x, double y, int boxid)
+{
if( x <= Bbox->ux[boxid] && x >= Bbox->lx[boxid] &&
y <= Bbox->uy[boxid] && y >= Bbox->ly[boxid])
return 1;
return 0;
}
-int point_in_box_3d(double x, double y, double z, int boxid){
-
+int point_in_box_3d(double x, double y, double z, int boxid)
+{
if( x <= Bbox->ux[boxid] && x >= Bbox->lx[boxid] &&
y <= Bbox->uy[boxid] && y >= Bbox->ly[boxid] &&
z <= Bbox->uz[boxid] && z >= Bbox->lz[boxid])
return 1;
@@ -372,43 +375,48 @@ int point_in_box_3d(double x, double y,
#define MAXIT 101 /* maximum iterations for convergence */
#define EPSILON 1.e-8 /* |x-xp| < EPSILON error tolerance */
#define GETRSDIV 15.0 /* |r,s| > GETRSDIV stop the search */
-
-
-Coord **meshX = NULL;
int meshX_nel;
int meshX_dim;
-void set_mesh_coord(Element_List *EL){
- Element *E;
- if(!meshX){
- meshX = (Coord**) calloc(EL->nel, sizeof(Coord*));
-
- for(E=EL->fhead;E;E=E->next){
- meshX[E->id] = (Coord*) calloc(1, sizeof(Coord));
- meshX[E->id]->x = dvector(0, E->qtot-1);
- meshX[E->id]->y = dvector(0, E->qtot-1);
- if(E->dim() == 3)
- meshX[E->id]->z = dvector(0, E->qtot-1);
- E->coord(meshX[E->id]);
- }
- meshX_nel = EL->nel;
- meshX_dim = EL->fhead->dim();
- }
- else{
- int k;
- for(k=0;k<EL->nel;++k){
- free(meshX[k]->x);
- free(meshX[k]->y);
+namespace {
+
+struct DestroyMeshX {
+ void operator()(Coord ** ptr)
+ {
+ for(int i = 0; i < meshX_nel; ++i) {
+ free(ptr[i]->x);
+ free(ptr[i]->y);
if(meshX_dim == 3)
- free(meshX[k]->z);
- free(meshX[k]);
+ free(ptr[i]->z);
+ free(ptr[i]);
}
- free(meshX);
- meshX = NULL;
- set_mesh_coord(EL);
+ free(ptr);
+ ptr = 0;
+ }
+};
+
+} // namespace anon
+
+nektar::scoped_ptr<Coord *, DestroyMeshX> meshX_ptr;
+
+void set_mesh_coord(Element_List *EL){
+ Element *E;
+
+ meshX_ptr.reset((Coord**) calloc(EL->nel, sizeof(Coord*)));
+ Coord ** meshX = meshX_ptr.get();
+
+ for(E=EL->fhead;E;E=E->next){
+ meshX[E->id] = (Coord*) calloc(1, sizeof(Coord));
+ meshX[E->id]->x = dvector(0, E->qtot-1);
+ meshX[E->id]->y = dvector(0, E->qtot-1);
+ if(E->dim() == 3)
+ meshX[E->id]->z = dvector(0, E->qtot-1);
+ E->coord(meshX[E->id]);
}
+ meshX_nel = EL->nel;
+ meshX_dim = EL->fhead->dim();
}
int Find_coords_2d(Element *E, double xo, double yo,
double *a, double *b, char trip){
@@ -440,8 +448,10 @@ int Find_coords_2d(Element *E, double xo
hr = dvector (0, QGmax-1);
hs = dvector (0, QGmax-1);
}
+ Coord ** meshX = meshX_ptr.get();
+
gX = meshX[E->id];
gridx = gX->x;
gridy = gX->y;
@@ -618,8 +628,10 @@ int Find_coords_3d(Element *E, double xo
hs = dvector (0, QGmax-1);
ht = dvector (0, QGmax-1);
}
+ Coord ** meshX = meshX_ptr.get();
+
gX = meshX[E->id];
gridx = gX->x;
gridy = gX->y;
gridz = gX->z;
@@ -911,9 +923,9 @@ void Prt_find_local_coords(Element_List
int *Eid, Coord *A, char status){
int i;
int dim = U->fhead->dim();
- if(!Bbox){
+ if(!Bbox.get()){
setup_elmt_boxes(U);
set_mesh_coord(U);
}
@@ -1162,9 +1174,9 @@ void Find_local_coords(Element_List *U,
int i;
int dim = U->fhead->dim();
- if(!Bbox){
+ if(!Bbox.get()){
setup_elmt_boxes(U);
set_mesh_coord(U);
}
@@ -1202,9 +1214,9 @@ void find_local_coords(Element_List *U,
int i;
int dim = U->fhead->dim();
- if(!Bbox){
+ if(!Bbox.get()){
setup_elmt_boxes(U);
set_mesh_coord(U);
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Memory.C clean2/Hlib/src/Memory.C
--- dirty/Hlib/src/Memory.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Memory.C Wed May 1 17:51:41 2002
@@ -18,11 +18,13 @@
#include <math.h>
#include <polylib.h>
#include "veclib.h"
#include "hotel.h"
+#include <smart_ptr.hpp>
- double *Tri_wk = (double*)0;
+
+nektar::scoped_c_ptr<double> Tri_wk;
double *Quad_wk = (double*)0;
double *Quad_Jbwd_wk = (double*)0;
double *Quad_Iprod_wk = (double*)0;
@@ -61,11 +63,9 @@
double *Hex_InterpToGaussFace_wk = (double*)0;
void Tri_work(){
- if(Tri_wk)
- free(Tri_wk);
- Tri_wk = dvector(0,5*QGmax*QGmax-1);
+ Tri_wk.reset(dvector(0,5*QGmax*QGmax-1));
}
void Quad_work(){
@@ -608,9 +608,14 @@ Function Notes:
*/
void Tri::Mem_shift(double *new_h, double *new_hj, char Trip){
+ if (h)
+ free(h);
h = (double**) calloc(qb,sizeof(double*));
+
+ if (face[0].hj)
+ free(face[0].hj);
face[0].hj = (double**) calloc(face[0].l,sizeof(double*));
// Physical memory
int i;
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Transforms.C clean2/Hlib/src/Transforms.C
--- dirty/Hlib/src/Transforms.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Transforms.C Wed May 1 17:51:41 2002
@@ -249,13 +249,15 @@ void Tri::Jbwd(Element *T, Basis *b){
for(i = 0; i < lmax-2; H+=lmax-2-i, ++i)
dgemv('N',Qb,lmax-2-i,1.0,b->edge[0][i].b,Qb,H,1,0.0,wk+2+i,lmax);
dgemm('N','N',Qa,Qb,lmax,1.0,b->vert[1].a,Qa,wk,lmax,0.0,T->h[0],Qa);
#else /* mxv and mxm vertions */
- mxva(b->vert[2].b,1,Qb,H,1,Tri_wk ,lmax,Qb,lmax); H+=lmax;
- mxva(b->vert[2].b,1,Qb,H,1,Tri_wk+1,lmax,Qb,lmax); H+=lmax;
+ double *wk = Tri_wk.get();
+
+ mxva(b->vert[2].b,1,Qb,H,1,wk ,lmax,Qb,lmax); H+=lmax;
+ mxva(b->vert[2].b,1,Qb,H,1,wk+1,lmax,Qb,lmax); H+=lmax;
for(i = 0; i < lmax-2; H+=lmax-2-i, ++i)
- mxva(b->edge[0][i].b,1,Qb,H,1,Tri_wk+2+i,lmax,Qb,lmax-2-i);
- mxm(Tri_wk,Qb,b->vert[1].a,lmax,T->h[0],Qa);
+ mxva(b->edge[0][i].b,1,Qb,H,1,wk+2+i,lmax,Qb,lmax-2-i);
+ mxm(wk,Qb,b->vert[1].a,lmax,T->h[0],Qa);
#endif
}
@@ -1295,21 +1297,22 @@ void Tri::Obwd(double *in, double *out,
register int i;
double **ba,***bb;
int mode;
+ double *wk = Tri_wk.get();
get_moda_GL (qa, &ba);
get_modb_GR (qb, &bb);
mode = 0;
/* bwd trans w.r.t. b modes */
for (i= 0; i < L; ++i){
- mxva (bb[i][0],1,qb,in+mode,1,Tri_wk+i*qb,1,qb,L-i);
+ mxva (bb[i][0],1,qb,in+mode,1,wk+i*qb,1,qb,L-i);
mode += L-i;
}
/* bwd trans w.r.t. a modes */
- dgemm('N','T',qa,qb,L,1.0,*ba,qa,Tri_wk,qb,0.0,out,qa);
+ dgemm('N','T',qa,qb,L,1.0,*ba,qa,wk,qb,0.0,out,qa);
}
@@ -1464,8 +1467,10 @@ void Tri::Ofwd(double *in, double *out,
register int i;
double **ba,***bb,*wa,*wb;
int mode;
+ double *wk = Tri_wk.get();
+
getzw (qa,&wa,&wa,'a');
getzw (qb,&wb,&wb,'b');
get_moda_GL (qa, &ba);
get_modb_GR (qb, &bb);
@@ -1484,20 +1489,20 @@ void Tri::Ofwd(double *in, double *out,
for(i = 0; i < qa; ++i)
dscal(qb,wa[i],in+i,qa);
/* Inner product trans w.r.t. a modes */
- dgemm('T','N',qb,L,qa,1.0,in,qa,*ba,qa,0.0,Tri_wk,qb);
+ dgemm('T','N',qb,L,qa,1.0,in,qa,*ba,qa,0.0,wk,qb);
#ifdef NONORM
for(i=0;i<L;++i)
- dscal(qb, 0.5*(2.*i+1.), Tri_wk+qb*i, 1);
+ dscal(qb, 0.5*(2.*i+1.), wk+qb*i, 1);
double fac = 1.;
#endif
mode = 0;
for (i= 0; i < L; ++i){ /* Inner product trans w.r.t. b modes */
- dvmul (qb,wb,1,Tri_wk+i*qb,1,Tri_wk+i*qb,1);
- mxva (bb[i][0],qb,1,Tri_wk+i*qb,1,out+mode,1,L-i,qb);
+ dvmul (qb,wb,1,wk+i*qb,1,wk+i*qb,1);
+ mxva (bb[i][0],qb,1,wk+i*qb,1,out+mode,1,L-i,qb);
#ifdef NONORM
for(j = 0; j < L-i; ++j)
out[mode+j] *= fac*(i+j+1.);
fac *= 0.25;
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Hlib/src/Tri.C clean2/Hlib/src/Tri.C
--- dirty/Hlib/src/Tri.C Fri Feb 1 13:26:29 2002
+++ clean2/Hlib/src/Tri.C Wed May 1 17:51:41 2002
@@ -30,9 +30,10 @@ Tri::Tri(){
}
Tri::Tri(Element *E){
- if(!Tri_wk)
+
+ if(!Tri_wk.get())
Tri_work();
id = E->id;
type = E->type;
state = 'p';
@@ -74,9 +75,9 @@ Tri::Tri(Element *E){
Tri::Tri(int i_d, char ty, int L, int Qa, int Qb, int Qc, Coord *X){
int i;
- if(!Tri_wk)
+ if(!Tri_wk.get())
Tri_work();
id = i_d;
type = ty;
@@ -414,11 +415,16 @@ void Tri_setup_iprodlap(){
free_mvec(gb1[0]); free((char *) gb1);
free((char*)w);
free(tmp);
+
+ // These should really, really be free-d in the Element d-tor
+ // but I don't know enough to know what side effects writing this
+ // d-tor would have.
+ // Angus 11 April 2002.
+ free(T->vert);
+ free(T->edge);
+ free(T->face);
+ free(T);
}
#endif
-
-
-
-
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/ChangeLog clean2/Nektar2d/src/ChangeLog
--- dirty/Nektar2d/src/ChangeLog Thu Jan 1 01:00:00 1970
+++ clean2/Nektar2d/src/ChangeLog Wed May 1 17:51:41 2002
@@ -0,0 +1,24 @@
+2002-10-04 Angus Leeming <[EMAIL PROTECTED]>
+
+ * nektar.C (main): initilise "step".
+ Valgrind reports "Use of uninitialised value of size 4" in
+ writeHeader(_IO_FILE *, Field *) (../src/Fieldfiles.C:1486)
+ that propogates from here!
+
+ * analyser.C (Analyser): free the pointer "V" (but not its contents).
+
+ * bwoptim.C (free_Fctcon): add the all-imporatant braces to the for-loop
+ so that we actually free some memory!
+
+ * prepost.C (PreProcess): Don't free already free'd memory; wrap the
+ free(P->base_h) call inside an if-statement. Ditto for P->base_hj.
+
+ * io.C (ReadMesh): free Coord X's member dvectors.
+
+ * io.C (bsort): Fix one memory leak caused by not free-ing the array,
+ tmp.
+ Attempt to fix memory leaks caused by not free-ing the original
+ bndry_list, but ultimately fail. I fail here because memory is
+ allocated in one of two different ways and we cannot resolve which in
+ order to free it cleanly. See the routine for details.
+
\ No newline at end of file
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/analyser.C clean2/Nektar2d/src/analyser.C
--- dirty/Nektar2d/src/analyser.C Tue Jan 9 18:38:14 2001
+++ clean2/Nektar2d/src/analyser.C Wed May 1 17:51:41 2002
@@ -145,8 +145,11 @@ void Analyser (Domain *omega, int step,
V[0]->Set_state('t');
V[1]->Set_state('t');
// V[2]->Set_state('t');
+ // free the pointer
+ free(V);
+
return;
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/bwoptim.C clean2/Nektar2d/src/bwoptim.C
--- dirty/Nektar2d/src/bwoptim.C Tue Jun 29 11:03:26 1999
+++ clean2/Nektar2d/src/bwoptim.C Wed May 1 17:51:41 2002
@@ -345,15 +345,16 @@ void addfct(Fctcon *con, int *pts, int n
void free_Fctcon(Fctcon *con, int n){
register int i;
Facet *f,*f1;
- for(i = 0; i < n; ++i)
+ for(i = 0; i < n; ++i) {
f = con[i].f;
while(f){
f1 = f->next;
free(f);
f = f1;
}
+ }
free(con);
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/io.C clean2/Nektar2d/src/io.C
--- dirty/Nektar2d/src/io.C Tue Feb 13 19:05:10 2001
+++ clean2/Nektar2d/src/io.C Wed May 1 17:51:41 2002
@@ -271,8 +271,11 @@ Element_List *ReadMesh(FILE *fp, char *s
new_E[k] = new Nodal_Tri (k,'u', L,qa,qb,0, &X);
#endif
}
+ free(X.x);
+ free(X.y);
+
for(k = 0; k < nel-1; ++k) new_E[k]->next = new_E[k+1];
new_E[k]->next = (Element*) NULL;
Element_List* E_List = new Element_List(new_E, nel);
@@ -979,30 +982,53 @@ Bndry *bsort(Bndry *bndry_list, int nbcs
tmp = (Bndry**) malloc(nbcs*sizeof(Bndry*));
/* Copy the linked list into a regular array and sort it */
-
for(i = 0 ; i < nbcs; ++i, bedg = bedg->next) tmp[i] = bedg;
qsort(tmp, nbcs, sizeof(Bndry*), bedgcmp);
/* Create a new version of the ordered list */
-
bedg = (Bndry*) malloc (nbcs * sizeof(Bndry));
for(i = 0; i < nbcs; ++i) {
-// memcpy (bedg + i, tmp[i], sizeof(Bndry));
bedg[i] = *(tmp[i]);
bedg[i].id = i;
bedg[i].next = bedg + i + 1;
}
bedg[nbcs-1].next = (Bndry*) NULL;
- /* Free up the space occupied by the old list */
- bndry_list = bedg;
-// for(i = 0; i < nbcs; ++i) free (tmp[i]);
+ /* Free up the memory occupied by the temporary array */
+ free (tmp);
+
+ // Free up the memory occupied by the old list.
+ // Ahhhhhhhhhhhhhh! (Scream for help!)
+
+ // The memory allocated to the linked list bndry_list is allocated in
+ // one of two different ways.
-// free (tmp);
+ // In ReadBCs, above, the memory is allocated in piecewise fashion to each
+ // new element of the linked list. The following therefore frees the memory:
+
+ // (Commented out for now until a proper solution is decided upon.)
+// Bndry * ptr = bndry_list;
+// Bndry * next;
+// while(ptr) {
+// next = ptr->next;
+// free(ptr);
+// ptr = next;
+// }
+
+ // In BuildPBCs in pressure.C, the memory is allocated in one lump at the
+ // start. The following, therefore, frees the memory in this case.
+ free(bndry_list);
+
+ // Three solutions present themselves.
+ // 1. Use the STL list and leave such worries behind you.
+ // 2. Decide on a consistent memory allocation strategy.
+ // 3. Pass an addition parameter to bsort that flags how the memory should
+ // be free-d.
- return bndry_list;
+ // return the new, sorted list
+ return bedg;
}
/*
* Boundaries are sorted by their ASCII character values except for
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/nektar.C clean2/Nektar2d/src/nektar.C
--- dirty/Nektar2d/src/nektar.C Fri Jan 19 21:12:13 2001
+++ clean2/Nektar2d/src/nektar.C Wed May 1 17:51:41 2002
@@ -12,12 +12,16 @@
#include <time.h>
#include "nektar.h"
+//void free_Domain(Domain * omega);
+
main(int argc, char *argv[]){
- Domain *Omega; /* Solution Domain */
- int step, nsteps; /* Discrete time (step) */
- double dt,time = 0.0; /* Continuous time */
+ Domain *Omega; /* Solution Domain */
+ int step = 0; /* Discrete time (step) */
+ int nsteps; /* Number of steps */
+ double time = 0.0; /* Continuous time */
+ double dt;
#ifdef DEBUG
/* mallopt(M_DEBUG,1); */
init_debug();
@@ -31,6 +35,7 @@ main(int argc, char *argv[]){
Navier_Stokes(Omega,time,time+dt*nsteps);
PostProcess(Omega, step, time+dt*nsteps);
+ //free_Domain(Omega);
exit(0);
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Nektar2d/src/prepost.C clean2/Nektar2d/src/prepost.C
--- dirty/Nektar2d/src/prepost.C Tue Apr 17 14:04:34 2001
+++ clean2/Nektar2d/src/prepost.C Wed May 1 17:51:41 2002
@@ -196,9 +196,9 @@ void parse_args(int argc, char *argv[]){
* ------------------------------------------------------------------------ */
void CorrectRobinFluxes(Bndry *B);
static void Set_Oseen(FILE *rea_file,Bsystem *PB, Bsystem *B, Element_List *U);
-
+
Domain *PreProcess(int argc, char **argv)
{
FILE *rea_file;
Element_List *U,*P;
@@ -287,9 +287,16 @@ Domain *PreProcess(int argc, char **argv
P = U->gen_aux_field('p');
option_set("variable",1);
/* set pressure field to L-1 */
- free(P->base_h); free(P->base_hj);
+ if (P->base_h) {
+ free(P->base_h);
+ P->base_h = 0;
+ }
+ if (P->base_hj) {
+ free(P->base_hj);
+ P->base_hj = 0;
+ }
for(E=P->fhead;E;E=E->next){
for(i=0;i<E->Nedges;++i)
size[i] = max(0,E->edge[i].l-2);
@@ -411,9 +418,9 @@ Domain *PreProcess(int argc, char **argv
}
if(slvtype == TurbCurr){
omega->Lfx = U->gen_aux_field('l');
- omega->Lfy = U->gen_aux_field('l');
+ omega->Lfy = U->gen_aux_field('l');
omega->lfx = (double**) malloc(Je*sizeof(double*));
for(k = 0; k < Je; ++k){
omega->lfx[k] = dvector(0,U->htot*U->nz-1);
@@ -580,8 +587,9 @@ void PostProcess(Domain *omega, int step
Writefld(fp, omega->name, step, time, nfields, V);
fflush(fp[0]);
+ free(V);
free(size);
//fclose(fp[0]);
//fclose(omega->his_file);
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/ChangeLog clean2/Particle/src/ChangeLog
--- dirty/Particle/src/ChangeLog Thu Jan 1 01:00:00 1970
+++ clean2/Particle/src/ChangeLog Wed May 1 17:51:41 2002
@@ -0,0 +1,16 @@
+2002-04-23 Angus Leeming <[EMAIL PROTECTED]>
+
+ * PC_makestep.C (PC_Make_step_hybrid): initialise double shear.
+ Initialise matrices h and J. Together these resolve valgrind "Use of
+ uninitialised value" reports.
+
+ * particle.h: give class Butcher_arrays a default c-tor and d-tor.
+
+ * PC_accessories.C (~Butcher_arrays): free memory.
+
+ * PC_accessories.C (~Particle_Cluster): Vcsi was free-d and then
+ delete []-ed!
+
+ * Pre_interpolation.C (Prt_Interp_field_in_time): move creation
+ of dmatrix data closer to where it's used and after those return
+ statements.
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/PC_accessories.C clean2/Particle/src/PC_accessories.C
--- dirty/Particle/src/PC_accessories.C Fri Feb 1 17:14:03 2002
+++ clean2/Particle/src/PC_accessories.C Wed May 1 17:51:41 2002
@@ -8,12 +8,24 @@
/* These functions are accessory functions for class Particle_Cluster */
/* Their puorpose is to access the member functions of this class */
+Butcher_arrays::~Butcher_arrays()
+{
+ if (ai+1)
+ free_dmatrix(ai, 1, 0);
+ if (bi)
+ free_dmatrix(bi, 0, 0);
+ if (ci+1)
+ free_dvector(ci, 1);
+}
+
/* Default constructor */
-Particle_Cluster::Particle_Cluster(){
+Particle_Cluster::Particle_Cluster()
+ : nel(0), mesh(0)
+{
dim = 0;
npts = 0;
xyz = NULL;
csi = NULL;
@@ -53,9 +65,11 @@ Coord *alloc_coord(int npts, int dim){
//#define dealloc_coord(c) {delete [] c->x; delete [] c;}
#define dealloc_coord(c) {delete [] c->x; free (c);}
-Particle_Cluster::Particle_Cluster(int dimin, int nptsin){
+Particle_Cluster::Particle_Cluster(int dimin, int nptsin)
+ : nel(0), mesh(0)
+{
int p;
dim = dimin;
npts = nptsin;
@@ -85,16 +99,21 @@ Particle_Cluster::Particle_Cluster(int d
}
/* Destructor: free memory */
-Particle_Cluster::~Particle_Cluster(){
+Particle_Cluster::~Particle_Cluster()
+{
dealloc_coord(xyz);
dealloc_coord(csi);
dealloc_coord(Vxyz);
dealloc_coord(xyzstart);
dealloc_coord(Vcsi);
- delete [] Vcsi;
+ for(int i = 0; i < nel; ++i){
+ free(mesh[i].x);
+ }
+
+ delete [] mesh;
delete [] k;
delete [] colour;
delete [] s;
delete [] restime;
@@ -104,9 +123,11 @@ Particle_Cluster::~Particle_Cluster(){
}
/* Copy constructor */
-Particle_Cluster::Particle_Cluster(const Particle_Cluster *cluster){
+Particle_Cluster::Particle_Cluster(const Particle_Cluster *cluster)
+ : nel(0), mesh(0)
+{
int p;
//sort out copy constructor
@@ -165,8 +186,11 @@ void Particle_Cluster::PC_Check_MaxTime(
void Particle_Cluster::Set_ai_Butcher_array(){
double **ai;
+ if (aibici.ai)
+ fprintf(stderr, "ai is set already!\n");
+
if(RKstg > 1){
ai = dmatrix(1,RKstg-1,0,RKstg-2);
dzero((RKstg-1)*(RKstg-1),&ai[1][0],1);
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/PC_makestep.C clean2/Particle/src/PC_makestep.C
--- dirty/Particle/src/PC_makestep.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/PC_makestep.C Wed May 1 17:51:41 2002
@@ -203,8 +203,9 @@ void Particle_Cluster::PC_Make_step_hybr
Element *E;
#ifdef SHEAR
+ shear = 0.0;
nvar = dim + 1;
#else
nvar = dim;
#endif
@@ -223,8 +224,11 @@ void Particle_Cluster::PC_Make_step_hybr
}
h = dmatrix(0,dim-1,0,QGmax-1);
J = dmatrix(0,dim-1,0,dim-1);
+
+ dzero(dim * QGmax, h[0], 1);
+ dzero(dim * dim, J[0], 1);
Tol = dparam("P_COOR_TOL");
init = 0;
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/Prt_interpolation.C clean2/Particle/src/Prt_interpolation.C
--- dirty/Particle/src/Prt_interpolation.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/Prt_interpolation.C Wed May 1 17:51:41 2002
@@ -4,8 +4,9 @@
#include <veclib.h>
#include <hotel.h>
#include <gen_utils.h>
#include "particle.h"
+#include <smart_ptr.hpp>
void Particle_Cluster::Prt_Interp_geom_fact(Element *E, Coord A, double **J,
double **h){
if(dim == 2) {
@@ -111,10 +112,13 @@ void Particle_Cluster::Prt_Interp_field_
int nekdumps = fdata->nekdumps;
int *dataskip = fdata->dataskip;
double *dtime = fdata->dumptime, **data,dt;
Element *E = Vel[0]->flist[eid];
- static double last_tlevel = 0, *time, *mode;
- static int *Interp_done;
+
+ static double last_tlevel = 0;
+ static nektar::scoped_c_ptr<double> scoped_time;
+ static nektar::scoped_c_ptr<double> scoped_mode;
+ static nektar::scoped_c_ptr<int> scoped_Interp_done;
#ifdef ANALYTIC_SOL
for(i = 0; i < E->qa*E->qb*E->qc; ++i){
Vel[0]->flist[eid]->h_3d[0][0][i] = -mesh[eid].x[i];
@@ -123,18 +127,20 @@ void Particle_Cluster::Prt_Interp_field_
}
return;
#endif
+ if(!scoped_Interp_done.get()){
+ scoped_Interp_done.reset(ivector(0,Vel[0]->nel-1));
+ izero(Vel[0]->nel,scoped_Interp_done.get(),1);
- if(!Interp_done){
- Interp_done = ivector(0,Vel[0]->nel-1);
- izero(Vel[0]->nel,Interp_done,1);
- time = dvector(0,nload-1);
- mode = dvector(0,nload-1);
+ scoped_time.reset(dvector(0,nload-1));
+ scoped_mode.reset(dvector(0,nload-1));
}
- data = dmatrix(0,dim-1,0,E->Nmodes-1);
-
+ int * Interp_done = scoped_Interp_done.get();
+ double * time = scoped_time.get();
+ double * mode = scoped_mode.get();
+
// make sure tlevel is in nektar field dump interval
if(fdata->periodic)
tlevel = fmod(fmod(tlevel,fdata->periodic) + fdata->periodic,
fdata->periodic) + dtime[0];
@@ -171,8 +177,10 @@ void Particle_Cluster::Prt_Interp_field_
Field_Dump = i;
break;
}
+ data = dmatrix(0,dim-1,0,E->Nmodes-1);
+
if(Field_Dump+1)
copyelfield(fld_dump[Field_Dump][0],dim,eid,Vel);
else{ // interpolate
for(f=0;f < dim;++f){
@@ -186,8 +194,10 @@ void Particle_Cluster::Prt_Interp_field_
}
copyelfield(fld_dump[0][0],dim,eid,data,Vel);
}
+ free_dmatrix(data,0,0);
+
// put data in physical space
for(f=0;f<dim;++f)
Vel[f]->flist[eid]->Trans(Vel[f]->flist[eid], J_to_Q);
@@ -196,10 +206,8 @@ void Particle_Cluster::Prt_Interp_field_
else{
izero(Vel[0]->nel,Interp_done,1);
last_tlevel = tlevel;
}
-
- free_dmatrix(data,0,0);
return;
}
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/misc/ChangeLog clean2/Particle/src/misc/ChangeLog
--- dirty/Particle/src/misc/ChangeLog Thu Jan 1 01:00:00 1970
+++ clean2/Particle/src/misc/ChangeLog Wed May 1 17:51:41 2002
@@ -0,0 +1,5 @@
+2002-04-23 Angus Leeming <[EMAIL PROTECTED]>
+
+ * io.C (ReadMesh): initialise components of Coord X to prevent
+ valgrind reporting a subsequent use of an uninitialised value.
+ On exitting loop, free these components.
\ No newline at end of file
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/misc/gen_utils.C clean2/Particle/src/misc/gen_utils.C
--- dirty/Particle/src/misc/gen_utils.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/misc/gen_utils.C Wed May 1 17:51:41 2002
@@ -114,8 +114,10 @@ int generic_args (int argc, char *argv[]
fprintf(stderr, "%s: unable to open the mesh file -- %s or %s\n",
prog, *argv, fname);
exit(1);
}
+ if (f->mesh.name)
+ free(f->mesh.name);
f->mesh.name = strdup(fname);
break;
case 'r': /* read the session file */
@@ -130,16 +132,20 @@ int generic_args (int argc, char *argv[]
fprintf(stderr,"%s: unable to open the session file -- %s or %s\n",
prog, *argv, fname);
exit(1);
}
+ if (f->rea.name)
+ free(f->rea.name);
f->rea.name = strdup(fname);
/* Try to open the connectivity file. This file is usually *
* optional, so the application is responsible for error *
* checking. */
strcpy(fname + strlen(fname)-3, "mor");
f->mor.fp = fopen (fname,"r");
+ if (f->mor.name)
+ free(f->mor.name);
f->mor.name = strdup(fname);
break;
case 'o': /* direct the output to a file */
@@ -153,8 +159,10 @@ int generic_args (int argc, char *argv[]
fprintf(stderr, "%s: unable to open the output file -- %s\n",
prog, fname);
exit(1);
}
+ if (f->out.name)
+ free(f->out.name);
f->out.name = strdup(fname);
break;
case 'v': { /* be verbose and echo the version */
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/misc/io.C clean2/Particle/src/misc/io.C
--- dirty/Particle/src/misc/io.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/misc/io.C Wed May 1 17:51:41 2002
@@ -229,8 +229,11 @@ Element_List *ReadMesh (FILE *fp, char*
Coord X;
X.x = dvector(0,Max_Nverts-1);
X.y = dvector(0,Max_Nverts-1);
X.z = dvector(0,Max_Nverts-1);
+ dzero(Max_Nverts-1, X.x, 1);
+ dzero(Max_Nverts-1, X.y, 1);
+ dzero(Max_Nverts-1, X.z, 1);
/* Read in mesh information */
for(k = 0; k < nel; k++) {
fgets(buf, BUFSIZ, fp); /* element header */
@@ -277,8 +280,12 @@ Element_List *ReadMesh (FILE *fp, char*
fgets(buf_a, BUFSIZ, fp); /* finish the line */
}
}
+ free_dvector(X.x, 0);
+ free_dvector(X.y, 0);
+ free_dvector(X.z, 0);
+
for(k = 0; k < nel-1; ++k) new_E[k]->next = new_E[k+1];
new_E[k]->next = (Element*) NULL;
Element_List* E_List = (Element_List*) new Element_List(new_E, nel);
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/misc/prepost.C clean2/Particle/src/misc/prepost.C
--- dirty/Particle/src/misc/prepost.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/misc/prepost.C Wed May 1 17:51:41 2002
@@ -696,8 +696,10 @@ Gmap *GlobalNumScheme(Element_List *UL,
for(E = U; E; E = E->next)
for(i = 0; i < E->Nedges; ++i)
E->edge[i].gid = gbmap[E->edge[i].gid];
+ free(gsolve); free(gbmap);
+
if(U->dim() == 3){
/*------------------*/
/* Face re-ordering */
/*------------------*/
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/particle.C clean2/Particle/src/particle.C
--- dirty/Particle/src/particle.C Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/particle.C Wed May 1 17:51:41 2002
@@ -56,8 +56,36 @@ static int PG_set_nmb_of_RK_stages();
#define Timing(s) \
/* Nothing */
#endif
+FieldData::FieldData()
+ : steady(0), nload(0), nekdumps(0), dmpid(0), backwards(0),
+ dataskip(0), periodic(0), dumptime(0), fld_dump(0)
+{}
+
+FieldData::~FieldData()
+{
+ delete [] dataskip;
+ delete [] dumptime;
+
+ for (int i = 0; i < nload; ++i) {
+ freeField(fld_dump[i]);
+ }
+
+ // The field data ptrs are shifted around in Transient_Field. However, we
+ // know that they were created in a block by operator new [], so we just
+ // need to check their addresses to ascertain which one to pass to
+ // operator delete [].
+ Field * allocated_field = fld_dump[0];
+ for (int i = 0; i < nload; ++i) {
+ if (allocated_field > fld_dump[i])
+ allocated_field = fld_dump[i];
+ }
+ delete [] allocated_field;
+
+ free(fld_dump);
+}
+
/**
This is the main routine for particle tracker
*/
main (int argc, char *argv[]){
@@ -75,8 +103,10 @@ main (int argc, char *argv[]){
outfile = file.out.fp;
setup_field (&file,&U,&fld);
+ int const dim = fld.dim;
+
setup_particles (&file,U,&pclust,fld.dim);
P_summary (pclust,file.in.name);
setup_fdata (*U,&fdata,&fld,pclust);
@@ -124,17 +154,16 @@ main (int argc, char *argv[]){
}
for(i = 0; i < pclust->npts; ++i) pclust->bnd[i] = BndInit;
pclust->PC_output(step,outfile);
+ fclose(outfile);
// calculate error
double err;
for(i = 0; i < 3*pclust->npts; ++i)
err += (pclust->xyz[0].x[i]- pclust->xyzstart[0].x[i])*
(pclust->xyz[0].x[i]- pclust->xyzstart[0].x[i]);
fprintf(stderr,"Error: %lg\n",sqrt(err));
-
-
#endif
// -----------------------------------------------------------------
@@ -143,8 +172,27 @@ main (int argc, char *argv[]){
fprintf(stderr,"%d particle paths have been traced in %d step(s)\n",
pclust->npts,step);
+ // -----------------------------------------------------------------
+ for(int i = 0; i < dim; ++i) {
+ delete U[i];
+ }
+ free(U);
+
+ delete pclust;
+
+ if (file.in.name)
+ free(file.in.name);
+ if (file.rea.name)
+ free(file.rea.name);
+ if (file.mesh.name)
+ free(file.mesh.name);
+ if (file.mor.name)
+ free(file.mor.name);
+ if (file.out.name)
+ free(file.out.name);
+
return 0;
}
static int PG_set_nmb_of_RK_stages(){
@@ -259,10 +307,47 @@ static void Calc_data_skip(Element_List
return;
}
-Element_List *GMesh;
-Gmap *gmap;
+// We create this store for gmesh and gmap so that they can be free-d
+// when the program ends.
+// In turn, we do this to make it easier to debug genuine memory leaks.
+class GlobalVars {
+public:
+ // This is a singleton class. Get the instance.
+ static GlobalVars & get();
+
+ // Some compilers don't like the d-tor being private.
+ ~GlobalVars();
+
+ Element_List *gmesh;
+ Gmap *gmap;
+
+private:
+ // Make the c-tor private so we can control how many objects
+ // are instantiated.
+ GlobalVars();
+};
+
+GlobalVars & GlobalVars::get()
+{
+ static GlobalVars singleton;
+ return singleton;
+}
+
+GlobalVars::GlobalVars()
+ : gmesh(0), gmap(0)
+{}
+
+GlobalVars::~GlobalVars()
+{
+ delete gmesh;
+
+ free(gmap->vgid);
+ free(gmap->egid);
+ free(gmap->fgid);
+ free(gmap);
+}
static void setup_field(FileList *file, Element_List ***U, Field *fld){
int i;
int nfields,dump;
@@ -280,11 +365,12 @@ static void setup_field(FileList *file,
iparam_set("MODES", fld->lmax);
/* Generate the list of elements */
- GMesh = ReadMesh(file->rea.fp, strtok(file->rea.name,"."));
- gmap = GlobalNumScheme(GMesh, (Bndry *)NULL);
- U[0][0] = LocalMesh(GMesh,strtok(file->rea.name,"."));
+ GlobalVars & gvars = GlobalVars::get();
+ gvars.gmesh = ReadMesh(file->rea.fp, strtok(file->rea.name,"."));
+ gvars.gmap = GlobalNumScheme(gvars.gmesh, (Bndry *)NULL);
+ U[0][0] = LocalMesh(gvars.gmesh,strtok(file->rea.name,"."));
init_ortho_basis();
/* load field structures to element_list */
@@ -353,11 +439,11 @@ static void setup_particles(FileList *fi
break;
}
// setup mesh coordinates
- nel = U[0]->nel;
- pclust[0]->mesh = new Coord[nel];
- for(i = 0; i < nel; ++i){
+ pclust[0]->nel = U[0]->nel;
+ pclust[0]->mesh = new Coord[pclust[0]->nel];
+ for(i = 0; i < pclust[0]->nel; ++i){
qt = U[0]->flist[i]->qa*U[0]->flist[i]->qb;
if(dim == 3) qt *= U[0]->flist[i]->qc;
pclust[0]->mesh[i].x = dvector(0,dim*qt-1);
pclust[0]->mesh[i].y = pclust[0]->mesh[i].x + qt;
@@ -616,15 +702,14 @@ static void setup_fdata(Element_List *U,
}
for(i = 0; i < fld_dump->nel; ++i) // ensure field is at correct time
pclust->Prt_Interp_field_in_time(fdata,i,pclust->abstime[0]);
-
+
// set up global velocities
pclust->PC_set_glob_vel();
freeField(startfld);
}
-
}
static void P_summary (Particle_Cluster *pc, char *name ){
fprintf(stderr,"Input Information\n");
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/Particle/src/particle.h clean2/Particle/src/particle.h
--- dirty/Particle/src/particle.h Fri Feb 1 13:26:29 2002
+++ clean2/Particle/src/particle.h Wed May 1 17:51:41 2002
@@ -1,8 +1,11 @@
#include <veclib.h>
#include <hotel.h>
struct Butcher_arrays{
+ Butcher_arrays() : ai(0), bi(0), ci(0) {}
+ ~Butcher_arrays();
+
double **ai;
double **bi;
double *ci;
};
@@ -15,8 +18,11 @@ enum BndAct{
BndInit = -10
};
struct FieldData {
+ FieldData();
+ ~FieldData();
+
int steady; // data is steady if set to 1
int nload; // number of fields to load
int nekdumps; // number of nektar fields dumps available
int dmpid; // id of starting field currently stored in fld_dump
@@ -70,8 +76,9 @@ class Particle_Cluster {
int dim; /* dimensions of the problem */
int npts; /* number of particles */
int kopt; /* use optimised time step */
char status; /* scheme type 'H', 'P' or 'E' */
+ int nel; /* The number of elements in the mesh */
Coord *mesh; /* Mesh coordinates */
Coord *xyz; /* physical coords */
Coord *csi; /* cartesian local coords */
Coord *xyzstart; /* initial coordinates of particles */
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/include/nekstruct.h clean2/include/nekstruct.h
--- dirty/include/nekstruct.h Fri Feb 1 13:26:29 2002
+++ clean2/include/nekstruct.h Wed May 1 17:51:41 2002
@@ -391,8 +391,11 @@ class Element{
Element(); // default constructor
Element(const Element&); // copy constructor
Element(Element*); // copy constructor
+ // A base class without a virtual d-tor? You're mad and bad.
+ virtual ~Element();
+
// Local Matrix builders
virtual void MassMat (LocMat *); // return mass-matri(ces)
virtual void MassMatC(LocMat *);
virtual void HelmMatC(LocMat *, Metric *lambda); // return Helmholtz op.
@@ -602,8 +605,9 @@ public:
int nztot;
Element_List();
Element_List(Element **hea, int n);
+ ~Element_List();
virtual Element *operator()(int i);
virtual Element *operator[](int i);
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/include/smart_ptr.hpp clean2/include/smart_ptr.hpp
--- dirty/include/smart_ptr.hpp Thu Jan 1 01:00:00 1970
+++ clean2/include/smart_ptr.hpp Wed May 1 17:51:41 2002
@@ -0,0 +1,104 @@
+/**
+ * \file smart_ptr.hpp
+ * \author Angus Leeming, [EMAIL PROTECTED]
+ * \date April 2002
+ *
+ * class nektar::scoped_ptr is a small extension of the boost::scoped_ptr idea.
+ * See www.boost.org for more about boost.
+ *
+ * The original, boost smart_ptr.hpp file contains the following notice:
+ *
+ * (C) Copyright Greg Colvin and Beman Dawes 1998, 1999. Permission to copy,
+ * use, modify, sell and distribute this software is granted provided this
+ * copyright notice appears in all copies. This software is provided "as is"
+ * without express or implied warranty, and with no claim as to its
+ * suitability for any purpose.
+ *
+ * nektar::scoped_ptr can be used in identical fashion to the boost class,
+ * but can also be used to invoke an arbitrary functor when it goes out of
+ * scope.
+ *
+ * A simple example:
+ *
+ * void example()
+ * {
+ * nektar::scoped_c_ptr<Foo> scoped_foo((Foo*)malloc(sizeof(Foo)));
+ * Foo * foo = scoped_foo.get();
+ * ...
+ * // scoped_foo calls free(foo) automatically as it goes out of scope
+ * }
+ *
+ * The additional power lies in being able to define arbitrary functors to
+ * perform the destruction. Great with c-style structs that have no d-tor.
+ */
+
+#ifndef NEKTAR_SMART_PTR_HPP
+#define NEKTAR_SMART_PTR_HPP
+
+namespace nektar {
+
+class noncopyable {
+protected:
+ noncopyable(){}
+ ~noncopyable(){}
+private:
+ noncopyable( const noncopyable& );
+ const noncopyable& operator=( const noncopyable& );
+};
+
+template< typename T >
+struct DeletePtr {
+ void operator()(T * x) { delete x; }
+};
+
+template< typename T, typename Destroyer=DeletePtr<T> >
+class scoped_ptr : noncopyable {
+
+ T* ptr;
+
+ void destroy_ptr() {
+ if (ptr) {
+ Destroyer destroy;
+ destroy(ptr);
+ }
+ }
+
+public:
+ typedef T element_type;
+
+ explicit scoped_ptr( T* p=0 ) : ptr(p) {}
+ ~scoped_ptr() { destroy_ptr(); }
+ void reset( T* p=0 ) { if ( ptr != p ) { destroy_ptr(); ptr = p; } }
+ T& operator*() const { return *ptr; }
+ T* operator->() const { return ptr; }
+ T* get() const { return ptr; }
+};
+
+// and some shorthand helpers...
+template< typename T >
+struct DeleteArrayPtr {
+ void operator()(T * x) { delete [] x; }
+};
+
+
+template< typename T >
+struct FreePtr {
+ void operator()(T * x) { free(x); }
+};
+
+
+template< typename T >
+struct scoped_array_ptr : public scoped_ptr<T, DeleteArrayPtr<T> > {
+ explicit scoped_array_ptr( T* p=0 )
+ : scoped_ptr<T, DeleteArrayPtr<T> >(p) {}
+};
+
+template< typename T >
+struct scoped_c_ptr : public scoped_ptr<T, FreePtr<T> > {
+ explicit scoped_c_ptr( T* p=0 )
+ : scoped_ptr<T, FreePtr<T> >(p) {}
+};
+
+} // namespace nektar
+
+#endif // NEKTAR_SMART_PTR_HPP
========================================================================
diff -p -N -r -U 4 -X excl.tmp dirty/include/work.h clean2/include/work.h
--- dirty/include/work.h Fri Feb 1 13:26:29 2002
+++ clean2/include/work.h Wed May 1 17:51:41 2002
@@ -1,6 +1,8 @@
#ifndef NEKWORK
-extern double *Tri_wk;
+#include <smart_ptr.hpp>
+
+extern nektar::scoped_c_ptr<double> Tri_wk;
extern double *Quad_wk;
extern double *Quad_Jbwd_wk;
extern double *Quad_Iprod_wk;
