[cig-commits] r7027 - in short/3D/PyLith/trunk: libsrc/feassemble
unittests/libtests/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu May 31 21:13:56 PDT 2007
Author: brad
Date: 2007-05-31 21:13:55 -0700 (Thu, 31 May 2007)
New Revision: 7027
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
Log:
Fixed forming Jacobian in 3-D for implicit elasticity integrator. Implemented Python unit tests for implicit elasticity integrator.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 03:24:33 UTC (rev 7026)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2007-06-01 04:13:55 UTC (rev 7027)
@@ -479,6 +479,10 @@
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 ki1j0 =
+ C1122 * Niq * Njp + C1112 * Nip * Njp + C1123 * Nir * Njp +
+ C2212 * Niq * Njq + C1212 * Nip * Njq + C1223 * Nir * Njq +
+ C2213 * Niq * Njr + C1213 * Nip * Njr + C2313 * Nir * Njr;
const double ki1j1 =
C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
@@ -487,11 +491,18 @@
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 ki2j0 =
+ C1133 * Nir * Njp + C1123 * Niq * Njp + C1113 * Nip * Njp +
+ C3312 * Nir * Njq + C1223 * Niq * Njq + C1213 * Nip * Njq +
+ C3313 * Nir * Njr + C2313 * Niq * Njr + C1313 * Nip * Njr;
+ const double ki2j1 =
+ C2233 * Nir * Njq + C2223 * Niq * Njq + C2213 * Nip * Njq +
+ C3312 * Nir * Njp + C1223 * Niq * Njp + C1213 * Nip * Njp +
+ C3323 * Nir * Njr + C2323 * Niq * Njr + C2313 * Nip * Njr;
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;
-
const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
const int iBlock2 = (iBasis*spaceDim+2) * (numBasis*spaceDim);
@@ -501,11 +512,11 @@
_cellMatrix[iBlock +jBlock ] += ki0j0;
_cellMatrix[iBlock +jBlock1] += ki0j1;
_cellMatrix[iBlock +jBlock2] += ki0j2;
- _cellMatrix[iBlock1+jBlock ] += ki0j1;
+ _cellMatrix[iBlock1+jBlock ] += ki1j0;
_cellMatrix[iBlock1+jBlock1] += ki1j1;
_cellMatrix[iBlock1+jBlock2] += ki1j2;
- _cellMatrix[iBlock2+jBlock ] += ki0j2;
- _cellMatrix[iBlock2+jBlock1] += ki1j2;
+ _cellMatrix[iBlock2+jBlock ] += ki2j0;
+ _cellMatrix[iBlock2+jBlock1] += ki2j1;
_cellMatrix[iBlock2+jBlock2] += ki2j2;
} // for
} // for
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2007-06-01 03:24:33 UTC (rev 7026)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestElasticityImplicit.cc 2007-06-01 04:13:55 UTC (rev 7027)
@@ -166,8 +166,6 @@
CPPUNIT_ASSERT_EQUAL(nrowsE, nrows);
CPPUNIT_ASSERT_EQUAL(ncolsE, ncols);
- MatView(jacobian, PETSC_VIEWER_STDOUT_WORLD);
-
PetscMat jDense;
PetscMat jSparseAIJ;
MatConvert(jacobian, MATSEQAIJ, MAT_INITIAL_MATRIX, &jSparseAIJ);
Modified: short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py
===================================================================
--- short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py 2007-06-01 03:24:33 UTC (rev 7026)
+++ short/3D/PyLith/trunk/unittests/pytests/feassemble/TestElasticityImplicit.py 2007-06-01 04:13:55 UTC (rev 7027)
@@ -17,6 +17,8 @@
import unittest
from pylith.feassemble.ElasticityImplicit import ElasticityImplicit
+from spatialdata.geocoords.CSCart import CSCart
+
# ----------------------------------------------------------------------
class TestElasticityImplicit(unittest.TestCase):
"""
@@ -43,4 +45,185 @@
return
+ def test_needNewJacobian(self):
+ """
+ Test needNewJacobian().
+ """
+ integrator = ElasticityImplicit()
+ self.assertEqual(False, integrator.needNewJacobian())
+ return
+
+
+ def test_setMesh(self):
+ """
+ Test setMesh().
+ """
+ cs = CSCart()
+ cs.spaceDim = 2
+
+ from pylith.meshio.MeshIOAscii import MeshIOAscii
+ importer = MeshIOAscii()
+ importer.filename = "data/tri3.mesh"
+ importer.coordsys = cs
+ mesh = importer.read(debug=False, interpolate=False)
+
+ integrator = ElasticityImplicit()
+ integrator.setMesh(mesh)
+ self.assertEqual(mesh, integrator.mesh)
+ return
+
+
+ def test_timeStep(self):
+ """
+ Test timeStep().
+ """
+ from pyre.units.time import second
+ dt = 2.4*second
+ integrator = ElasticityImplicit()
+ integrator.timeStep = dt
+ self.assertEqual(dt, integrator.timeStep)
+ return
+
+
+ def test_stableTimeStep(self):
+ """
+ Test stableTimeStep().
+ """
+ from pyre.units.time import second
+ dt = 2.3*second
+ integrator = ElasticityImplicit()
+ integrator.timeStep(dt)
+ self.assertEqual(dt, integrator.stableTimeStep())
+ return
+
+
+ def test_needNewJacobian(self):
+ """
+ Test needNewJacobian().
+ """
+ integrator = ElasticityImplicit()
+ self.assertEqual(False, integrator.needNewJacobian())
+ return
+
+
+ def test_integrateResidual(self):
+ """
+ Test integrateResidual().
+
+ WARNING: This is not a rigorous test of integrateResidual() because we
+ neither set the input fields or verify the results.
+ """
+ (mesh, integrator, fields) = self._initialize()
+
+ residual = fields.getReal("residual")
+ integrator.integrateResidual(residual, fields)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_integrateJacobian(self):
+ """
+ Test integrateJacobian().
+
+ WARNING: This is not a rigorous test of integrateJacobian() because we
+ neither set the input fields or verify the results.
+ """
+ (mesh, integrator, fields) = self._initialize()
+
+ jacobian = mesh.createMatrix(fields.getReal("residual"))
+ import pylith.utils.petsc as petsc
+ petsc.mat_setzero(jacobian)
+ integrator.integrateJacobian(jacobian, fields)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ def test_updateState(self):
+ """
+ Test updateState().
+
+ WARNING: This is not a rigorous test of updateState() because we
+ neither set the input fields or verify the results.
+ """
+ (mesh, integrator, fields) = self._initialize()
+
+ dispT = fields.getReal("dispT")
+ integrator.updateState(dispT)
+
+ # We should really add something here to check to make sure things
+ # actually initialized correctly
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _initialize(self):
+ """
+ Initialize integrator.
+ """
+ from pyre.units.time import second
+ dt = 1.0*second
+
+ # Setup mesh
+ cs = CSCart()
+ cs.spaceDim = 2
+ from pylith.meshio.MeshIOAscii import MeshIOAscii
+ importer = MeshIOAscii()
+ importer.filename = "data/tri3.mesh"
+ importer.coordsys = cs
+ mesh = importer.read(debug=False, interpolate=False)
+
+ # Setup material
+ from pylith.feassemble.FIATSimplex import FIATSimplex
+ cell = FIATSimplex()
+ cell.shape = "triangle"
+ cell.degree = 1
+ cell.order = 1
+ from pylith.feassemble.quadrature.Quadrature2D import Quadrature2D
+ quadrature = Quadrature2D()
+ quadrature.cell = cell
+
+ from spatialdata.spatialdb.SimpleDB import SimpleDB
+ from spatialdata.spatialdb.SimpleIOAscii import SimpleIOAscii
+ iohandler = SimpleIOAscii()
+ iohandler.filename = "data/elasticplanestrain.spatialdb"
+ db = SimpleDB()
+ db.label = "elastic plane strain"
+ db.iohandler = iohandler
+
+ from pylith.materials.ElasticPlaneStrain import ElasticPlaneStrain
+ material = ElasticPlaneStrain()
+ material.id = 0
+ material.label = "elastic plane strain"
+ material.db = db
+ material.quadrature = quadrature
+
+ # Setup integrator
+ integrator = ElasticityImplicit()
+ integrator.setMesh(mesh)
+ integrator.initQuadrature(material.quadrature)
+ integrator.initMaterial(mesh, material)
+ integrator.timeStep(dt)
+
+ # Setup fields
+ from pylith.topology.FieldsManager import FieldsManager
+ fields = FieldsManager(mesh)
+ fields.addReal("residual")
+ fields.addReal("dispBCTpdt")
+ fields.addReal("dispT")
+ fields.createHistory(["dispBCTpdt", "dispT"])
+ fields.setFiberDimension("residual", cs.spaceDim)
+ fields.allocate("residual")
+ fields.copyLayout("residual")
+
+ import pylith.topology.topology as bindings
+ bindings.zeroRealSection(fields.getReal("residual"))
+
+ return (mesh, integrator, fields)
+
+
# End of file
More information about the cig-commits
mailing list