[cig-commits] r11189 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Tue Feb 19 20:30:42 PST 2008
Author: luis
Date: 2008-02-19 20:30:42 -0800 (Tue, 19 Feb 2008)
New Revision: 11189
Modified:
cs/benchmark/cigma/trunk/src/CompareCmd.cpp
cs/benchmark/cigma/trunk/src/FE_Field.cpp
cs/benchmark/cigma/trunk/src/FE_Field.h
Log:
Improvements to tabulation of shape functions (65% faster comparisons)
Modified: cs/benchmark/cigma/trunk/src/CompareCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 04:30:41 UTC (rev 11188)
+++ cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 04:30:42 UTC (rev 11189)
@@ -264,36 +264,54 @@
// indices
int e,q;
+ int i,j;
// dimensions
int nel = mesh->nel;
int nq = quadrature->n_points();
//int celldim = cell_a->n_celldim();
+ int ndofs = cell_a->n_nodes();
int valdim = field_a->n_rank();
// local data;
//double *jxw = new double[nq];
+ double *jxw = field_a->fe->jxw;
+ double *dofs_a = new double[ndofs * valdim];
double *phi_a = new double[nq * valdim];
double *phi_b = new double[nq * valdim];
- double *jxw = field_a->fe->jxw;
// norm
double L2 = 0.0;
double *epsilon = new double[nel];
+ // XXX: what about case (mesh != field_a->meshPart)?
//FE *fe;
Cell *cell = cell_a;
+ // start timer
Timer timer;
if (verbose)
{
- std::cout << std::setprecision(4);
+ std::cout << std::setprecision(5);
timer.print_header(std::cout, "elts");
timer.start(nel);
timer.update(0);
std::cout << timer;
}
+
+ if (true)
+ {
+ // XXX: move this block into FE::tabulate(...)
+
+ // get shape function values at known quadrature points
+ cell->shape(nq, quadrature->qpts, field_a->fe->basis_tab);
+
+ // get shape function derivatives at known quadrature points
+ cell->grad_shape(nq, quadrature->qpts, field_a->fe->basis_jet);
+ }
+
+
for (e = 0; e < nel; e++)
{
// update cell data
@@ -304,13 +322,31 @@
// obtain global points from current quadrature rule
quadrature->apply_refmap(cell);
- field_a->tabulate();
// ... calculate phi_a[]
- // XXX: use tabulation instead of calling eval repeatedly!
+ // XXX: using eval()
+ //for (q = 0; q < nq; q++)
+ //{
+ // field_a->eval((*quadrature)[q], &phi_a[valdim*q]);
+ //}
+
+ // ... calculate phi_a[]
+ // XXX: using tabulation
+
+ field_a->get_cell_dofs(e, dofs_a);
+ //field_a->tabulate();
for (q = 0; q < nq; q++)
{
- field_a->eval((*quadrature)[q], &phi_a[valdim*q]);
+ double *N = &(field_a->fe->basis_tab[ndofs*q]);
+ for (i = 0; i < valdim; i++)
+ {
+ double valsum = 0.0;
+ for (j = 0; j < ndofs; j++)
+ {
+ valsum += dofs_a[i + valdim*j] * N[j];
+ }
+ phi_a[valdim*q + i] = valsum;
+ }
}
// ... calculate phi_b[]
@@ -319,11 +355,14 @@
field_b->eval((*quadrature)[q], &phi_b[valdim*q]);
}
+ // evaluate jacobian at known quadrature points
+ field_a->fe->update_jxw();
+
// ... apply quadrature rule
double err = 0.0;
for (q = 0; q < nq; q++)
{
- for (int i = 0; i < valdim; i++)
+ for (i = 0; i < valdim; i++)
{
int qi = valdim*q + i;
err += jxw[q] * SQR(phi_a[qi] - phi_b[qi]);
@@ -408,3 +447,6 @@
return 0;
}
+
+
+// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/FE_Field.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 04:30:41 UTC (rev 11188)
+++ cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 04:30:42 UTC (rev 11189)
@@ -73,7 +73,21 @@
void cigma::FE_Field::tabulate()
{
+ // XXX: move this function to FE::tabulate()
+
assert(fe != 0);
+
+ // get shape function values at known quadrature points
+ fe->cell->shape(fe->quadrature->n_points(), fe->quadrature->qpts, fe->basis_tab);
+
+ // get shape function derivatives at known quadrature points
+ fe->cell->grad_shape(fe->quadrature->n_points(), fe->quadrature->qpts, fe->basis_jet);
+}
+
+
+void cigma::FE_Field::eval(int e, double *values)
+{
+ assert(fe != 0);
assert(meshPart != 0);
Cell *cell = fe->cell;
@@ -84,16 +98,39 @@
assert(quadrature != 0);
int nq = quadrature->n_points();
double *qpts = quadrature->qpts;
- double *qwts = quadrature->qwts;
+ //double *qwts = quadrature->qwts;
// get shape function values at known quadrature points
- cell->shape(nq, qpts, fe->basis_tab);
+ //cell->shape(nq, qpts, fe->basis_tab);
// get shape function derivatives at known quadrature points
- cell->grad_shape(nq, qpts, fe->basis_jet);
+ //cell->grad_shape(nq, qpts, fe->basis_jet);
+ // XXX: move this step out of tabulate()
// evaluate jacobian at known quadrature points and calculate jxw
- fe->update_jxw();
+ //fe->update_jxw();
+
+ // tabulate the function values
+ int i,j;
+ const int valdim = n_rank();
+ const int ndofs = cell->n_nodes();
+ double dofs[ndofs * valdim];
+
+ get_cell_dofs(e, dofs);
+
+ for (int q = 0; q < nq; q++)
+ {
+ double *N = &(fe->basis_tab[ndofs*q]);
+ for (i = 0; i < valdim; i++)
+ {
+ double sum = 0.0;
+ for (j = 0; j < ndofs; j++)
+ {
+ sum += dofs[i + valdim*j] * N[j];
+ }
+ values[valdim*q + i] = sum;
+ }
+ }
}
// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/FE_Field.h
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 04:30:41 UTC (rev 11188)
+++ cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 04:30:42 UTC (rev 11189)
@@ -29,7 +29,10 @@
public:
void eval(double *point, double *value);
+
+public:
void tabulate();
+ void eval(int e, double *values);
public:
void get_cell_dofs(int cellIndex, double *cellDofs);
More information about the cig-commits
mailing list