[cig-commits] r7030 - in short/3D/PyLith/trunk: libsrc/feassemble
pylith/problems unittests/libtests/feassemble
unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Fri Jun 1 09:38:34 PDT 2007
Author: brad
Date: 2007-06-01 09:38:34 -0700 (Fri, 01 Jun 2007)
New Revision: 7030
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/pylith/problems/Implicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DLinear.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DQuadratic.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Solution2DLinear.py
Log:
Updated implicit time integration formulation, removing superfluous dispT field. Updated corresponding unit tests.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -86,8 +86,9 @@
const ALE::Obj<real_section_type>& coordinates =
mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
- assert(!dispT.isNull());
+ const ALE::Obj<real_section_type>& dispTBctpdt =
+ fields->getReal("dispTBctpdt");
+ assert(!dispTBctpdt.isNull());
// Get cell geometry information that doesn't depend on cell
const int numQuadPts = _quadrature->numQuadPts();
@@ -100,7 +101,7 @@
// Allocate vector for cell values.
_initCellVector();
const int cellVecSize = numBasis*spaceDim;
- double_array dispTCell(cellVecSize);
+ double_array dispTBctpdtCell(cellVecSize);
//double_array gravCell(cellVecSize);
// Allocate vector for total strain
@@ -133,7 +134,7 @@
_resetCellVector();
// Restrict input fields to cell
- mesh->restrict(dispT, *c_iter, &dispTCell[0], cellVecSize);
+ mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -174,7 +175,7 @@
if (1 == cellDim) {
// Compute total strains and then use these to compute stresses
Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -194,7 +195,7 @@
} else if (2 == cellDim) {
// Compute total strains and then use these to compute stresses
Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -220,7 +221,7 @@
} else if (3 == cellDim) {
// Compute total strains and then use these to compute stresses
Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
const std::vector<double_array>& stress =
_material->calcStress(totalStrain);
@@ -286,8 +287,9 @@
const ALE::Obj<real_section_type>& coordinates =
mesh->getRealSection("coordinates");
assert(!coordinates.isNull());
- const ALE::Obj<real_section_type>& dispT = fields->getHistoryItem(1);
- assert(!dispT.isNull());
+ const ALE::Obj<real_section_type>& dispTBctpdt =
+ fields->getReal("dispTBctpdt");
+ assert(!dispTBctpdt.isNull());
// Get parameters used in integration.
const double dt = _dt;
@@ -306,7 +308,7 @@
// Allocate matrix and vectors for cell values.
_initCellMatrix();
const int cellVecSize = numBasis*spaceDim;
- double_array dispTCell(cellVecSize);
+ double_array dispTBctpdtCell(cellVecSize);
// Allocate vector for total strain
int tensorSize = 0;
@@ -338,7 +340,7 @@
_resetCellMatrix();
// Restrict input fields to cell
- mesh->restrict(dispT, *c_iter, &dispTCell[0], cellVecSize);
+ mesh->restrict(dispTBctpdt, *c_iter, &dispTBctpdtCell[0], cellVecSize);
// Get cell geometry information that depends on cell
const double_array& basis = _quadrature->basis();
@@ -348,7 +350,7 @@
if (1 == cellDim) { // 1-D case
// Compute strains
Elasticity::calcTotalStrain1D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const std::vector<double_array>& elasticConsts =
@@ -375,7 +377,7 @@
} else if (2 == cellDim) { // 2-D case
// Compute strains
Elasticity::calcTotalStrain2D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const std::vector<double_array>& elasticConsts =
@@ -428,7 +430,7 @@
} else if (3 == cellDim) { // 3-D case
// Compute strains
Elasticity::calcTotalStrain3D(&totalStrain, basisDeriv,
- dispTCell, numBasis);
+ dispTBctpdtCell, numBasis);
// Get "elasticity" matrix at quadrature points for this cell
const std::vector<double_array>& elasticConsts =
_material->calcDerivElastic(totalStrain);
@@ -533,8 +535,8 @@
// Assemble cell contribution into field. Not sure if this is correct for
// global stiffness matrix.
const ALE::Obj<Mesh::order_type>& globalOrder =
- mesh->getFactory()->getGlobalOrder(mesh, "default", dispT);
- err = updateOperator(*mat, mesh, dispT, globalOrder,
+ mesh->getFactory()->getGlobalOrder(mesh, "default", dispTBctpdt);
+ err = updateOperator(*mat, mesh, dispTBctpdt, globalOrder,
*c_iter, _cellMatrix, ADD_VALUES);
if (err)
throw std::runtime_error("Update to PETSc Mat failed.");
Modified: short/3D/PyLith/trunk/pylith/problems/Implicit.py
===================================================================
--- short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/pylith/problems/Implicit.py 2007-06-01 16:38:34 UTC (rev 7030)
@@ -88,19 +88,17 @@
interfaceConditions, dimension, dt)
self._info.log("Creating fields and matrices.")
- self.fields.addReal("dispT")
self.fields.addReal("dispTBctpdt")
self.fields.addReal("dispIncr")
self.fields.addReal("residual")
- self.fields.createHistory(["dispTBctpdt", "dispT"])
- self.fields.setFiberDimension("dispT", dimension)
+ self.fields.setFiberDimension("dispTBctpdt", dimension)
for constraint in self.constraints:
- constraint.setConstraintSizes(self.fields.getReal("dispT"))
- self.fields.allocate("dispT")
+ constraint.setConstraintSizes(self.fields.getReal("dispTBctpdt"))
+ self.fields.allocate("dispTBctpdt")
for constraint in self.constraints:
- constraint.setConstraints(self.fields.getReal("dispT"))
- self.fields.copyLayout("dispT")
- self.jacobian = mesh.createMatrix(self.fields.getReal("dispT"))
+ constraint.setConstraints(self.fields.getReal("dispTBctpdt"))
+ self.fields.copyLayout("dispTBctpdt")
+ self.jacobian = mesh.createMatrix(self.fields.getReal("dispTBctpdt"))
from pyre.units.time import s
self._solveElastic(mesh, materials, t=0.0*s, dt=dt)
@@ -167,22 +165,19 @@
"""
Hook for doing stuff after advancing time step.
"""
- # This should give us the total displacements after stepping
- # forward from t to t+dt, which is renamed as time step t for the
- # next solve. The field dispTBctpdt contains the displacements
- # from time step t along with the displacement BC from time step
- # t+dt. The displacement increments computed from the residual
- # are then added to this to give us the total displacement field
- # at time t+dt.
+ # After solving, dispTBctpdt contains the displacements at time t
+ # for unconstrained DOF and displacements at time t+dt at
+ # constrained DOF. We add in the displacement increments (only
+ # nonzero at unconstrained DOF) so that after poststrp(),
+ # dispTBctpdt constains the solution at time t+dt.
import pylith.topology.topology as bindings
- dispT = self.fields.getReal("dispT")
dispTBctpdt = self.fields.getReal("dispTBctpdt")
- bindings.addRealSections(dispT, dispTBctpdt,
+ bindings.addRealSections(dispTBctpdt, dispTBctpdt,
self.fields.getReal("dispIncr"))
self._info.log("Updating integrators states.")
for integrator in self.integrators:
- integrator.updateState(dispT)
+ integrator.updateState(dispTBctpdt)
return
@@ -230,13 +225,12 @@
self.solver.solve(dispIncr, self.jacobian, residual)
import pylith.topology.topology as bindings
- dispT = self.fields.getReal("dispT")
dispTBctpdt = self.fields.getReal("dispTBctpdt")
- bindings.addRealSections(dispT, dispTBctpdt, dispIncr)
+ bindings.addRealSections(dispTBctpdt, dispTBctpdt, dispIncr)
self._info.log("Updating integrators states.")
for integrator in self.integrators:
- integrator.updateState(dispT)
+ integrator.updateState(dispTBctpdt)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -95,9 +95,10 @@
topology::FieldsManager fields(mesh);
_initialize(&mesh, &integrator, &fields);
- const ALE::Obj<real_section_type>& dispT = fields.getReal("dispT");
- CPPUNIT_ASSERT(!dispT.isNull());
- integrator.updateState(dispT, mesh);
+ const ALE::Obj<real_section_type>& dispTBctpdt =
+ fields.getReal("dispTBctpdt");
+ CPPUNIT_ASSERT(!dispTBctpdt.isNull());
+ integrator.updateState(dispTBctpdt, mesh);
} // testUpdateState
// ----------------------------------------------------------------------
@@ -143,11 +144,12 @@
topology::FieldsManager fields(mesh);
_initialize(&mesh, &integrator, &fields);
- const ALE::Obj<pylith::real_section_type>& dispT = fields.getReal("dispT");
- CPPUNIT_ASSERT(!dispT.isNull());
+ const ALE::Obj<pylith::real_section_type>& dispTBctpdt =
+ fields.getReal("dispTBctpdt");
+ CPPUNIT_ASSERT(!dispTBctpdt.isNull());
PetscMat jacobian;
- PetscErrorCode err = MeshCreateMatrix(mesh, dispT, MATMPIBAIJ, &jacobian);
+ PetscErrorCode err = MeshCreateMatrix(mesh, dispTBctpdt, MATMPIBAIJ, &jacobian);
CPPUNIT_ASSERT(0 == err);
integrator.integrateJacobian(&jacobian, &fields, mesh);
@@ -257,11 +259,7 @@
// Setup fields
CPPUNIT_ASSERT(0 != fields);
fields->addReal("residual");
- fields->addReal("dispBCTpdt");
- fields->addReal("dispT");
- const char* history[] = { "dispBCTpdt", "dispT" };
- const int historySize = 2;
- fields->createHistory(history, historySize);
+ fields->addReal("dispTBctpdt");
const ALE::Obj<real_section_type>& residual = fields->getReal("residual");
CPPUNIT_ASSERT(!residual.isNull());
@@ -271,17 +269,13 @@
fields->copyLayout("residual");
const int fieldSize = _data->spaceDim * _data->numVertices;
- const ALE::Obj<real_section_type>& dispBCTpdt =
- fields->getReal("dispBCTpdt");
- const ALE::Obj<real_section_type>& dispT = fields->getReal("dispT");
- CPPUNIT_ASSERT(!dispBCTpdt.isNull());
- CPPUNIT_ASSERT(!dispT.isNull());
+ const ALE::Obj<real_section_type>& dispTBctpdt =
+ fields->getReal("dispTBctpdt");
+ CPPUNIT_ASSERT(!dispTBctpdt.isNull());
const int offset = _data->numCells;
for (int iVertex=0; iVertex < _data->numVertices; ++iVertex) {
- dispBCTpdt->updatePoint(iVertex+offset,
- &_data->fieldTpdt[iVertex*_data->spaceDim]);
- dispT->updatePoint(iVertex+offset,
- &_data->fieldT[iVertex*_data->spaceDim]);
+ dispTBctpdt->updatePoint(iVertex+offset,
+ &_data->fieldTpdt[iVertex*_data->spaceDim]);
} // for
} // _initialize
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicit.py 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicit.py 2007-06-01 16:38:34 UTC (rev 7030)
@@ -48,7 +48,7 @@
"""
K = self._calculateStiffnessMat()
- self.valsResidual = -numpy.dot(K, self.fieldT)
+ self.valsResidual = -numpy.dot(K, self.fieldTpdt)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DLinear.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DLinear.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -80,8 +80,8 @@
};
const double pylith::feassemble::ElasticityImplicitData1DLinear::_valsResidual[] = {
- 2.02500000e+10,
- -2.02500000e+10,
+ 2.53125000e+10,
+ -2.53125000e+10,
};
const double pylith::feassemble::ElasticityImplicitData1DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DQuadratic.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData1DQuadratic.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -93,9 +93,9 @@
};
const double pylith::feassemble::ElasticityImplicitData1DQuadratic::_valsResidual[] = {
- -1.41750000e+11,
- 3.24000000e+11,
- -1.82250000e+11,
+ -1.70437500e+11,
+ 3.91500000e+11,
+ -2.21062500e+11,
};
const double pylith::feassemble::ElasticityImplicitData1DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DLinear.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DLinear.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -66,7 +66,7 @@
};
const double pylith::feassemble::ElasticityImplicitData2DLinear::_fieldTpdt[] = {
- 1.90000000e+00, -9.00000000e-01,
+ 1.30000000e+00, -9.00000000e-01,
1.40000000e+00, 1.50000000e+00,
5.00000000e-01, -9.00000000e-01,
};
@@ -84,9 +84,9 @@
};
const double pylith::feassemble::ElasticityImplicitData2DLinear::_valsResidual[] = {
- -4.20750000e+10, -2.88750000e+10,
- 3.96000000e+10, 2.47500000e+09,
- 2.47500000e+09, 2.64000000e+10,
+ 1.81500000e+10, 1.48500000e+10,
+ -4.95000000e+09, -1.32000000e+10,
+ -1.32000000e+10, -1.65000000e+09,
};
const double pylith::feassemble::ElasticityImplicitData2DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DQuadratic.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData2DQuadratic.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -121,12 +121,12 @@
};
const double pylith::feassemble::ElasticityImplicitData2DQuadratic::_valsResidual[] = {
- -2.15074074e+10, 4.58148148e+09,
- -1.24762963e+11, -3.03740741e+10,
- -2.17574074e+10, -9.72851852e+10,
- 2.25390741e+11, 3.72085185e+11,
- -3.39309259e+11, -4.12214815e+11,
- 2.81946296e+11, 1.63207407e+11,
+ -1.48468519e+11, -1.75705556e+11,
+ -2.40611111e+11, -4.21759259e+10,
+ -1.05142593e+11, -3.85807407e+11,
+ 5.56192593e+11, 7.79718519e+11,
+ -7.80303704e+11, -9.67622222e+11,
+ 7.18333333e+11, 7.91592593e+11,
};
const double pylith::feassemble::ElasticityImplicitData2DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DLinear.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DLinear.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DLinear.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -89,10 +89,10 @@
};
const double pylith::feassemble::ElasticityImplicitData3DLinear::_valsResidual[] = {
- -4.51920000e+10, 1.85610000e+10, 1.04910000e+10,
- 2.74380000e+10, 8.07000000e+09, 9.68400000e+09,
- 8.07000000e+09, -2.09820000e+10, -5.64900000e+09,
- 9.68400000e+09, -5.64900000e+09, -1.45260000e+10,
+ -3.30870000e+10, 1.12980000e+10, 2.42100000e+09,
+ 1.93680000e+10, 5.64900000e+09, 8.07000000e+09,
+ 5.64900000e+09, -1.29120000e+10, -4.03500000e+09,
+ 8.07000000e+09, -4.03500000e+09, -6.45600000e+09,
};
const double pylith::feassemble::ElasticityImplicitData3DLinear::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DQuadratic.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DQuadratic.cc 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityImplicitData3DQuadratic.cc 2007-06-01 16:38:34 UTC (rev 7030)
@@ -164,16 +164,16 @@
};
const double pylith::feassemble::ElasticityImplicitData3DQuadratic::_valsResidual[] = {
- 1.37302827e+12, 3.59474617e+11, 2.01771825e+12,
- -1.31076938e+11, -1.24488420e+11, 7.73353086e+09,
- -3.68697654e+11, -2.72246232e+12, -1.48102395e+11,
- 5.48027062e+11, 3.15720222e+11, 2.06778262e+12,
- -1.21553814e+12, -2.00599957e+12, -3.41367532e+12,
- -1.41822963e+11, 2.78297909e+12, -2.96592593e+11,
- -6.77812870e+12, -6.14554194e+12, -4.48612599e+12,
- 4.50985444e+11, 1.95704638e+12, 3.90063222e+11,
- 9.75448605e+11, 4.73812614e+12, 2.71904605e+11,
- 5.28777501e+12, 8.45145802e+11, 3.58929407e+12,
+ -2.82689704e+11, -1.77165531e+11, 4.81077062e+11,
+ -1.89584198e+11, -2.62414988e+11, -2.05662222e+10,
+ -2.64410790e+11, -2.99364405e+12, -1.43747136e+11,
+ 4.64419852e+11, 9.60363704e+10, 1.15315605e+12,
+ 8.39069593e+11, -1.89835947e+12, -3.19266342e+12,
+ 7.92577481e+11, 3.50778306e+12, 1.54368667e+12,
+ -4.45565456e+12, -5.73392127e+12, -1.56481211e+12,
+ -4.92241198e+11, 2.04747312e+12, 1.67797247e+11,
+ 5.70022407e+11, 4.69695932e+12, -1.55162542e+12,
+ 3.01849111e+12, 7.17253432e+11, 3.12769728e+12,
};
const double pylith::feassemble::ElasticityImplicitData3DQuadratic::_valsJacobian[] = {
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Solution2DLinear.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Solution2DLinear.py 2007-06-01 15:43:29 UTC (rev 7029)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/Solution2DLinear.py 2007-06-01 16:38:34 UTC (rev 7030)
@@ -38,7 +38,7 @@
# Input fields
self.dt = 0.01
- self.fieldTpdt = numpy.array([ 1.9, -0.9,
+ self.fieldTpdt = numpy.array([ 1.3, -0.9,
1.4, 1.5,
0.5, -0.9 ], dtype=numpy.float64)
self.fieldT = numpy.array([ 1.6, -0.8,
More information about the cig-commits
mailing list