[cig-commits] r8945 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Dec 19 12:06:32 PST 2007
Author: luis
Date: 2007-12-19 12:06:32 -0800 (Wed, 19 Dec 2007)
New Revision: 8945
Modified:
cs/benchmark/cigma/trunk/src/Cell.cpp
cs/benchmark/cigma/trunk/src/Cell.h
cs/benchmark/cigma/trunk/src/Hex.cpp
cs/benchmark/cigma/trunk/src/Hex.h
cs/benchmark/cigma/trunk/src/Tet.cpp
cs/benchmark/cigma/trunk/src/Tet.h
Log:
Finally fixed the virtual table issue when calling interpolation routine!!
Modified: cs/benchmark/cigma/trunk/src/Cell.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Cell.cpp 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Cell.cpp 2007-12-19 20:06:32 UTC (rev 8945)
@@ -1,80 +1,132 @@
-#include "Cell.h"
+//#include <iostream>
#include <cstdlib>
#include <cassert>
+#include "Cell.h"
cigma::Cell::Cell()
{
- nsd = -1;
- celldim = -1;
+ //std::cout << "Calling cigma::Cell::Cell()\n";
- nno = 0;
+ /*nno = 0;
+ nsd = 0;
+ celldim = 0;*/
+
globverts = NULL;
refverts = NULL;
- ndofs = 0;
+ /*ndofs = 0;
basis_tab = NULL;
- basis_jet = NULL;
+ basis_jet = NULL;*/
- nq = 0;
- jxw = NULL;
+ /*nq = 0;
qpts = NULL;
qwts = NULL;
gqpts = NULL;
+ jxw = NULL;*/
- _tabulate = false;
+ //_tabulate = false;
}
cigma::Cell::~Cell()
{
-
+ //std::cout << "Calling cigma::Cell::~Cell()\n";
}
-//----------------------------------------------------------------------------
-
-void cigma::Cell::clear(void)
+/*
+int cigma::Cell::n_nodes() const
{
- globverts = NULL;
+ return nno;
}
+int cigma::Cell::n_dim() const
+{
+ return nsd;
+}
+
+int cigma::Cell::n_celldim() const
+{
+ return celldim;
+}*/
+
+//----------------------------------------------------------------------------
+
-void cigma::Cell::set(int nsd, int celldim, int ndofs)
+/*
+void cigma::Cell::set_dims(int ndofs, int celldim, int nsd)
{
+ this->nno = ndofs;
+ this->celldim = celldim;
this->nsd = nsd;
- this->celldim = celldim;
- this->ndofs = ndofs;
-}
+}*/
-void cigma::Cell::set_reference_vertices(const double *vertices, int num_vertices)
+void cigma::Cell::set_reference_vertices(double *vertices, int num_vertices)
{
- assert(nsd > 0);
- assert(celldim > 0);
+ int i,j;
- nno = num_vertices;
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ //const int nsd = n_dim();
+
+ assert(nno == num_vertices);
+
+ if (refverts != 0) delete [] refverts;
+ //if (globverts != 0) delete [] globverts;
+
refverts = new double[nno*celldim];
- for (int i = 0; i < nno; i++)
+ //globverts = new double[nno*nsd];
+
+ for (i = 0; i < nno; i++)
{
- for(int j = 0; j < celldim; j++)
+ /*
+ for (j = 0; j < nsd; j++)
{
+ globverts[nsd*i + j] = 0.0;
+ } // */
+ for (j = 0; j < celldim; j++)
+ {
int n = celldim*i + j;
refverts[n] = vertices[n];
+ //globverts[nsd*i+j] = vertices[n];
}
}
globverts = refverts;
}
-void cigma::Cell::update_vertices(double *vertices, int num_vertices, int nsd)
+void cigma::Cell::update_vertices(double *vertices, int num, int dim)
{
- assert((this->nsd) == nsd);
- assert(nno == num_vertices);
- this->globverts = vertices;
+ assert(num == n_nodes());
+ assert(dim == n_dim());
+
+ //* // Copy reference
+ globverts = vertices;
+ // */
+
+ /* // Copy data
+ int i,j;
+ for (i = 0; i < nno; i++)
+ {
+ // zero out globverts
+ for (j = 0; j < nsd; j++)
+ {
+ globverts[nsd*i+j] = 0.0;
+ }
+ // copy globverts
+ for (j = 0; j < dim; j++)
+ {
+ globverts[nsd*i + j] = vertices[dim*i + j];
+ }
+ } // */
+
+ return;
}
+/*
void cigma::Cell::set_quadrature(const double *quadpts, const double *quadwts, int num_points)
{
assert(nsd > 0);
@@ -96,8 +148,6 @@
update_quadrature();
}
}
-
-
void cigma::Cell::update_quadrature(void)
{
assert(nsd > 0);
@@ -105,6 +155,7 @@
for (int q = 0; q < nq; q++)
uvw2xyz(&qpts[celldim*q], &gqpts[nsd*q]);
}
+*/
//----------------------------------------------------------------------------
@@ -112,7 +163,7 @@
double cigma::Cell::jacobian(double *point, double jac[3][3])
{
- assert(celldim > 0);
+ const int celldim = n_celldim();
switch (celldim)
{
case 3: return jacobian(point[0], point[1], point[2], jac);
@@ -129,8 +180,11 @@
jac[1][0] = jac[1][1] = jac[1][2] = 0.0;
jac[2][0] = jac[2][1] = jac[2][2] = 0.0;
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+
double uvw[3] = {u,v,w};
- double *grad = new double[ndofs*3];
+ double *grad = new double[nno*3];
double *s;
#define X(i) globverts[3*(i) + 0]
@@ -265,7 +319,8 @@
*/
void cigma::Cell::interpolate(double *dofs, double *point, double *value, int valdim)
{
- double N[ndofs];
+ const int nno = n_nodes();
+ double N[nno];
shape(1, point, N);
@@ -275,7 +330,7 @@
for (int i = 0; i < valdim; i++)
{
double sum = 0.0;
- for (int j = 0; j < ndofs; j++)
+ for (int j = 0; j < nno; j++)
sum += dofs[i + valdim*j] * N[j];
value[i] = sum;
}
@@ -289,20 +344,23 @@
*/
void cigma::Cell::interpolate_grad(double *dofs, double *point, double *value, int stride, double invjac[3][3])
{
- assert(celldim > 0);
- assert(nno > 0);
- assert(ndofs > 0);
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+
+ //assert(nno > 0);
+ //assert(celldim > 0);
+ //assert(ndofs > 0); <-- ndofs is (cell->nno)
int i, j, k;
double dfdu[3] = {0.0, 0.0, 0.0};
- double grad[ndofs*3];
+ double grad[nno*3];
double *s;
grad_shape(1, point, grad);
k = 0;
- for (i = 0; i < ndofs; i++)
+ for (i = 0; i < nno; i++)
{
s = &grad[3*i];
for (j = 0; j < celldim; j++)
@@ -334,6 +392,7 @@
// general newton routine for the nonlinear case
uvw[0] = uvw[1] = uvw[2] = 0.0;
+ const int nno = n_nodes();
int iter = 1, maxiter = 20;
double error = 1.0, tol = 1.0e-6;
@@ -346,7 +405,7 @@
double xn = 0.0, yn = 0.0, zn = 0.0;
- double s[ndofs];
+ double s[nno];
shape(1, uvw, s);
@@ -396,18 +455,23 @@
/* isoparametric reference map */
void cigma::Cell::uvw2xyz(double uvw[3], double xyz[3])
{
+ assert(globverts != 0);
interpolate(globverts, uvw, xyz, 3);
}
void cigma::Cell::bbox(double *min, double *max)
{
+ const int nno = n_nodes();
+ const int nsd = n_dim();
minmax(globverts, nno, nsd, min, max);
}
void cigma::Cell::centroid(double c[3])
{
+ const int nno = n_nodes();
+ const int nsd = n_dim();
cigma::centroid(globverts, nno, nsd, c);
}
Modified: cs/benchmark/cigma/trunk/src/Cell.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Cell.h 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Cell.h 2007-12-19 20:06:32 UTC (rev 8945)
@@ -21,25 +21,25 @@
public:
- int n_dims() const;
- int n_celldims() const;
- int n_nodes() const;
+ virtual int n_nodes() = 0;
+ virtual int n_celldim() = 0;
+ virtual int n_dim() = 0;
-
public:
- void set(int nsd, int celldim, int ndofs);
+ //void set_dims(int ndofs, int celldim, int nsd);
- void set_reference_vertices(const double *vertices, int num_vertices);
+ void set_reference_vertices(double *vertices, int num_vertices);
+
void update_vertices(double *vertices, int num_vertices, int nsd);
- void set_quadrature(const double *quadpts, const double *quadwts, int num_points);
- void update_quadrature(void);
+ //void set_quadrature(const double *quadpts, const double *quadwts, int num_points);
+ //void update_quadrature(void);
- void set_tabulation(const double *basis_tab, const double *basis_jet);
- void update_tabulation(void);
+ //void set_tabulation(const double *basis_tab, const double *basis_jet);
+ //void update_tabulation(void);
- void clear();
+ //void clear();
public:
@@ -47,8 +47,6 @@
virtual void shape(int num, double *points, double *values) = 0;
virtual void grad_shape(int num, double *points, double *values) = 0;
- void tabulate(void);
-
double jacobian(double *point, double jac[3][3]);
double jacobian(double u, double v, double w, double jac[3][3]);
@@ -68,45 +66,31 @@
public:
- int nsd;
+ /*int nno;
int celldim;
+ int nsd;*/
- int nno;
double *refverts; // [nno x celldim]
double *globverts; // [nno x nsd]
- int nq;
- double *jxw; // [nq x 1]
- double *qwts; // [nq x 1]
- double *qpts; // [nq x celldim]
- double *gqpts; // [nq x nsd]
+ //int nq;
+ //double *jxw; // [nq x 1]
+ //double *qwts; // [nq x 1]
+ //double *qpts; // [nq x celldim]
+ //double *gqpts; // [nq x nsd]
- int ndofs;
- double *basis_tab; // [nq x ndofs]
- double *basis_jet; // [nq x ndofs x celldim]
+ //int ndofs;
+ //double *basis_tab; // [nq x ndofs]
+ //double *basis_jet; // [nq x ndofs x celldim]
- bool _tabulate;
+ //bool _tabulate;
};
-/*
- * Inline methods
- */
-inline int cigma::Cell::n_dims() const
-{
- return nsd;
-}
+// ---------------------------------------------------------------------------
-inline int cigma::Cell::n_celldims() const
-{
- return celldim;
-}
-inline int cigma::Cell::n_nodes() const
-{
- return nno;
-}
+// ---------------------------------------------------------------------------
-
#endif
Modified: cs/benchmark/cigma/trunk/src/Hex.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.cpp 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Hex.cpp 2007-12-19 20:06:32 UTC (rev 8945)
@@ -1,6 +1,8 @@
#include "Hex.h"
#include <cassert>
+// ---------------------------------------------------------------------------
+
static void hex_shape(double u, double v, double w, double s[8])
{
s[0] = (1.0 - u) * (1.0 - v) * (1.0 - w) * 0.125;
@@ -52,12 +54,15 @@
#undef S
}
+
+// ---------------------------------------------------------------------------
+
cigma::Hex::Hex()
{
- const int nno = 8;
- const int nsd = 3;
- const int celldim = 3;
- double verts[nno * celldim] = {
+ const int hex_nno = 8;
+ const int hex_celldim = 3;
+ //const int hex_nsd = 3;
+ double verts[hex_nno * hex_celldim] = {
-1.0, -1.0, -1.0,
+1.0, -1.0, -1.0,
+1.0, +1.0, -1.0,
@@ -67,9 +72,8 @@
+1.0, +1.0, +1.0,
-1.0, +1.0, +1.0
};
-
- set(nsd, celldim, nno);
- set_reference_vertices(verts, nno);
+ //set_dims(hex_nno, hex_celldim, hex_nsd);
+ set_reference_vertices(verts, hex_nno);
}
cigma::Hex::~Hex()
@@ -84,13 +88,13 @@
*/
void cigma::Hex::shape(int num, double *points, double *values)
{
- assert(ndofs > 0);
+ const int nno = n_nodes();
for (int i = 0; i < num; i++)
{
double u = points[3*i + 0];
double v = points[3*i + 1];
double w = points[3*i + 2];
- hex_shape(u, v, w, &values[ndofs*i]);
+ hex_shape(u, v, w, &values[nno*i]);
}
}
@@ -101,13 +105,15 @@
*/
void cigma::Hex::grad_shape(int num, double *points, double *values)
{
- assert(ndofs > 0);
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ const int stride = nno * celldim;
for (int i = 0; i < num; i++)
{
double u = points[3*i + 0];
double v = points[3*i + 1];
double w = points[3*i + 2];
- hex_grad_shape(u, v, w, &values[ndofs*celldim*i]);
+ hex_grad_shape(u, v, w, &values[stride*i]);
}
}
Modified: cs/benchmark/cigma/trunk/src/Hex.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.h 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Hex.h 2007-12-19 20:06:32 UTC (rev 8945)
@@ -15,6 +15,11 @@
~Hex();
public:
+ int n_nodes() { return 8; }
+ int n_celldim() { return 3; }
+ int n_dim() { return 3; }
+
+public:
void shape(int num, double *points, double *values);
void grad_shape(int num, double *points, double *values);
bool interior(double u, double v, double w);
Modified: cs/benchmark/cigma/trunk/src/Tet.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.cpp 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Tet.cpp 2007-12-19 20:06:32 UTC (rev 8945)
@@ -1,6 +1,9 @@
+//#include <iostream>
+#include <cassert>
#include "Tet.h"
-#include <cassert>
+// ---------------------------------------------------------------------------
+
static void tet_shape(double u, double v, double w, double s[4])
{
s[0] = 1.0 - u - v - w;
@@ -28,46 +31,53 @@
s[3*3+2] = +1.0;
}
+
+// ---------------------------------------------------------------------------
+
cigma::Tet::Tet()
{
- const int nno = 4;
- const int nsd = 3;
- const int celldim = 3;
- double verts[nno * celldim] = {
+ //std::cout << "Calling cigma::Tet::Tet()\n";
+ const int tet_nno = 4;
+ const int tet_celldim = 3;
+ //const int tet_nsd = 3;
+ double verts[tet_nno * tet_celldim] = {
0.0, 0.0, 0.0,
1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0
};
- set(nsd, celldim, nno);
- set_reference_vertices(verts, nno);
+ //set_dims(tet_nno, tet_celldim, tet_nsd);
+ set_reference_vertices(verts, tet_nno);
}
cigma::Tet::~Tet()
{
+ //std::cout << "Calling cigma::Tet::~Tet()\n";
}
void cigma::Tet::shape(int num, double *points, double *values)
{
- assert(ndofs > 0);
+ const int nno = n_nodes();
for (int i = 0; i < num; i++)
{
double u = points[3*i + 0];
double v = points[3*i + 1];
double w = points[3*i + 2];
- tet_shape(u, v, w, &values[ndofs*i]);
+ tet_shape(u, v, w, &values[nno*i]);
}
}
void cigma::Tet::grad_shape(int num, double *points, double *values)
{
- assert(ndofs > 0);
+ const int nno = n_nodes();
+ const int celldim = n_celldim();
+ const int stride = nno * celldim;
for (int i = 0; i < num; i++)
{
double u = points[3*i + 0];
double v = points[3*i + 1];
double w = points[3*i + 2];
- tet_grad_shape(u, v, w, &values[ndofs*celldim*i]);
+ tet_grad_shape(u, v, w, &values[stride*i]);
}
}
Modified: cs/benchmark/cigma/trunk/src/Tet.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.h 2007-12-19 20:06:25 UTC (rev 8944)
+++ cs/benchmark/cigma/trunk/src/Tet.h 2007-12-19 20:06:32 UTC (rev 8945)
@@ -15,6 +15,11 @@
~Tet();
public:
+ int n_nodes() { return 4; }
+ int n_celldim() { return 3; }
+ int n_dim() { return 3; }
+
+public:
void shape(int num, double *points, double *values);
void grad_shape(int num, double *points, double *values);
void xyz2uvw(double xyz[3], double uvw[3]);
More information about the cig-commits
mailing list