[cig-commits] r5068 - short/3D/PyLith/trunk/playpen/integrate/src
knepley at geodynamics.org
knepley at geodynamics.org
Wed Oct 18 07:55:28 PDT 2006
Author: knepley
Date: 2006-10-18 07:55:28 -0700 (Wed, 18 Oct 2006)
New Revision: 5068
Modified:
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
Log:
Added elastic operator formation to playpen
Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-10-17 18:44:53 UTC (rev 5067)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc 2006-10-18 14:55:28 UTC (rev 5068)
@@ -58,7 +58,7 @@
}
};
-void pylith::feassemble::Integrator::integrateFunction(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const value_type []))
+void pylith::feassemble::Integrator::integrateFunction_2d(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const value_type []))
{
const topology_type::patch_type patch = 0;
const Obj<topology_type>& topology = field->getTopology();
@@ -92,9 +92,48 @@
}
};
-void pylith::feassemble::Integrator::integrateLaplacianAction(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+void pylith::feassemble::Integrator::integrateFunction_3d(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const 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*3+0] + 1.0;
+ value_type eta = _points[q*3+1] + 1.0;
+ value_type zeta = _points[q*3+2] + 1.0;
+ _realCoords[0] = _Jac[0]*xi + _Jac[1]*eta + _Jac[2]*zeta + _v0[0];
+ _realCoords[1] = _Jac[3]*xi + _Jac[4]*eta + _Jac[5]*zeta + _v0[1];
+ _realCoords[2] = _Jac[6]*xi + _Jac[7]*eta + _Jac[8]*zeta + _v0[2];
+ 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);
+ }
+};
+
+void pylith::feassemble::Integrator::integrateFunction(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const value_type []))
+{
+ if (_dim == 2) {
+ this->integrateFunction_2d(field, coordinates, f);
+ } else if (_dim == 2) {
+ this->integrateFunction_3d(field, coordinates, f);
+ } else {
+ throw ALE::Exception("No integration function for this dimension");
+ }
+};
+
+void pylith::feassemble::Integrator::integrateLaplacianAction_2d(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ const topology_type::patch_type patch = 0;
const Obj<topology_type>& topology = X->getTopology();
const Obj<topology_type::label_sequence>& elements = topology->heightStratum(patch, 0);
const topology_type::label_sequence::iterator end = elements->end();
@@ -126,3 +165,113 @@
F->updateAdd(patch, *e_iter, _elemVector);
}
};
+
+void pylith::feassemble::Integrator::integrateLaplacianAction_3d(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ const topology_type::patch_type patch = 0;
+ const Obj<topology_type>& topology = X->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, _invJac, detJ);
+ // Element integral
+ PetscMemzero(_elemVector, NUM_BASIS_FUNCTIONS*sizeof(value_type));
+ PetscMemzero(_elemMatrix, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*sizeof(value_type));
+ for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ _testWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[6]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ _testWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[4]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[7]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ _testWork[2] = _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[5]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[8]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ _basisWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[6]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ _basisWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[4]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[7]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ _basisWork[2] = _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[5]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[8]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ _elemMatrix[i*NUM_BASIS_FUNCTIONS+j] += (_testWork[0]*_basisWork[0] + _testWork[1]*_basisWork[1] + _testWork[2]*_basisWork[2])*_weights[q]*detJ;
+ }
+ }
+ }
+ // Assembly
+ const value_type *xWork = X->restrict(patch, *e_iter);
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ _elemVector[i] += _elemMatrix[i*NUM_BASIS_FUNCTIONS+j]*xWork[j];
+ }
+ }
+ F->updateAdd(patch, *e_iter, _elemVector);
+ }
+};
+
+void pylith::feassemble::Integrator::integrateLaplacianAction(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ if (_dim == 2) {
+ this->integrateLaplacianAction_2d(X, F, coordinates);
+ } else if (_dim == 2) {
+ this->integrateLaplacianAction_3d(X, F, coordinates);
+ } else {
+ throw ALE::Exception("No integration function for this dimension");
+ }
+};
+
+void pylith::feassemble::Integrator::integrateElasticAction_3d(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ const topology_type::patch_type patch = 0;
+ const Obj<topology_type>& topology = X->getTopology();
+ const Obj<topology_type::label_sequence>& elements = topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator end = elements->end();
+ const double C = 1.0;
+ value_type detJ;
+
+ for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != end; ++e_iter) {
+ computeElementGeometry(coordinates, *e_iter, _v0, _Jac, _invJac, detJ);
+ // Element integral
+ PetscMemzero(_elemVector, NUM_BASIS_FUNCTIONS*3*sizeof(value_type));
+ PetscMemzero(_elemMatrix, NUM_BASIS_FUNCTIONS*NUM_BASIS_FUNCTIONS*9*sizeof(value_type));
+ for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ _testWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[6]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ _testWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[4]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[7]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ _testWork[2] = _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+0] + _invJac[5]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+1] + _invJac[8]*_basisDer[(q*NUM_BASIS_FUNCTIONS+i)*3+2];
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ _basisWork[0] = _invJac[0]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[3]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[6]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ _basisWork[1] = _invJac[1]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[4]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[7]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ _basisWork[2] = _invJac[2]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+0] + _invJac[5]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+1] + _invJac[8]*_basisDer[(q*NUM_BASIS_FUNCTIONS+j)*3+2];
+ for(int c = 0; c < 3; c++) {
+ for(int d = 0; d < 3; d++) {
+ for(int e = 0; e < 3; e++) {
+ for(int f = 0; f < 3; f++) {
+ _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += 0.25*C*(_testWork[d] + _testWork[c])*(_basisWork[f] +_basisWork[e])*_weights[q]*detJ;
+ }
+ }
+ _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS + j)*3+d] += _basis[q*NUM_BASIS_FUNCTIONS+i]*_basis[q*NUM_BASIS_FUNCTIONS+j]*_weights[q]*detJ;
+ }
+ }
+ }
+ }
+ }
+ // Assembly
+ const value_type *xWork = X->restrict(patch, *e_iter);
+ for(int i = 0; i < NUM_BASIS_FUNCTIONS; i++) {
+ for(int c = 0; c < 3; c++) {
+ for(int j = 0; j < NUM_BASIS_FUNCTIONS; j++) {
+ for(int d = 0; d < 3; d++) {
+ _elemVector[i*3+c] += _elemMatrix[((i*3+c)*NUM_BASIS_FUNCTIONS+j)*3+d]*xWork[j*3+d];
+ }
+ }
+ }
+ }
+ F->updateAdd(patch, *e_iter, _elemVector);
+ }
+};
+
+void pylith::feassemble::Integrator::integrateElasticAction(const Obj<section_type>& X, const Obj<section_type>& F, const Obj<section_type>& coordinates)
+{
+ if (_dim == 2) {
+ //this->integrateElasticAction_2d(X, F, coordinates);
+ } else if (_dim == 2) {
+ this->integrateElasticAction_3d(X, F, coordinates);
+ } else {
+ throw ALE::Exception("No integration function for this dimension");
+ }
+};
Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-10-17 18:44:53 UTC (rev 5067)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh 2006-10-18 14:55:28 UTC (rev 5068)
@@ -9,7 +9,7 @@
typedef ALE::Mesh Mesh;
typedef Mesh::topology_type topology_type;
typedef topology_type::point_type point_type;
- typedef Mesh::section_type section_type;
+ typedef Mesh::real_section_type section_type;
typedef section_type::value_type value_type;
protected:
const int NUM_QUADRATURE_POINTS;
@@ -37,8 +37,8 @@
_Jac = new value_type[dim*dim];
_invJac = new value_type[dim*dim];
_realCoords = new value_type[dim];
- _elemVector = new value_type[numBasisFuncs];
- _elemMatrix = new value_type[numBasisFuncs*numBasisFuncs];
+ _elemVector = new value_type[numBasisFuncs*dim];
+ _elemMatrix = new value_type[numBasisFuncs*numBasisFuncs*dim*dim];
_testWork = new value_type[numBasisFuncs];
_basisWork = new value_type[numBasisFuncs];
@@ -68,11 +68,18 @@
delete [] _testWork;
delete [] _basisWork;
};
+ protected:
+ void integrateFunction_2d(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
+ void integrateFunction_3d(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
+ void integrateLaplacianAction_2d(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
+ void integrateLaplacianAction_3d(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
+ void integrateElasticAction_3d(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
public:
void computeElementGeometry(const Obj<section_type>&, const point_type&, value_type [], value_type [], value_type [], value_type&);
void fillSection(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
void integrateFunction(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
void integrateLaplacianAction(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
+ void integrateElasticAction(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
};
}
}
Modified: short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc 2006-10-17 18:44:53 UTC (rev 5067)
+++ short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc 2006-10-18 14:55:28 UTC (rev 5068)
@@ -17,17 +17,17 @@
#include "Integration.hh"
// ----------------------------------------------------------------------
-ALE::Mesh::section_type::value_type zeroF(const ALE::Mesh::section_type::value_type coords[])
+ALE::Mesh::real_section_type::value_type zeroF(const ALE::Mesh::real_section_type::value_type coords[])
{
return 0.0;
}
-ALE::Mesh::section_type::value_type constantF(const ALE::Mesh::section_type::value_type coords[])
+ALE::Mesh::real_section_type::value_type constantF(const ALE::Mesh::real_section_type::value_type coords[])
{
return 1.0;
}
-ALE::Mesh::section_type::value_type linearF(const ALE::Mesh::section_type::value_type coords[])
+ALE::Mesh::real_section_type::value_type linearF(const ALE::Mesh::real_section_type::value_type coords[])
{
return coords[0];
}
@@ -53,12 +53,12 @@
iohandler.read(mesh, false);
const ALE::Mesh::topology_type::patch_type patch = 0;
- const Obj<ALE::Mesh::section_type>& coords = mesh->getSection("coordinates");
- const Obj<ALE::Mesh::section_type>& X = mesh->getSection("X");
- const Obj<ALE::Mesh::section_type>& F = mesh->getSection("F");
- pylith::feassemble::Integrator integrator(mesh->getDimension(),
- NUM_QUADRATURE_POINTS, points, weights,
- NUM_BASIS_FUNCTIONS, Basis, BasisDerivatives);
+ const Obj<ALE::Mesh::real_section_type>& coords = mesh->getRealSection("coordinates");
+ const Obj<ALE::Mesh::real_section_type>& X = mesh->getRealSection("X");
+ const Obj<ALE::Mesh::real_section_type>& F = mesh->getRealSection("F");
+ pylith::feassemble::Integrator integrator(mesh->getDimension(),
+ NUM_QUADRATURE_POINTS, points, weights,
+ NUM_BASIS_FUNCTIONS, Basis, BasisDerivatives);
coords->view("Coordinates");
X->setFiberDimensionByDepth(patch, 0, 1);
More information about the cig-commits
mailing list