[cig-commits] r11712 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Apr 2 03:24:41 PDT 2008
Author: luis
Date: 2008-04-02 03:24:41 -0700 (Wed, 02 Apr 2008)
New Revision: 11712
Modified:
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/Quad.cpp
cs/benchmark/cigma/trunk/src/Quad.h
cs/benchmark/cigma/trunk/src/QuadraturePoints.cpp
cs/benchmark/cigma/trunk/src/QuadraturePoints.h
cs/benchmark/cigma/trunk/src/QuadratureReader.cpp
cs/benchmark/cigma/trunk/src/Tet.cpp
cs/benchmark/cigma/trunk/src/Tet.h
cs/benchmark/cigma/trunk/src/Tri.cpp
cs/benchmark/cigma/trunk/src/Tri.h
Log:
Move default quadrature rules into its respective cell subclass
Modified: cs/benchmark/cigma/trunk/src/Cell.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Cell.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Cell.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -28,6 +28,9 @@
void set_global_vertices(double *vertices, int num_vertices, int nsd);
public:
+ virtual void qr_default(double **wts, double **pts, int *npts, int *ndim) = 0;
+
+public:
virtual void shape(int num, double *points, double *values) = 0;
virtual void grad_shape(int num, double *points, double *values) = 0;
Modified: cs/benchmark/cigma/trunk/src/Hex.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Hex.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -82,6 +82,42 @@
}
+void Hex::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // hex_qr(3)
+ const int hex_nno = 8;
+ const int hex_celldim = 3;
+ double hex_qpts[hex_nno * hex_celldim] = {
+ -0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, -0.57735027, -0.57735027,
+ 0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, 0.57735027, -0.57735027,
+ -0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, -0.57735027, 0.57735027,
+ 0.57735027, 0.57735027, 0.57735027,
+ -0.57735027, 0.57735027, 0.57735027
+ };
+ double hex_qwts[hex_nno] = { 1., 1., 1., 1., 1., 1., 1., 1. };
+
+ int i,j,n;
+ *npts = hex_nno;
+ *ndim = hex_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < hex_nno; i++)
+ {
+ (*wts)[i] = hex_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = hex_qpts[n];
+ }
+ }
+}
+
+
+// ---------------------------------------------------------------------------
+
/*
* @param points points is an [num x celldim] array (in)
* @param values values is an [num x ndofs] array (out)
Modified: cs/benchmark/cigma/trunk/src/Hex.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Hex.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Hex.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -13,6 +13,7 @@
public:
Hex();
~Hex();
+ void qr_default(double **wts, double **pts, int *npts, int *ndim);
public:
int n_nodes() { return 8; }
Modified: cs/benchmark/cigma/trunk/src/Quad.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Quad.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Quad.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -50,7 +50,52 @@
{
}
+void Quad::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // quad_qr(7)
+ const int quad_nno = 16;
+ const int quad_celldim = 2;
+ double quad_qpts[quad_nno * quad_celldim] = {
+ -0.86113631, -0.86113631,
+ -0.33998104, -0.86113631,
+ 0.33998104, -0.86113631,
+ 0.86113631, -0.86113631,
+ 0.86113631, -0.33998104,
+ 0.86113631, 0.33998104,
+ 0.86113631, 0.86113631,
+ 0.33998104, 0.86113631,
+ -0.33998104, 0.86113631,
+ -0.86113631, 0.86113631,
+ -0.86113631, 0.33998104,
+ -0.86113631, -0.33998104,
+ -0.33998104, -0.33998104,
+ 0.33998104, -0.33998104,
+ -0.33998104, 0.33998104,
+ 0.33998104, 0.33998104
+ };
+ double quad_qwts[quad_nno] = {
+ 0.12100299, 0.22685185, 0.22685185, 0.12100299, 0.22685185,
+ 0.22685185, 0.12100299, 0.22685185, 0.22685185, 0.12100299,
+ 0.22685185, 0.22685185, 0.4252933 , 0.4252933 , 0.4252933 ,
+ 0.4252933
+ };
+ int i,j,n;
+ *npts = quad_nno;
+ *ndim = quad_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < quad_nno; i++)
+ {
+ (*wts)[i] = quad_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = quad_qpts[n];
+ }
+ }
+}
+
// ---------------------------------------------------------------------------
/*
Modified: cs/benchmark/cigma/trunk/src/Quad.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Quad.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Quad.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -13,6 +13,7 @@
public:
Quad();
~Quad();
+ void qr_default(double **wts, double **pts, int *npts, int *ndim);
public:
int n_nodes() { return 4; }
Modified: cs/benchmark/cigma/trunk/src/QuadraturePoints.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadraturePoints.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/QuadraturePoints.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -38,7 +38,7 @@
// ---------------------------------------------------------------------------
-void QuadraturePoints::set_quadrature(double *qpts, double *qwts, int npts, int qdim)
+void QuadraturePoints::set_quadrature(double *qwts, double *qpts, int npts, int qdim)
{
/* some basic assertions */
assert(qpts != 0);
Modified: cs/benchmark/cigma/trunk/src/QuadraturePoints.h
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadraturePoints.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/QuadraturePoints.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -22,7 +22,7 @@
void clear();
public:
- void set_quadrature(double *qpts, double *qwts, int npts, int qdim);
+ void set_quadrature(double *qwts, double *qpts, int npts, int qdim);
void apply_refmap(Cell *cell);
public:
Modified: cs/benchmark/cigma/trunk/src/QuadratureReader.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadratureReader.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/QuadratureReader.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -111,124 +111,16 @@
void QuadratureReader::load_quadrature(Cell *cell)
{
- //cout << "Calling QuadratureReader::load_quadrature()" << endl;
-
assert(cell != 0);
-
int i,j;
- //
- // XXX: move these default rules into the appropriate Cell subclasses
- //
-
- // tri_qr(5)
- const int tri_nno = 9;
- const int tri_celldim = 2;
- double tri_qpts[tri_nno * tri_celldim] = {
- -0.79456469, -0.82282408,
- -0.86689186, -0.18106627,
- -0.95213774, 0.57531892,
- -0.08858796, -0.82282408,
- -0.40946686, -0.18106627,
- -0.78765946, 0.57531892,
- 0.61738877, -0.82282408,
- 0.04795814, -0.18106627,
- -0.62318119, 0.57531892
- };
- double tri_qwts[tri_nno] = {
- 0.22325768, 0.25471234, 0.07758553, 0.35721229, 0.40753974,
- 0.12413685, 0.22325768, 0.25471234, 0.07758553
- };
- for (i = 0; i < tri_nno; i++)
- {
- for (j = 0; j < tri_celldim; j++)
- {
- int n = tri_celldim * i + j;
- tri_qpts[n] += 1;
- tri_qpts[n] *= 0.5;
- }
- }
-
- // quad_qr(7)
- const int quad_nno = 16;
- const int quad_celldim = 2;
- double quad_qpts[quad_nno * quad_celldim] = {
- -0.86113631, -0.86113631,
- -0.33998104, -0.86113631,
- 0.33998104, -0.86113631,
- 0.86113631, -0.86113631,
- 0.86113631, -0.33998104,
- 0.86113631, 0.33998104,
- 0.86113631, 0.86113631,
- 0.33998104, 0.86113631,
- -0.33998104, 0.86113631,
- -0.86113631, 0.86113631,
- -0.86113631, 0.33998104,
- -0.86113631, -0.33998104,
- -0.33998104, -0.33998104,
- 0.33998104, -0.33998104,
- -0.33998104, 0.33998104,
- 0.33998104, 0.33998104
- };
- double quad_qwts[quad_nno] = {
- 0.12100299, 0.22685185, 0.22685185, 0.12100299, 0.22685185,
- 0.22685185, 0.12100299, 0.22685185, 0.22685185, 0.12100299,
- 0.22685185, 0.22685185, 0.4252933 , 0.4252933 , 0.4252933 ,
- 0.4252933
- };
-
-
- // tet_qr(3)
- const int tet_nno = 8;
- const int tet_celldim = 3;
- double tet_qpts[tet_nno * tet_celldim] = {
- -0.68663473, -0.72789005, -0.75497035,
- -0.83720867, -0.85864055, 0.08830369,
- -0.86832263, 0.13186633, -0.75497035,
- -0.93159441, -0.4120024 , 0.08830369,
- 0.16949513, -0.72789005, -0.75497035,
- -0.39245447, -0.85864055, 0.08830369,
- -0.50857335, 0.13186633, -0.75497035,
- -0.74470688, -0.4120024 , 0.08830369
- };
- double tet_qwts[tet_nno] = {
- 0.29583885, 0.12821632, 0.16925605, 0.07335544, 0.29583885,
- 0.12821632, 0.16925605, 0.07335544
- };
- for (i = 0; i < tet_nno; i++)
- {
- for (j = 0; j < tet_celldim; j++)
- {
- int n = tet_celldim * i + j;
- tet_qpts[n] += 1;
- tet_qpts[n] *= 0.5;
- }
- }
-
- // hex_qr(3)
- const int hex_nno = 8;
- const int hex_celldim = 3;
- double hex_qpts[hex_nno * hex_celldim] = {
- -0.57735027, -0.57735027, -0.57735027,
- 0.57735027, -0.57735027, -0.57735027,
- 0.57735027, 0.57735027, -0.57735027,
- -0.57735027, 0.57735027, -0.57735027,
- -0.57735027, -0.57735027, 0.57735027,
- 0.57735027, -0.57735027, 0.57735027,
- 0.57735027, 0.57735027, 0.57735027,
- -0.57735027, 0.57735027, 0.57735027
- };
- double hex_qwts[hex_nno] = { 1., 1., 1., 1., 1., 1., 1., 1. };
-
-
int nq, nd;
double *qx, *qw;
nq = nd = 0;
qx = qw = 0;
-
if ((pointsPath != "") && (weightsPath != ""))
{
int ierr;
@@ -348,7 +240,7 @@
int npts,dim;
double *qx,*qw;
fiat->quadrature(geometry, order, &qx, &qw, &npts, &dim);
- quadrature->set_quadrature(qx, qw, nno, nsd);
+ quadrature->set_quadrature(qw, qx, nno, nsd);
delete [] qx;
delete [] qw;
// */
@@ -357,33 +249,15 @@
else
{
// assign reasonable defaults
- switch (cell->geometry())
- {
- case Cell::TRIANGLE:
- quadrature = new QuadraturePoints();
- quadrature->set_quadrature(tri_qpts, tri_qwts, tri_nno, tri_celldim);
- break;
- case Cell::QUADRANGLE:
- quadrature = new QuadraturePoints();
- quadrature->set_quadrature(quad_qpts, quad_qwts, quad_nno, quad_celldim);
- break;
- case Cell::TETRAHEDRON:
- quadrature = new QuadraturePoints();
- quadrature->set_quadrature(tet_qpts, tet_qwts, tet_nno, tet_celldim);
- break;
- case Cell::HEXAHEDRON:
- quadrature = new QuadraturePoints();
- quadrature->set_quadrature(hex_qpts, hex_qwts, hex_nno, hex_celldim);
- break;
- default:
- break;
- }
+ quadrature = new QuadraturePoints();
+ cell->qr_default(&qw, &qx, &nq, &nd);
+ quadrature->set_quadrature(qw, qx, nq, nd);
}
if ((qx != 0) && (qw != 0))
{
quadrature = new QuadraturePoints();
- quadrature->set_quadrature(qx, qw, nq, nd);
+ quadrature->set_quadrature(qw, qx, nq, nd);
}
assert(quadrature != 0);
Modified: cs/benchmark/cigma/trunk/src/Tet.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Tet.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -58,6 +58,53 @@
{
}
+void Tet::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // tet_qr(3)
+ const int tet_nno = 8;
+ const int tet_celldim = 3;
+ double tet_qpts[tet_nno * tet_celldim] = {
+ -0.68663473, -0.72789005, -0.75497035,
+ -0.83720867, -0.85864055, 0.08830369,
+ -0.86832263, 0.13186633, -0.75497035,
+ -0.93159441, -0.4120024 , 0.08830369,
+ 0.16949513, -0.72789005, -0.75497035,
+ -0.39245447, -0.85864055, 0.08830369,
+ -0.50857335, 0.13186633, -0.75497035,
+ -0.74470688, -0.4120024 , 0.08830369
+ };
+ double tet_qwts[tet_nno] = {
+ 0.29583885, 0.12821632, 0.16925605, 0.07335544, 0.29583885,
+ 0.12821632, 0.16925605, 0.07335544
+ };
+
+ int i,j,n;
+
+ for (i = 0; i < tet_nno; i++)
+ {
+ for (j = 0; j < tet_celldim; j++)
+ {
+ n = tet_celldim * i + j;
+ tet_qpts[n] += 1;
+ tet_qpts[n] *= 0.5;
+ }
+ }
+
+ *npts = tet_nno;
+ *ndim = tet_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < tet_nno; i++)
+ {
+ (*wts)[i] = tet_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = tet_qpts[n];
+ }
+ }
+}
+
void Tet::shape(int num, double *points, double *values)
{
const int nno = n_nodes();
Modified: cs/benchmark/cigma/trunk/src/Tet.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tet.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Tet.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -13,6 +13,7 @@
public:
Tet();
~Tet();
+ void qr_default(double **wts, double **pts, int *npts, int *ndim);
public:
int n_nodes() { return 4; }
Modified: cs/benchmark/cigma/trunk/src/Tri.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.cpp 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Tri.cpp 2008-04-02 10:24:41 UTC (rev 11712)
@@ -42,10 +42,59 @@
set_reference_vertices(verts, tri_nno);
}
-cigma::Tri::~Tri()
+Tri::~Tri()
{
}
+void Tri::qr_default(double **wts, double **pts, int *npts, int *ndim)
+{
+ // tri_qr(5)
+ const int tri_nno = 9;
+ const int tri_celldim = 2;
+ double tri_qpts[tri_nno * tri_celldim] = {
+ -0.79456469, -0.82282408,
+ -0.86689186, -0.18106627,
+ -0.95213774, 0.57531892,
+ -0.08858796, -0.82282408,
+ -0.40946686, -0.18106627,
+ -0.78765946, 0.57531892,
+ 0.61738877, -0.82282408,
+ 0.04795814, -0.18106627,
+ -0.62318119, 0.57531892
+ };
+ double tri_qwts[tri_nno] = {
+ 0.22325768, 0.25471234, 0.07758553, 0.35721229, 0.40753974,
+ 0.12413685, 0.22325768, 0.25471234, 0.07758553
+ };
+
+ int i,j,n;
+
+ for (i = 0; i < tri_nno; i++)
+ {
+ for (j = 0; j < tri_celldim; j++)
+ {
+ n = tri_celldim * i + j;
+ tri_qpts[n] += 1;
+ tri_qpts[n] *= 0.5;
+ }
+ }
+
+ *npts = tri_nno;
+ *ndim = tri_celldim;
+ *wts = new double[(*npts)];
+ *pts = new double[(*npts) * (*ndim)];
+ for (i = 0; i < tri_nno; i++)
+ {
+ (*wts)[i] = tri_qwts[i];
+ for (j = 0; j < (*ndim); j++)
+ {
+ n = (*ndim) * i + j;
+ (*pts)[n] = tri_qpts[n];
+ }
+ }
+}
+
+
// ---------------------------------------------------------------------------
void Tri::shape(int num, double *points, double *values)
Modified: cs/benchmark/cigma/trunk/src/Tri.h
===================================================================
--- cs/benchmark/cigma/trunk/src/Tri.h 2008-04-01 20:35:59 UTC (rev 11711)
+++ cs/benchmark/cigma/trunk/src/Tri.h 2008-04-02 10:24:41 UTC (rev 11712)
@@ -13,6 +13,7 @@
public:
Tri();
~Tri();
+ void qr_default(double **wts, double **pts, int *npts, int *ndim);
public:
int n_nodes() { return 3; }
More information about the cig-commits
mailing list