[cig-commits] r11191 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 20 01:27:52 PST 2008
Author: luis
Date: 2008-02-20 01:27:51 -0800 (Wed, 20 Feb 2008)
New Revision: 11191
Modified:
cs/benchmark/cigma/trunk/src/CompareCmd.cpp
cs/benchmark/cigma/trunk/src/FE.cpp
cs/benchmark/cigma/trunk/src/FE_Field.cpp
cs/benchmark/cigma/trunk/src/FE_Field.h
Log:
Cleaning up CompareCmd::run()
* Shape fn tabulations now happen earlier in FE::set_quadrature()
* Added FE_Field::tabulate_element()
* Calling field_a->tabulate_element() from main loop
* Added a few more reminders
Modified: cs/benchmark/cigma/trunk/src/CompareCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 05:47:32 UTC (rev 11190)
+++ cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 09:27:51 UTC (rev 11191)
@@ -204,16 +204,14 @@
assert(mesh != 0);
+ /* Now load the quadrature rule. If no rule is specified on the
+ * command line, a default rule is assigned based on the type of
+ * the cell. Also, an exception is thrown if the specified a rule
+ * does not conform geometrically to the interior of the cell.
+ */
quadratureIO.load(mesh->cell);
quadrature = quadratureIO.quadrature;
- if (quadrature == 0)
- {
- quadrature = field_a->fe->quadrature;
- }
assert(quadrature != 0);
-
- // XXX: perhaps quadrature data should be loaded first!
- // move this line back into FieldIO load() method
field_a->fe->set_quadrature(quadrature);
@@ -274,10 +272,10 @@
int valdim = field_a->n_rank();
// local data;
- //double *jxw = new double[nq];
+ //double *jxw = new double[nq]; // XXX: instead, use memory allocated by field_a->fe
double *jxw = field_a->fe->jxw;
double *dofs_a = new double[ndofs * valdim];
- double *phi_a = new double[nq * valdim];
+ double *phi_a = new double[nq * valdim]; // XXX: not needed when using tabulation
double *phi_b = new double[nq * valdim];
// norm
@@ -299,19 +297,6 @@
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
@@ -324,30 +309,19 @@
// ... calculate phi_a[]
- // XXX: using eval()
+ //
+ // XXX: time to move this main loop into its own function
+ // so we can use polymorphic dispatch on the argument types
+ //
+ // XXX: using eval() -- applicable in general (Field obj)
//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++)
- {
- 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;
- }
- }
+ // XXX: using tabulation -- applicable to FE_Field obj
+ field_a->tabulate_element(e, phi_a);
// ... calculate phi_b[]
for (q = 0; q < nq; q++)
@@ -356,7 +330,9 @@
}
// evaluate jacobian at known quadrature points
- field_a->fe->update_jxw();
+ field_a->fe->update_jxw(); // XXX: somehow, this method needs to be attached
+ // to the MeshPart object, so maybe it needs its
+ // own FE object after all?
// ... apply quadrature rule
double err = 0.0;
Modified: cs/benchmark/cigma/trunk/src/FE.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/FE.cpp 2008-02-20 05:47:32 UTC (rev 11190)
+++ cs/benchmark/cigma/trunk/src/FE.cpp 2008-02-20 09:27:51 UTC (rev 11191)
@@ -42,8 +42,15 @@
int dim = cell->n_celldim();
jxw = new double[nq];
+
basis_tab = new double[nq * ndofs];
basis_jet = new double[nq * ndofs * dim];
+
+ // get shape function values at known quadrature points
+ cell->shape(nq, quadrature->qpts, basis_tab);
+
+ // get shape function values at known quadrature points
+ cell->grad_shape(nq, quadrature->qpts, basis_jet);
}
// ---------------------------------------------------------------------------
Modified: cs/benchmark/cigma/trunk/src/FE_Field.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 05:47:32 UTC (rev 11190)
+++ cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 09:27:51 UTC (rev 11191)
@@ -71,23 +71,9 @@
cell->interpolate(field_dofs, uvw, value, valdim);
}
-void cigma::FE_Field::tabulate()
+void cigma::FE_Field::tabulate_element(int e, double *values)
{
- // 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;
Modified: cs/benchmark/cigma/trunk/src/FE_Field.h
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 05:47:32 UTC (rev 11190)
+++ cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 09:27:51 UTC (rev 11191)
@@ -31,8 +31,7 @@
void eval(double *point, double *value);
public:
- void tabulate();
- void eval(int e, double *values);
+ void tabulate_element(int e, double *values);
public:
void get_cell_dofs(int cellIndex, double *cellDofs);
More information about the cig-commits
mailing list