[cig-commits] r4545 - short/3D/PyLith/trunk/playpen/integrate/src
knepley at geodynamics.org
knepley at geodynamics.org
Thu Sep 14 19:28:25 PDT 2006
Author: knepley
Date: 2006-09-14 19:28:24 -0700 (Thu, 14 Sep 2006)
New Revision: 4545
Added:
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
Log:
Added source for integrator
Added: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-09-15 02:27:46 UTC (rev 4544)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-09-15 02:28:24 UTC (rev 4545)
@@ -0,0 +1,72 @@
+#include "elemVector.hh"
+
+void pylith::feassemble::Integrator::computeElementGeometry(const Obj<section_type>& coordinates, const point_type& e, value_type v0[], value_type J[], value_type invJ[], value_type& detJ)
+{
+ const topology_type::patch_type patch = 0;
+ const value_type *coords = coordinates->restrict(patch, e);
+ value_type invDet;
+
+ for(int d = 0; d < _dim; d++) {
+ v0[d] = coords[d];
+ }
+ for(int d = 0; d < _dim; d++) {
+ for(int f = 0; f < _dim; f++) {
+ J[d*_dim+f] = 0.5*(coords[(f+1)*_dim+d] - coords[0*_dim+d]);
+ }
+ }
+ if (_dim == 1) {
+ detJ = J[0];
+ } else if (_dim == 2) {
+ detJ = J[0]*J[3] - J[1]*J[2];
+ } else if (_dim == 3) {
+ detJ = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
+ J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
+ J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
+ }
+ if (invJ) {
+ invDet = 1.0/detJ;
+ if (_dim == 2) {
+ invJ[0] = invDet*J[3];
+ invJ[1] = -invDet*J[1];
+ invJ[2] = -invDet*J[2];
+ invJ[3] = invDet*J[0];
+ } else if (_dim == 3) {
+ // FIX: This may be wrong
+ invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
+ invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
+ invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
+ invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
+ invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
+ invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
+ invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
+ invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
+ invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
+ }
+ }
+};
+
+void pylith::feassemble::Integrator::integrateFunction(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(value_type []))
+{
+ const topology_type::patch_type patch = 0;
+ const Obj<topology_type>& topology = field->getTopology();
+ const Obj<topology_type::label_sequence>& elements = topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator end = elements->end();
+ value_type detJ;
+
+ for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != end; ++e_iter) {
+ computeElementGeometry(coordinates, *e_iter, _v0, _Jac, PETSC_NULL, detJ);
+ // Element integral
+ PetscMemzero(_elemVector, NUM_BASIS_FUNCTIONS*sizeof(value_type));
+ for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
+ value_type xi = _points[q*2+0] + 1.0;
+ value_type eta = _points[q*2+1] + 1.0;
+ _realCoords[0] = _Jac[0]*xi + _Jac[1]*eta + _v0[0];
+ _realCoords[1] = _Jac[2]*xi + _Jac[3]*eta + _v0[1];
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ _elemVector[i] += _basis[q*NUM_BASIS_FUNCTIONS+i]*f(_realCoords)*_weights[q]*detJ;
+ }
+ }
+ /* Assembly */
+ field->updateAdd(patch, *e_iter, _elemVector);
+ }
+};
Added: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-09-15 02:27:46 UTC (rev 4544)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-09-15 02:28:24 UTC (rev 4545)
@@ -0,0 +1,62 @@
+#include <Mesh.hh>
+
+using ALE::Obj;
+
+namespace pylith {
+ namespace feassemble {
+ class Integrator {
+ public:
+ typedef ALE::Mesh Mesh;
+ typedef Mesh::topology_type topology_type;
+ typedef topology_type::point_type point_type;
+ typedef Mesh::section_type section_type;
+ typedef section_type::value_type value_type;
+ protected:
+ const int NUM_QUADRATURE_POINTS;
+ const int NUM_BASIS_FUNCTIONS;
+ const int _dim;
+ value_type *_points;
+ value_type *_weights;
+ value_type *_basis;
+ value_type *_v0;
+ value_type *_Jac;
+ value_type *_invJac;
+ value_type *_realCoords;
+ value_type *_elemVector;
+ public:
+ Integrator(int dim, int numQuadPoints, value_type points[], value_type weights[], int numBasisFuncs, value_type basis[]) : NUM_QUADRATURE_POINTS(numQuadPoints), NUM_BASIS_FUNCTIONS(numBasisFuncs), _dim(dim) {
+ _points = new value_type[numQuadPoints*dim];
+ _weights = new value_type[numQuadPoints];
+ _basis = new value_type[numQuadPoints*numBasisFuncs];
+ _v0 = new value_type[dim];
+ _Jac = new value_type[dim*dim];
+ _invJac = new value_type[dim*dim];
+ _realCoords = new value_type[dim];
+ _elemVector = new value_type[numBasisFuncs];
+
+ for(int q = 0; q < numQuadPoints; ++q) {
+ for(int d = 0; d < dim; ++d) {
+ _points[q*dim+d] = points[q*dim+d];
+ }
+ _weights[q] = weights[q];
+ for(int b = 0; b < numBasisFuncs; ++b) {
+ _basis[q*numBasisFuncs+b] = basis[q*numBasisFuncs+b];
+ }
+ }
+ };
+ ~Integrator() {
+ delete [] _points;
+ delete [] _weights;
+ delete [] _basis;
+ delete [] _v0;
+ delete [] _Jac;
+ delete [] _invJac;
+ delete [] _realCoords;
+ delete [] _elemVector;
+ };
+ public:
+ void computeElementGeometry(const Obj<section_type>&, const point_type&, value_type [], value_type [], value_type [], value_type&);
+ void integrateFunction(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(value_type []));
+ };
+ }
+}
More information about the cig-commits
mailing list