[cig-commits] r11208 - cs/benchmark/cigma/trunk/src
luis at geodynamics.org
luis at geodynamics.org
Wed Feb 20 13:12:21 PST 2008
Author: luis
Date: 2008-02-20 13:12:20 -0800 (Wed, 20 Feb 2008)
New Revision: 11208
Modified:
cs/benchmark/cigma/trunk/src/CompareCmd.cpp
cs/benchmark/cigma/trunk/src/FE.cpp
cs/benchmark/cigma/trunk/src/FE.h
cs/benchmark/cigma/trunk/src/FE_Field.cpp
cs/benchmark/cigma/trunk/src/FE_Field.h
cs/benchmark/cigma/trunk/src/QuadratureRule.cpp
cs/benchmark/cigma/trunk/src/QuadratureRule.h
Log:
Changed cigma::FE to inherit from cigma::QuadratureRule
* Adds tabulation functionality to QuadratureRule
* Intended for use with FE_Field::tabulate_element()
Modified: cs/benchmark/cigma/trunk/src/CompareCmd.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/CompareCmd.cpp 2008-02-20 21:12:20 UTC (rev 11208)
@@ -47,7 +47,12 @@
cigma::CompareCmd::~CompareCmd()
{
- if (qr != 0) delete qr;
+ /* XXX: for now, don't delete
+ if (qr != 0)
+ {
+ delete qr; // XXX
+ } // */
+
if (residuals != 0) delete residuals;
}
@@ -214,25 +219,25 @@
firstIO.load();
field_a = firstIO.field;
assert(field_a != 0);
- field_a->fe = new FE();
- field_a->fe->set_cell(field_a->meshPart->cell);
+ //field_a->fe = new FE();
+ //field_a->fe->set_cell(field_a->meshPart->cell);
cout << "first field path = " << firstIO.field_path << endl;
cout << "first field dimensions = "
<< field_a->meshPart->nel << " cells, "
<< field_a->meshPart->nno << " nodes, "
- << field_a->fe->cell->n_nodes() << " dofs/cell, "
+ << field_a->meshPart->cell->n_nodes() << " dofs/cell, "
<< "rank " << field_a->n_rank() << endl;
secondIO.load();
field_b = secondIO.field;
assert(field_b != 0);
- field_b->fe = new FE();
- field_b->fe->set_cell(field_b->meshPart->cell);
+ //field_b->fe = new FE();
+ //field_b->fe->set_cell(field_b->meshPart->cell);
cout << "second field path = " << secondIO.field_path << endl;
cout << "second field dimensions = "
<< field_b->meshPart->nel << " cells, "
<< field_b->meshPart->nno << " nodes, "
- << field_b->fe->cell->n_nodes() << " dofs/cell, "
+ << field_b->meshPart->cell->n_nodes() << " dofs/cell, "
<< "rank " << field_b->n_rank() << endl;
@@ -261,12 +266,37 @@
quadratureIO.load(mesh->cell);
quadrature = quadratureIO.quadrature;
assert(quadrature != 0);
- field_a->fe->set_quadrature(quadrature); // XXX: eliminate need for this statement
+ //field_a->fe->set_quadrature(quadrature); // XXX: eliminate need for this statement
+
+
+ field_a->fe = new FE();
+ field_a->fe->set_mesh(field_a->meshPart);
+ //field_a->fe->set_quadrature_points(field_a->meshPart->cell->default_quadrature_points);
+ field_a->fe->set_quadrature_points(quadrature); // XXX
+
+
+ field_b->fe = new FE();
+ field_b->fe->set_mesh(field_b->meshPart);
+ //field_b->fe->set_quadrature_points(field_b->meshPart->cell->default_quadrature_points);
+
+
+ // XXX: if field_a has a quadrature rule, reuse it
+ // XXX: however, this complicates our destructor...
+ qr = field_a->fe;
+ assert(field_a->meshPart == mesh); // XXX
+
+
+ /* XXX: when neither field_a or field_b have a mesh, need own QuadratureRule instance
+ assert(mesh != 0);
+ assert(field_a->meshPart == 0);
+ assert(field_b->meshPart == 0);
qr = new QuadratureRule();
qr->set_mesh(mesh);
qr->set_quadrature_points(quadrature);
+ // */
+
return;
}
@@ -293,6 +323,7 @@
//qr->set_quadrature_points(quadrature);
QuadratureRule *qr = env->qr; // XXX: replace MeshPart fn arg by QuadratureRule
assert(qr != 0);
+ assert(qr->points != 0);
// dimensions
int nel = qr->meshPart->nel;
@@ -407,8 +438,9 @@
// start with the obvious checks
+ assert(quadrature != 0);
assert(mesh != 0);
- assert(quadrature != 0);
+ assert(qr != 0);
assert(field_a != 0);
assert(field_b != 0);
Modified: cs/benchmark/cigma/trunk/src/FE.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/FE.cpp 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/FE.cpp 2008-02-20 21:12:20 UTC (rev 11208)
@@ -5,65 +5,36 @@
cigma::FE::FE()
{
- quadrature = 0;
- cell = 0;
-
- jxw = 0;
basis_tab = 0;
basis_jet = 0;
}
cigma::FE::~FE()
{
- if (jxw != 0) delete [] jxw;
if (basis_tab != 0) delete [] basis_tab;
if (basis_jet != 0) delete [] basis_jet;
}
// ---------------------------------------------------------------------------
-void cigma::FE::set_cell(Cell *cell)
+void cigma::FE::set_quadrature_points(QuadraturePoints *pts)
{
- assert(cell != 0);
- this->cell = cell;
-}
+ assert(meshPart != 0);
+ assert(meshPart->cell != 0);
-void cigma::FE::set_quadrature(QuadraturePoints *quadrature)
-{
- assert(cell != 0);
+ this->QuadratureRule::set_quadrature_points(pts);
- this->quadrature = quadrature;
+ assert(points != 0);
+ const int nq = points->n_points();
+ const int ndofs = meshPart->cell->n_nodes();
- assert(quadrature->n_points() > 0);
- assert(quadrature->n_dim() == cell->n_celldim());
-
- int nq = quadrature->n_points();
- int ndofs = cell->n_nodes();
- int dim = cell->n_celldim();
-
- jxw = new double[nq];
-
// get shape function values at known quadrature points
basis_tab = new double[nq * ndofs];
- cell->shape(nq, quadrature->qpts, basis_tab);
+ meshPart->cell->shape(nq, points->qpts, basis_tab);
// get shape function values at known quadrature points
//basis_jet = new double[nq * ndofs * dim];
- //cell->grad_shape(nq, quadrature->qpts, basis_jet);
+ //meshPart->cell->grad_shape(nq, quadrature->qpts, basis_jet);
}
// ---------------------------------------------------------------------------
-
-void cigma::FE::update_jxw()
-{
- const int nq = quadrature->n_points();
- const int celldim = cell->n_celldim();
- for (int q = 0; q < nq; q++)
- {
- double jac[3][3];
- jxw[q] = quadrature->weight(q) * cell->jacobian((*quadrature)[q], jac);
- }
-}
-
-// ---------------------------------------------------------------------------
-
Modified: cs/benchmark/cigma/trunk/src/FE.h
===================================================================
--- cs/benchmark/cigma/trunk/src/FE.h 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/FE.h 2008-02-20 21:12:20 UTC (rev 11208)
@@ -4,33 +4,27 @@
#include "Cell.h"
#include "Points.h"
#include "QuadraturePoints.h"
+#include "QuadratureRule.h"
namespace cigma
{
- class FE;
+ class FE; // XXX: rename to FE_QR since now we derive from QuadratureRule
}
/**
* @brief Tabulation for Finite Elements;
*
*/
-class cigma::FE
+class cigma::FE : public QuadratureRule
{
public:
FE();
~FE();
public:
- void set_cell(Cell *cell);
- void set_quadrature(QuadraturePoints *quadrature);
+ void set_quadrature_points(QuadraturePoints *pts);
public:
- void update_jxw();
-
-public:
- Cell *cell;
- QuadraturePoints *quadrature;
- double *jxw;
double *basis_tab; // [nq x ndofs]
double *basis_jet; // [nq x ndofs x celldim]
};
Modified: cs/benchmark/cigma/trunk/src/FE_Field.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/FE_Field.cpp 2008-02-20 21:12:20 UTC (rev 11208)
@@ -20,13 +20,13 @@
void cigma::FE_Field::get_cell_dofs(int cellIndex, double *cellDofs)
{
- assert(fe != 0);
assert(dofHandler != 0);
assert(meshPart != 0);
+
assert(0 <= cellIndex);
assert(cellIndex < (meshPart->nel));
- const int ndofs = fe->cell->n_nodes();
+ const int ndofs = meshPart->cell->n_nodes();
const int valdim = n_rank();
int *n = &(meshPart->connect[ndofs * cellIndex]);
@@ -47,7 +47,7 @@
assert(meshPart != 0);
// get reference cell object from fe
- Cell *cell = fe->cell;
+ Cell *cell = meshPart->cell;
assert(cell != 0);
// find the cell which contains given point
@@ -72,17 +72,17 @@
{
assert(fe != 0);
assert(meshPart != 0);
+ assert(fe->meshPart == meshPart);
+ assert(fe->points != 0);
- Cell *cell = fe->cell;
+ Cell *cell = meshPart->cell;
assert(cell != 0);
- QuadraturePoints *quadrature = fe->quadrature;
- assert(quadrature != 0);
- int nq = quadrature->n_points();
- double *qpts = quadrature->qpts;
+ QuadraturePoints *points = fe->points;
+ assert(points != 0);
+ int nq = points->n_points();
// tabulate the function values
- int i,j;
const int valdim = n_rank();
const int ndofs = cell->n_nodes();
double dofs[ndofs * valdim]; // XXX
@@ -93,10 +93,10 @@
for (int q = 0; q < nq; q++)
{
double *N = &(fe->basis_tab[ndofs*q]);
- for (i = 0; i < valdim; i++)
+ for (int i = 0; i < valdim; i++)
{
double sum = 0.0;
- for (j = 0; j < ndofs; j++)
+ for (int j = 0; j < ndofs; j++)
{
sum += dofs[i + valdim*j] * N[j];
}
Modified: cs/benchmark/cigma/trunk/src/FE_Field.h
===================================================================
--- cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/FE_Field.h 2008-02-20 21:12:20 UTC (rev 11208)
@@ -2,9 +2,9 @@
#define __FE_FIELD_H__
#include "Field.h"
+#include "DofHandler.h"
+#include "MeshPart.h"
#include "FE.h"
-#include "MeshPart.h"
-#include "DofHandler.h"
namespace cigma
{
@@ -31,16 +31,18 @@
void eval(double *point, double *value);
public:
+ void set_quadrature_rule(QuadratureRule *rule);
void tabulate_element(int e, double *values);
public:
void get_cell_dofs(int cellIndex, double *cellDofs);
public:
+
int dim;
int rank;
+ DofHandler *dofHandler;
MeshPart *meshPart;
- DofHandler *dofHandler;
FE *fe;
};
Modified: cs/benchmark/cigma/trunk/src/QuadratureRule.cpp
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadratureRule.cpp 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/QuadratureRule.cpp 2008-02-20 21:12:20 UTC (rev 11208)
@@ -10,10 +10,12 @@
{
meshPart = 0;
points = 0;
+ jxw = 0;
}
cigma::QuadratureRule::~QuadratureRule()
{
+ if (jxw != 0) delete [] jxw;
}
@@ -32,9 +34,8 @@
assert(meshPart != 0);
assert(meshPart->cell != 0);
+ this->points = points;
assert(points != 0);
-
- this->points = points;
assert(points->n_points() > 0);
assert(points->n_dim() == meshPart->cell->n_celldim());
Modified: cs/benchmark/cigma/trunk/src/QuadratureRule.h
===================================================================
--- cs/benchmark/cigma/trunk/src/QuadratureRule.h 2008-02-20 21:12:17 UTC (rev 11207)
+++ cs/benchmark/cigma/trunk/src/QuadratureRule.h 2008-02-20 21:12:20 UTC (rev 11208)
@@ -20,7 +20,7 @@
public:
void set_mesh(MeshPart *mesh);
- void set_quadrature_points(QuadraturePoints *pts);
+ virtual void set_quadrature_points(QuadraturePoints *pts);
public:
void select_cell(int e);
More information about the cig-commits
mailing list