[cig-commits] r6455 - in short/3D/PyLith/trunk: libsrc/feassemble
modulesrc/feassemble pylith/feassemble pylith/problems
brad at geodynamics.org
brad at geodynamics.org
Wed Mar 28 17:50:21 PDT 2007
Author: brad
Date: 2007-03-28 17:50:21 -0700 (Wed, 28 Mar 2007)
New Revision: 6455
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py
short/3D/PyLith/trunk/pylith/problems/Explicit.py
Log:
Worked on implementing integration actions in ExplicitElasticity.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.cc 2007-03-29 00:50:21 UTC (rev 6455)
@@ -57,15 +57,15 @@
} // material
// ----------------------------------------------------------------------
-// Integrate residual term (b) for dynamic elasticity term for 3-D
+// Integrate constant term (b) for dynamic elasticity term for 3-D
// finite elements.
void
-pylith::feassemble::ExplicitElasticity::integrateResidual(
- const ALE::Obj<real_section_type>& residual,
+pylith::feassemble::ExplicitElasticity::integrateConstant(
+ const ALE::Obj<real_section_type>& b,
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<real_section_type>& dispTmdt,
const ALE::Obj<real_section_type>& coordinates)
-{ // integrateResidual
+{ // integrateConstant
assert(0 != _quadrature);
// Get information about section
@@ -77,10 +77,18 @@
// Get parameters used in integration.
const double dt = _dt;
+ const double dt2 = dt*dt;
// Allocate vector for cell values (if necessary)
_initCellVector();
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double* quadWts = _quadrature->quadWts();
+ const int numBasis = _quadrature->numCorners();
+ const int spaceDim = _quadrature->spaceDim();
+ const int cellDim = _quadrature->cellDim();
+
for (topology_type::label_sequence::iterator cellIter=cells->begin();
cellIter != cellsEnd;
++cellIter) {
@@ -96,24 +104,20 @@
const real_section_type::value_type* dispTmdtCell =
dispTmdt->restrict(patch, *cellIter);
- // Get cell geometry information
- const int numQuadPts = _quadrature->numQuadPts();
+ // Get cell geometry information that depends on cell
const double* basis = _quadrature->basis();
- const double* quadPts = _quadrature->quadPts();
- const double* quadWts = _quadrature->quadWts();
+ const double* basisDeriv = _quadrature->basisDeriv();
const double* jacobianDet = _quadrature->jacobianDet();
- const int numBasis = _quadrature->numCorners();
- const int spaceDim = _quadrature->spaceDim();
// Get material physical properties at quadrature points for this cell
_material->calcProperties(*cellIter, patch, numQuadPts);
const double* density = _material->density();
const double* elasticConsts = _material->elasticConsts();
+ const int numElasticConsts = _material->numElasticConsts();
// Compute action for cell
// Compute action for inertial terms
- const double dt2 = dt*dt;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt =
quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
@@ -125,24 +129,171 @@
const double valIJ = valI * basis[iQ+jBasis];
for (int iDim=0; iDim < spaceDim; ++iDim)
_cellVector[iBlock+iDim] +=
- valIJ * 2.0 * (dispTCell[jBlock+iDim] -
- dispTmdtCell[jBlock+iDim]);
+ valIJ * (2.0 * dispTCell[jBlock+iDim] -
+ dispTmdtCell[jBlock+iDim]);
} // for
} // for
} // for
-
- // Compute action for elastic terms
- // ADD STUFF HERE
-
PetscErrorCode err =
PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+4*spaceDim))));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
+ /** :TODO:
+ *
+ * If cellDim and spaceDim are different, we need to transform
+ * displacements into cellDim, compute action, and transform
+ * result back into spaceDim. Can we get this from the inverse of
+ * the Jacobian?
+ */
+ if (cellDim != spaceDim)
+ throw std::logic_error("Not implemented yet.");
+
+ // Compute action for elastic terms
+ if (1 == cellDim) {
+ assert(1 == numElasticConsts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const double C1111 = elasticConsts[iQuad];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double valI = wt*basisDeriv[iQ+iBasis]*C1111;
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double valIJ = valI * basisDeriv[iQ+jBasis];
+ _cellVector[iBlock] -= valIJ * dispTCell[jBlock];
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*3)));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } else if (2 == cellDim) {
+ assert(6 == numElasticConsts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const int iConst = iQuad*numElasticConsts;
+ const double C1111 = elasticConsts[iConst+0];
+ const double C1122 = elasticConsts[iConst+1];
+ const double C1112 = elasticConsts[iConst+2];
+ const double C2222 = elasticConsts[iConst+3];
+ const double C2212 = elasticConsts[iConst+4];
+ const double C1212 = elasticConsts[iConst+5];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp;
+ _cellVector[iBlock ] -=
+ ki0j0 * dispTCell[jBlock ] + ki0j1 * dispTCell[jBlock+1];
+ _cellVector[iBlock+1] -=
+ ki0j1 * dispTCell[jBlock ] + ki1j1 * dispTCell[jBlock+1];
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(1+numBasis*(2+numBasis*(2+3*11+2*4))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } else if (3 == cellDim) {
+ assert(21 == numElasticConsts);
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ const int iConst = iQuad*numElasticConsts;
+ const double C1111 = elasticConsts[iConst+ 0];
+ const double C1122 = elasticConsts[iConst+ 1];
+ const double C1133 = elasticConsts[iConst+ 2];
+ const double C1112 = elasticConsts[iConst+ 3];
+ const double C1123 = elasticConsts[iConst+ 4];
+ const double C1113 = elasticConsts[iConst+ 5];
+ const double C2222 = elasticConsts[iConst+ 6];
+ const double C2233 = elasticConsts[iConst+ 7];
+ const double C2212 = elasticConsts[iConst+ 8];
+ const double C2223 = elasticConsts[iConst+ 9];
+ const double C2213 = elasticConsts[iConst+10];
+ const double C3333 = elasticConsts[iConst+11];
+ const double C3312 = elasticConsts[iConst+12];
+ const double C3323 = elasticConsts[iConst+13];
+ const double C3313 = elasticConsts[iConst+14];
+ const double C1212 = elasticConsts[iConst+15];
+ const double C1223 = elasticConsts[iConst+16];
+ const double C1213 = elasticConsts[iConst+17];
+ const double C2323 = elasticConsts[iConst+18];
+ const double C2313 = elasticConsts[iConst+19];
+ const double C1313 = elasticConsts[iConst+20];
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = iBasis * spaceDim;
+ const double Nip = wt*basisDeriv[iQ+iBasis*cellDim+0];
+ const double Niq = wt*basisDeriv[iQ+iBasis*cellDim+1];
+ const double Nir = wt*basisDeriv[iQ+iBasis*cellDim+2];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = jBasis * spaceDim;
+ const double Njp = basisDeriv[iQ+jBasis*cellDim+0];
+ const double Njq = basisDeriv[iQ+jBasis*cellDim+1];
+ const double Njr = basisDeriv[iQ+jBasis*cellDim+2];
+ const double ki0j0 =
+ C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
+ C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
+ const double ki0j1 =
+ C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
+ C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
+ const double ki0j2 =
+ C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
+ C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
+ C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
+ const double ki1j1 =
+ C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
+ C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
+ const double ki1j2 =
+ C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
+ C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
+ C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
+ const double ki2j2 =
+ C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
+ C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
+ C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
+
+ _cellVector[iBlock ] -=
+ ki0j0 * dispTCell[jBlock ] +
+ ki0j1 * dispTCell[jBlock+1] +
+ ki0j2 * dispTCell[jBlock+2];
+ _cellVector[iBlock+1] -=
+ ki0j1 * dispTCell[jBlock ] +
+ ki1j1 * dispTCell[jBlock+1] +
+ ki1j2 * dispTCell[jBlock+2];
+ _cellVector[iBlock+2] -=
+ ki0j2 * dispTCell[jBlock ] +
+ ki1j2 * dispTCell[jBlock+1] +
+ ki2j2 * dispTCell[jBlock+2];
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(1+numBasis*(3+numBasis*(3+6*26+3*6))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+ } // if/else
+
// Assemble cell contribution into field
- residual->updateAdd(patch, *cellIter, _cellVector);
+ b->updateAdd(patch, *cellIter, _cellVector);
} // for
-} // integrateResidual
+} // integrateConstant
// ----------------------------------------------------------------------
// Compute matrix associated with operator.
@@ -152,27 +303,70 @@
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<real_section_type>& coordinates)
{ // integrateJacobian
-} // integrateJacobian
+ assert(0 != _quadrature);
-// ----------------------------------------------------------------------
-// Compute lumped matrix associated with operator.
-void
-pylith::feassemble::ExplicitElasticity::integrateJacobian(
- const ALE::Obj<real_section_type>& fieldOut,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& coordinates)
-{ // integrateJacobian
+ // Get information about section
+ const topology_type::patch_type patch = 0;
+ const ALE::Obj<topology_type>& topology = dispT->getTopology();
+ const ALE::Obj<topology_type::label_sequence>& cells =
+ topology->heightStratum(patch, 0);
+ const topology_type::label_sequence::iterator cellsEnd = cells->end();
+
+ // Get parameters used in integration.
+ const double dt = _dt;
+ const double dt2 = dt*dt;
+
+ // Allocate vector for cell values (if necessary)
+ _initCellMatrix();
+
+ // Get cell geometry information that doesn't depend on cell
+ const int numQuadPts = _quadrature->numQuadPts();
+ const double* quadWts = _quadrature->quadWts();
+ const int numBasis = _quadrature->numCorners();
+ const int spaceDim = _quadrature->spaceDim();
+ for (topology_type::label_sequence::iterator cellIter=cells->begin();
+ cellIter != cellsEnd;
+ ++cellIter) {
+ // Compute geometry information for current cell
+ _quadrature->computeGeometry(coordinates, *cellIter);
+
+ // Reset element vector to zero
+ _resetCellMatrix();
+
+ // Get cell geometry information that depends on cell
+ const double* basis = _quadrature->basis();
+ const double* jacobianDet = _quadrature->jacobianDet();
+
+ // Get material physical properties at quadrature points for this cell
+ _material->calcProperties(*cellIter, patch, numQuadPts);
+ const double* density = _material->density();
+
+ // Compute Jacobian for cell
+
+ // Compute Jacobian for inertial terms
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt =
+ quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad] / dt2;
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ const int iBlock = (iBasis * spaceDim) * (spaceDim * numBasis);
+ const double valI = wt*basis[iQ+iBasis];
+ for (int jBasis=0; jBasis < numBasis; ++jBasis) {
+ const int jBlock = (jBasis * spaceDim);
+ const double valIJ = valI * basis[iQ+jBasis];
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ _cellMatrix[iBlock+jBlock+iDim*spaceDim+iDim] += valIJ;
+ } // for
+ } // for
+ } // for
+ PetscErrorCode err =
+ PetscLogFlops(numQuadPts*(3+numBasis*(1+numBasis*(1+spaceDim))));
+ if (err)
+ throw std::runtime_error("Logging PETSc flops failed.");
+
+ // Assemble cell contribution into field
+ //mat->updateAdd(patch, *cellIter, _cellVector);
+ } // for
} // integrateJacobian
-// ----------------------------------------------------------------------
-// Setup material property parameters by querying database.
-void
-pylith::feassemble::ExplicitElasticity::initialize(ALE::Obj<ALE::Mesh>& mesh,
- spatialdata::geocoords::CoordSys* cs)
-{ // initialize
- assert(0 != _material);
- _material->initialize(mesh, cs, _quadrature);
-} // initialize
-
// End of file
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ExplicitElasticity.hh 2007-03-29 00:50:21 UTC (rev 6455)
@@ -94,24 +94,24 @@
*/
void material(const materials::ElasticMaterial* m);
- /** Integrate residual term (b) for dynamic elasticity term
+ /** Integrate constant term (b) for dynamic elasticity term
* for 3-D finite elements.
*
* Compute b = 2/(dt*dt)[M]{u(t) - 1/(dt*dt)[M]{u(t-dt)} - [K]{u(t)}, where
* [M] = density * [N]^T [N]
*
*
- * @param residual Residual field (output)
+ * @param b Constant field (output)
* @param dispT Displacement field at time t
* @param dispTmdt Displacement field at time t-dt
* @param coordinates Field of cell vertex coordinates
*/
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ void integrateConstant(const ALE::Obj<real_section_type>& b,
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<real_section_type>& dispTmdt,
const ALE::Obj<real_section_type>& coordinates);
- /** Compute matrix (A) associated with operator.
+ /** Compute Jacobian matrix (A) associated with operator.
*
* @param mat Sparse matrix
* @param dispT Displacement at time t
@@ -121,24 +121,6 @@
const ALE::Obj<real_section_type>& dispT,
const ALE::Obj<real_section_type>& coordinates);
- /** Compute field (A) associated with lumped operator.
- *
- * @param fieldOut Output Jacobian field
- * @param dispT Displacement at time t
- * @param coordinates Field of cell vertex coordinates
- */
- void integrateJacobian(const ALE::Obj<real_section_type>& fieldOut,
- const ALE::Obj<real_section_type>& dispT,
- const ALE::Obj<real_section_type>& coordinates);
-
- /** Setup material property parameters by querying database.
- *
- * @param mesh PETSc mesh
- * @param cs Pointer to coordinate system of vertices
- */
- void initialize(ALE::Obj<ALE::Mesh>& mesh,
- spatialdata::geocoords::CoordSys* cs);
-
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorExplicit.hh 2007-03-29 00:50:21 UTC (rev 6455)
@@ -73,21 +73,21 @@
virtual
double stableTimeStep(void) const;
- /** Integrate residual term (b) for dynamic elasticity term
+ /** Integrate constant term (b) for dynamic elasticity term
* for finite elements.
*
- * @param residual Residual field (output)
+ * @param fieldOut Constant field (output)
* @param fieldInT Input field at time t
* @param fieldInTmdt Input field at time t-dt
* @param coordinates Field of cell vertex coordinates
*/
virtual
- void integrateResidual(const ALE::Obj<real_section_type>& residual,
+ void integrateConstant(const ALE::Obj<real_section_type>& fieldOut,
const ALE::Obj<real_section_type>& fieldInT,
const ALE::Obj<real_section_type>& fieldInTmdt,
const ALE::Obj<real_section_type>& coordinates) = 0;
- /** Compute matrix (A) associated with operator.
+ /** Compute Jacobian matrix (A) associated with operator.
*
* @param mat Sparse matrix
* @param fieldIn Input field at time t
@@ -97,18 +97,7 @@
void integrateJacobian(PetscMat* mat,
const ALE::Obj<real_section_type>& fieldIn,
const ALE::Obj<real_section_type>& coordinates) = 0;
-
- /** Compute field (A) associated with lumped operator.
- *
- * @param fieldOut Output Jacobian field
- * @param fieldInT Input field at time t
- * @param coordinates Field of cell vertex coordinates
- */
- virtual
- void integrateJacobian(const ALE::Obj<real_section_type>& fieldOut,
- const ALE::Obj<real_section_type>& fieldInT,
- const ALE::Obj<real_section_type>& coordinates) = 0;
-
+
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2007-03-29 00:50:21 UTC (rev 6455)
@@ -222,9 +222,9 @@
* size = numQuadPts * numCorners
* index = iQuadPt*numCorners + iBasis
*/
- double* _basis; ///< Array of basis fns evaluated at quad pts
+ double* _basis;
- /** Array of basis functions evaluated at the quadrature points.
+ /** Array of basis function derivatives evaluated at the quadrature points.
*
* N0xQp0, N0yQp0, N0zQp0, N1xQp0, N1yQp0, N1zQp0, ...
* N0xQp1, N0yQp1, N0zQp1, N1xQp1, N1yQp1, N1zQp1, ...
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2007-03-29 00:50:21 UTC (rev 6455)
@@ -508,42 +508,6 @@
return
- def initialize(self, mesh, cs):
- """
- Initialize integrator.
-
- @param mesh PETSc mesh
- @param cs Coordinate system associated with mesh vertices
- """
- # create shim for method 'initialize'
- #embed{ void ExplicitElasticity_initialize(void* pObj, void* meshObj, void* csObj)
- try {
- ALE::Obj<ALE::Mesh>* mesh = (ALE::Obj<ALE::Mesh>*) meshObj;
- spatialdata::geocoords::CoordSys* cs =
- (spatialdata::geocoords::CoordSys*) csObj;
- ((pylith::feassemble::ExplicitElasticity*) pObj)->initialize(*mesh, cs);
- } catch (const std::exception& err) {
- PyErr_SetString(PyExc_RuntimeError,
- const_cast<char*>(err.what()));
- } catch (...) {
- PyErr_SetString(PyExc_RuntimeError,
- "Caught unknown C++ exception.");
- } // try/catch
- #}embed
-
- if not mesh.name == "pylith_topology_Mesh":
- raise TypeError, \
- "Argument must be extension module type " \
- "'pylith::topology::Mesh'."
- if not cs.name == "spatialdata_geocoords_CoordSys":
- raise TypeError, \
- "Argument must be extension module type " \
- "'spatialdata::geocoords::CoordSys'."
- ExplicitElasticity_initialize(self.thisptr, ptrFromHandle(mesh),
- ptrFromHandle(cs))
- return
-
-
property material:
def __set__(self, m):
"""
Modified: short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/pylith/feassemble/ExplicitElasticity.py 2007-03-29 00:50:21 UTC (rev 6455)
@@ -46,10 +46,8 @@
self._info.log("Initializing integrator for material '%s'." % \
material.label)
material.initialize(mesh)
-
self.material = material
self.cppHandle.material = self.material.cppHandle
- #self.cppHandle.createParameters(mesh.cppHandle)
return
Modified: short/3D/PyLith/trunk/pylith/problems/Explicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-03-29 00:39:37 UTC (rev 6454)
+++ short/3D/PyLith/trunk/pylith/problems/Explicit.py 2007-03-29 00:50:21 UTC (rev 6455)
@@ -65,15 +65,22 @@
"""
Create explicit integrators for each element family.
"""
- self._info.log("Initializing integrators.")
from pylith.feassemble.ExplicitElasticity import ExplicitElasticity
+ self._info.log("Initializing integrators.")
self.integrators = []
for material in materialsBin.materials:
integrator = ExplicitElasticity()
integrator.initQuadrature(material.quadrature)
integrator.initMaterial(mesh, material)
self.integrators.append(integrator)
+
+ self._info.log("Creating fields and matrices.")
+ # ADD STUFF HERE
+
+ self._info.log("Integrating Jacobian of operator.")
+ #for integrator in integrators:
+ # integrator.integrateJacobian(jacobian, dispT, coords)
return
@@ -99,7 +106,12 @@
"""
Advance to next time step.
"""
- self._info.log("WARNING: Explicit::step() not implemented.")
+ self._info.log("Integrating constant term in operator.")
+ #for integrator in self.integrators:
+ # integrator.integrateConstant(constant, dispT, dispTmdt, coords)
+
+ self._info.log("Solving equations.")
+ # solve
return
More information about the cig-commits
mailing list