[cig-commits] r6762 - short/3D/PyLith/trunk/libsrc/feassemble
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed May 2 20:23:45 PDT 2007
Author: willic3
Date: 2007-05-02 20:23:45 -0700 (Wed, 02 May 2007)
New Revision: 6762
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh
Log:
Fixed some compile problems. More need fixing.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-05-03 02:56:55 UTC (rev 6761)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ImplicitElasticity.cc 2007-05-03 03:23:45 UTC (rev 6762)
@@ -15,6 +15,7 @@
#include "ImplicitElasticity.hh" // implementation of class methods
#include "Quadrature.hh" // USES Quadrature
+#include "IntegratorElasticity.hh" // USES IntegratorElasticity
#include "pylith/materials/ElasticMaterial.hh" // USES ElasticMaterial
#include "pylith/utils/array.hh" // USES double_array
@@ -267,7 +268,7 @@
_cellVector[iBlock+1] -= N1*s12 + N2*s22;
} // for
} // for
- err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(8+2+9)));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
@@ -298,7 +299,7 @@
_cellVector[iBlock+2] -= N1*s13 + N2*s23 + N3*s33;
} // for
} // for
- err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
+ PetscErrorCode err = PetscLogFlops(numQuadPts*(1+numBasis*(3+12)));
if (err)
throw std::runtime_error("Logging PETSc flops failed.");
} // if/else
@@ -324,7 +325,7 @@
assert(!mesh.isNull());
// Clear stiffness matrix. Not sure if this is the correct way.
- PetscErrorCode err = MatZeroEntries(mat);
+ PetscErrorCode err = MatZeroEntries(*mat);
// Get information about section
const ALE::Obj<Mesh::label_sequence>& cells = mesh->heightStratum(0);
@@ -342,12 +343,12 @@
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
const double_array& quadWts = _quadrature->quadWts();
- const int numBasis = _quadrature->numCorners();
+ const int numBasis = _quadrature->numBasis();
const int spaceDim = _quadrature->spaceDim();
const int cellDim = _quadrature->cellDim();
if (cellDim != spaceDim)
- throw std::logic_error("Not implemented yet.")
+ throw std::logic_error("Not implemented yet.");
// Allocate vector for cell values (if necessary)
_initCellMatrix();
@@ -381,6 +382,10 @@
// Reset element vector to zero
_resetCellMatrix();
+ // Restrict input fields to cell
+ const real_section_type::value_type* dispTCell =
+ mesh->restrict(dispT, *cellIter);
+
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basisQuad();
const double_array& basisDeriv = _quadrature->basisDerivQuad();
@@ -388,14 +393,14 @@
// 1D Case
if (1 == cellDim) {
- assert(1 == numElasticConsts);
// Compute strains
IntegratorElasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
- const int numElasticConsts = _material->numElasticConsts();
+ const int numElasticConsts = _material->_numElasticConsts();
+ assert(1 == numElasticConsts);
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -418,14 +423,14 @@
// 2D Case
} else if (2 == cellDim) {
- assert(6 == numElasticConsts);
// Compute strains
IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
- const int numElasticConsts = _material->numElasticConsts();
+ const int numElasticConsts = _material->_numElasticConsts();
+ assert(6 == numElasticConsts);
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -467,14 +472,14 @@
// 3D Case
} else if (3 == cellDim) {
- assert(21 == numElasticConsts);
// Compute strains
IntegratorElasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
dispTCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const double_array& elasticConsts = _material->calcDerivElastic(&totalStrain);
- const int numElasticConsts = _material->numElasticConsts();
+ const int numElasticConsts = _material->_numElasticConsts();
+ assert(21 == numElasticConsts);
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -535,15 +540,15 @@
C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
- _cellMatrix[iblock +jBlock ] += ki0j0;
- _cellMatrix[iblock+1+jBlock ] += ki0j1;
- _cellMatrix[iblock+2+jBlock ] += ki0j2;
- _cellMatrix[iblock +jBlock+1] += ki0j1;
- _cellMatrix[iblock+1+jBlock+1] += ki1j1;
- _cellMatrix[iblock+2+jBlock+1] += ki1j2;
- _cellMatrix[iblock +jBlock+2] += ki0j2;
- _cellMatrix[iblock+1+jBlock+2] += ki1j2;
- _cellMatrix[iblock+2+jBlock+2] += ki2j2;
+ _cellMatrix[iBlock +jBlock ] += ki0j0;
+ _cellMatrix[iBlock+1+jBlock ] += ki0j1;
+ _cellMatrix[iBlock+2+jBlock ] += ki0j2;
+ _cellMatrix[iBlock +jBlock+1] += ki0j1;
+ _cellMatrix[iBlock+1+jBlock+1] += ki1j1;
+ _cellMatrix[iBlock+2+jBlock+1] += ki1j2;
+ _cellMatrix[iBlock +jBlock+2] += ki0j2;
+ _cellMatrix[iBlock+1+jBlock+2] += ki1j2;
+ _cellMatrix[iBlock+2+jBlock+2] += ki2j2;
} // for
} // for
} // for
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc 2007-05-03 02:56:55 UTC (rev 6761)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.cc 2007-05-03 03:23:45 UTC (rev 6762)
@@ -20,8 +20,8 @@
// Constructor
pylith::feassemble::IntegratorImplicit::IntegratorImplicit(void) :
Integrator(),
- _dt(0.0),
- _dtm1(0.0)
+ _dt(-1.0),
+ _dtm1(-1.0)
{ // constructor
} // constructor
@@ -45,7 +45,10 @@
void
pylith::feassemble::IntegratorImplicit::timeStep(const double dt)
{ // timeStep
- _dtm1 = _dt;
+ if (_dt != -1.0)
+ _dtm1 = _dt;
+ else
+ _dtm1 = _dt;
_dt = dt;
assert(_dt == _dtm1); // For now, don't allow variable time step
} // timeStep
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh 2007-05-03 02:56:55 UTC (rev 6761)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorImplicit.hh 2007-05-03 03:23:45 UTC (rev 6762)
@@ -85,7 +85,7 @@
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;
+ const ALE::Obj<Mesh>& mesh) = 0;
/** Compute Jacobian matrix (A) associated with operator.
*
@@ -96,7 +96,7 @@
virtual
void integrateJacobian(PetscMat* mat,
const ALE::Obj<real_section_type>& fieldIn,
- const ALE::Obj<real_section_type>& coordinates) = 0;
+ const ALE::Obj<Mesh>& mesh) = 0;
// PROTECTED METHODS ////////////////////////////////////////////////////
protected :
More information about the cig-commits
mailing list