[cig-commits] r6984 -
short/3D/PyLith/trunk/unittests/libtests/feassemble/data
brad at geodynamics.org
brad at geodynamics.org
Mon May 28 16:27:41 PDT 2007
Author: brad
Date: 2007-05-28 16:27:41 -0700 (Mon, 28 May 2007)
New Revision: 6984
Added:
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py
Removed:
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py
Modified:
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh
Log:
Started work implementing C++ unit tests for elasticity integrators.
Copied: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py (from rev 6980, short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py)
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py 2007-05-27 11:25:19 UTC (rev 6980)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/ElasticityExplicit.py 2007-05-28 23:27:41 UTC (rev 6984)
@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/ElasticityExplicit.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ ElasticityExplicit object.
+
+from IntegratorElasticity import IntegratorElasticity
+
+import numpy
+
+# ----------------------------------------------------------------------
+
+# ElasticityExplicit class
+class ElasticityExplicit(IntegratorElasticity):
+ """
+ Python application for generating C++ data files for testing C++
+ ElasticityExplicit object.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="elasticityexplicit"):
+ """
+ Constructor.
+ """
+ IntegratorElasticity.__init__(self, name)
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _calculateResidual(self):
+ """
+ Calculate contribution to residual of operator for integrator.
+
+ {r} = (1/dt**2)[M](-{u(t+dt)} + 2 {u(t)} - {u(t-dt)}) -
+ [K]{u(t)}
+ """
+ dispResult = -self.fieldTpdt + 2.0*self.fieldT - self.fieldTmdt
+ self.valsResidual = 1.0/self.dt**2 * numpy.dot(self.M, dispResult) - \
+ numpy.dot(self.K, fieldT)
+ return
+
+
+ def _calculateJacobian(self):
+ """
+ Calculate contribution to Jacobian matrix of operator for integrator.
+
+ [A] = (1/dt**2)[M]
+ """
+ self.valsJacobian = 1.0/self.dt**2 * self.M
+ return
+
+
+# End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2007-05-28 23:26:31 UTC (rev 6983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorApp.py 2007-05-28 23:27:41 UTC (rev 6984)
@@ -55,24 +55,34 @@
"""
Script.__init__(self, name)
- self.numVertices = None
+ # Mesh information
+ self.meshFilename = None
+
+ # Quadrature information
self.spaceDim = None
- self.numCells = None
self.cellDim = None
self.numBasis = None
self.numQuadPts = None
- self.fiberDim = None
-
self.quadPts = None
self.quadWts = None
- self.vertices = None
- self.cells = None
-
self.basis = None
self.basisDeriv = None
- self.fieldIn = None
- self.valsActions = None
- self.valsMatrix = None
+
+ # Material information
+ self.matType = None
+ self.matDBFilename = None
+ self.matId = None
+ self.matLabel = None
+
+ # Input fields
+ self.dt = None
+ self.fieldTpdt = None
+ self.fieldT = None
+ self.fieldTmdt = None
+
+ # Calculated values
+ self.valsResidual = None
+ self.valsJacobian = None
return
@@ -81,8 +91,8 @@
Run the application.
"""
self._initialize()
- self._calculateMatrix()
- self._calculateAction()
+ self._calculateResidual()
+ self._calculateJacobian()
self._initData()
self.data.write(self.name)
return
@@ -117,48 +127,37 @@
return
- def _calculateMatrix(self):
+ def _calculateResidual(self):
"""
- Calculate matrix associated with integration.
+ Calculate contribution to residual of operator for integrator.
"""
raise NotImplementedError("Not implemented in abstract base class.")
- def _calculateAction(self):
+ def _calculateJacobian(self):
"""
- Calculate integration action using matrix.
+ Calculate contribution to Jacobian matrix of operator for integrator.
"""
self.valsAction = numpy.dot(self.valsMatrix, self.fieldIn)[:]
return
def _initData(self):
- self.data.addScalar(vtype="int", name="_numVertices",
- value=self.numVertices,
- format="%d")
+ # Mesh information
+ self.data.addScalar(vtype="char*", name="_meshFilename",
+ value=self.meshFilename,
+ format="%s")
+
+ # Quadrature information
self.data.addScalar(vtype="int", name="_spaceDim", value=self.spaceDim,
format="%d")
- self.data.addScalar(vtype="int", name="_numCells", value=self.numCells,
- format="%d")
self.data.addScalar(vtype="int", name="_cellDim", value=self.cellDim,
format="%d")
- self.data.addScalar(vtype="int", name="_numBasis", value=
- self.numBasis,
+ self.data.addScalar(vtype="int", name="_numBasis", value=self.numBasis,
format="%d")
- self.data.addScalar(vtype="int", name="_numQuadPts",
- value=self.numQuadPts,
+ self.data.addScalar(vtype="int", name="_numQuadPts", value=self.numQuadPts,
format="%d")
- self.data.addScalar(vtype="int", name="_fiberDim",
- value=self.fiberDim,
- format="%d")
-
- self.data.addArray(vtype="double", name="_vertices", values=self.vertices,
- format="%16.8e", ncols=self.spaceDim)
- self.data.addArray(vtype="int", name="_cells", values=self.cells,
- format="%8d", ncols=self.numVertices)
-
- self.data.addArray(vtype="double", name="_quadPts",
- values=self.quadPts,
+ self.data.addArray(vtype="double", name="_quadPts", values=self.quadPts,
format="%16.8e", ncols=self.cellDim)
self.data.addArray(vtype="double", name="_quadWts", values=self.quadWts,
format="%16.8e", ncols=self.numQuadPts)
@@ -168,16 +167,37 @@
self.data.addArray(vtype="double", name="_basisDeriv",
values=self.basisDeriv,
format="%16.8e", ncols=self.cellDim)
-
- self.data.addArray(vtype="double", name="_fieldIn",
- values=self.fieldIn,
+
+ # Material information
+ self.data.addScalar(vtype="char*", name="_matType", value=self.matType,
+ format="%s")
+ self.data.addScalar(vtype="char*", name="_matDBFilename",
+ value=self.matDBFilename, format="%s")
+ self.data.addScalar(vtype="int", name="_matId",
+ value=self.matId, format="%d")
+ self.data.addScalar(vtype="char*", name="_matLabel",
+ value=self.matLabel, format="%s")
+
+ # Input files
+ self.data.addScalar(vtype="double", name="_dt", values=self.dt,
+ format="%16.8e")
+ self.data.addArray(vtype="double", name="_fieldTpdt",
+ values=self.fieldTpdt,
format="%16.8e", ncols=self.spaceDim)
- self.data.addArray(vtype="double", name="_valsAction",
- values=self.valsAction,
+ self.data.addArray(vtype="double", name="_fieldT",
+ values=self.fieldTpdt,
format="%16.8e", ncols=self.spaceDim)
- self.data.addArray(vtype="double", name="_valsMatrix",
- values=self.valsMatrix,
+ self.data.addArray(vtype="double", name="_fieldTmdt",
+ values=self.fieldTpdt,
format="%16.8e", ncols=self.spaceDim)
+
+ # Calculated values
+ self.data.addArray(vtype="double", name="_valsResidual",
+ values=self.valsResidual,
+ format="%16.8e", ncols=self.spaceDim)
+ self.data.addArray(vtype="double", name="_valsJacobian",
+ values=self.valsJacobian,
+ format="%16.8e", ncols=self.spaceDim)
return
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc 2007-05-28 23:26:31 UTC (rev 6983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.cc 2007-05-28 23:27:41 UTC (rev 6984)
@@ -15,23 +15,24 @@
// ----------------------------------------------------------------------
// Constructor
pylith::feassemble::IntegratorData::IntegratorData(void) :
- numVertices(0),
+ meshFilename(0),
spaceDim(0),
- numCells(0),
cellDim(0),
numBasis(0),
numQuadPts(0),
- fiberDim(0),
- vertices(0),
- cells(0),
quadPts(0),
quadWts(0),
basis(0),
basisDeriv(0),
- fieldIn(0),
- valsAction(0),
- valsMatrix(0),
- valsLumped(0)
+ matType(0),
+ matDBFilename(0),
+ matId(0),
+ matlabel(0),
+ fieldTpdt(0),
+ fieldT(0),
+ fieldTmdt(0),
+ valsResidual(0),
+ valsJacobian(0)
{ // constructor
} // constructor
@@ -41,4 +42,5 @@
{ // destructor
} // destructor
+
// End of file
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh 2007-05-28 23:26:31 UTC (rev 6983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorData.hh 2007-05-28 23:27:41 UTC (rev 6984)
@@ -34,35 +34,43 @@
// PUBLIC MEMBERS ///////////////////////////////////////////////////////
public:
- int numVertices; ///< Number of vertices
- int spaceDim; ///< Number of dimensions in vertex coordinates
- int numCells; ///< Number of cells
- int cellDim; ///< Number of dimensions associated with cell
- int numBasis; ///< Number of vertices in cell
- int numQuadPts; ///< Number of quadrature points
- int fiberDim; ///< Number of values per vertex in field
-
/// @name Mesh information
//@{
- double* vertices; ///< Pointer to coordinates of vertices
- int* cells; ///< Pointer to indices of vertices in cells
+ char* meshFilename; ///< Name of mesh file.
//@}
- /// @name Discretization information
+ /// @name Quadrature information
//@{
+ int spaceDim; ///< Number of dimensions in vertex coordinates
+ int cellDim; ///< Number of dimensions associated with cell
+ int numBasis; ///< Number of vertices in cell
+ int numQuadPts; ///< Number of quadrature points
double* quadPts; ///< Coordinates of quad pts in ref cell
double* quadWts; ///< Weights of quadrature points
double* basis; ///< Basis fns at quadrature points
double* basisDeriv; ///< Derivatives of basis fns at quad pts
//@}
- /// @name Integration information
+ /// @name Material information
//@{
- double* fieldIn; ///< Input field for integration action
- double* valsAction; ///< Expected output for integration action
- double* valsMatrix; ///< Expected output for integration
- double* valsLumped; ///< Expected output for lumped integration
+ double* matType; ///< String corresponding to material type.
+ double* matDBFilename; ///< Filename for database of material properties.
+ double* matId; ///< Material identifier.
+ double* matLabel; ///< Label of material.
//@}
+
+ /// @name Input fields
+ //@{
+ double* fieldTpdt; ///< Input field at time t+dt.
+ double* fieldT; ///< Input field at time t.
+ double* fieldTmdt; ///< Input field at time t-dt.
+ //@}
+
+ /// @name Calculated values.
+ //@{
+ double* valsResidual; ///< Expected values from residual calculation.
+ double* valsJacobian; ///< Expected values from Jacobian calculation.
+ //@}
};
#endif // pylith_feassemble_integratordata_hh
Added: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py 2007-05-28 23:26:31 UTC (rev 6983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorElasticity.py 2007-05-28 23:27:41 UTC (rev 6984)
@@ -0,0 +1,129 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+## @file unittests/libtests/feassemble/data/IntegratorElasticity.py
+
+## @brief Python application for generating C++ data files for testing
+## C++ elasticity integrator objects.
+
+from IntegratorApp import IntegratorApp
+
+import numpy
+import feutils
+
+# ----------------------------------------------------------------------
+
+# IntegratorElasticity class
+class IntegratorElasticity(IntegratorApp):
+ """
+ Python application for generating C++ data files for testing C++
+ elasticity integrator objects.
+ """
+
+ # PUBLIC METHODS /////////////////////////////////////////////////////
+
+ def __init__(self, name="integratorelasticity"):
+ """
+ Constructor.
+ """
+ IntegratorApp.__init__(self, name)
+
+ self.density = None
+ self.lameMu = None
+ self.lameLambda = None
+
+ self.K = self._calculateStiffnessMat()
+ self.M = self._calculateMassMat()
+ return
+
+
+ # PRIVATE METHODS ////////////////////////////////////////////////////
+
+ def _calculateStiffnessMat(self):
+ """
+ Calculate stiffness matrix.
+
+ """
+ import feutils
+
+ self.K = numpy.zeros ( (self.spaceDim*self.numVertices,
+ self.spaceDim*self.numVertices),
+ dtype=numpy.float64)
+
+ # Matrix of elasticity values
+ D = self._calculateElasticityMat()
+
+ for cell in self.cells:
+ vertices = self.vertices[cell, :]
+
+ # Matrix of shape function derivatives
+ B = self._calculateBasisDerivMat(vertices)
+
+ cellK = numpy.dot(numpy.dot(B.transpose(), D), B)
+
+ feutils.assembleMat(self.K, cellK, cell, self.spaceDim)
+ return
+
+
+ def _calculateMassMat(self):
+ """
+ Calculate mass matrix.
+ """
+
+ self.M = numpy.zeros ( (self.spaceDim*self.numVertices,
+ self.spaceDim*self.numVertices),
+ dtype=numpy.float64)
+
+
+ for cell in self.cells:
+ vertices = self.vertices[cell, :]
+ (jacobian, jacobianInv, jacobianDet) = \
+ feutils.calculateJacobianQuad(self.quadrature, vertices)
+ for iQuad in xrange(self.numQuadPts):
+ wt = self.quadWts[iQuad] * jacobianDet[iQuad]
+ N = self._calculateBasisMat(vertices, iQuad)
+ cellM[:] += self.density * wt * numpy.dot(N.transpose(), N)
+ feutils.assembleMat(self.M, cellM, cell, self.spaceDim)
+ return
+
+
+ def _calculateElasticityMat(self):
+ """
+ Calculate elasticity matrix.
+ """
+ return
+
+ def _calculateBasisMat(self, iQuad):
+ """
+ Calculate matrix of basis functions.
+ """
+ N = numpy.zeros( (self.spaceDim, self.spaceDim*self.numBasis),
+ dtype=numpy.float64)
+ for iBasis in xrange(self.numBasis):
+ for iDim in xrange(self.spaceDim):
+ N[iDim, iBasis*self.spaceDim+iDim] = self.basis[iQuad, iBasis]
+ return N
+
+
+ def _calculateBasisDerivMat(self):
+ """
+ Calculate matrix of derivatives of basis functions.
+ """
+ B = numpy.zeros( (self.spaceDim, self.spaceDim*self.numBasis),
+ dtype=numpy.float64)
+ for iBasis in xrange(self.numBasis):
+ for iDim in xrange(self.spaceDim):
+ N[iDim, iBasis*self.spaceDim+iDim] = self.basis[iQuad, iBasis]
+ return
+
+
+# End of file
Deleted: short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py 2007-05-28 23:26:31 UTC (rev 6983)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/data/IntegratorInertia.py 2007-05-28 23:27:41 UTC (rev 6984)
@@ -1,109 +0,0 @@
-#!/usr/bin/env python
-#
-# ----------------------------------------------------------------------
-#
-# Brad T. Aagaard
-# U.S. Geological Survey
-#
-# <LicenseText>
-#
-# ----------------------------------------------------------------------
-#
-
-## @file unittests/libtests/feassemble/data/IntegratorInertia.py
-
-## @brief Python application for generating C++ data files for testing
-## C++ IntegratorInertia object.
-
-from IntegratorApp import IntegratorApp
-
-import numpy
-
-# ----------------------------------------------------------------------
-
-# IntegratorInertia class
-class IntegratorInertia(IntegratorApp):
- """
- Python application for generating C++ data files for testing C++
- IntegratorInertia object.
- """
-
- # PUBLIC METHODS /////////////////////////////////////////////////////
-
- def __init__(self, name="integratorinertia"):
- """
- Constructor.
- """
- IntegratorApp.__init__(self, name)
-
- self.valsLumped = None
- return
-
-
- def _initData(self):
- IntegratorApp._initData(self)
- self.data.addArray(vtype="double", name="_valsLumped",
- values=self.valsLumped,
- format="%16.8e", ncols=self.spaceDim)
- return
-
-
- # PRIVATE METHODS ////////////////////////////////////////////////////
-
- def _calculateMatrix(self):
- """
- Calculate matrix associated with integration.
- """
- import feutils
-
- self.valsMatrix = numpy.zeros( (self.fiberDim*self.numVertices,
- self.fiberDim*self.numVertices),
- dtype=numpy.float64)
- self.valsLumped = numpy.zeros( (self.fiberDim*self.numVertices,),
- dtype=numpy.float64)
-
- n = numpy.zeros( (self.fiberDim, self.fiberDim*self.numBasis),
- dtype=numpy.float64)
-
- for cell in self.cells:
- cellMatrix = numpy.zeros( (self.fiberDim*self.numBasis,
- self.fiberDim*self.numBasis),
- dtype=numpy.float64)
-
- vertices = self.vertices[cell, :]
- (jacobian, jacobianInv, jacobianDet) = \
- feutils.calculateJacobianQuad(self.quadrature, vertices)
- density = 1.0
- for iQuad in xrange(self.numQuadPts):
- n *= 0.0
- for iBasis in xrange(self.numBasis):
- for iDim in xrange(self.fiberDim):
- n[iDim, iBasis*self.fiberDim+iDim] = self.basis[iQuad, iBasis]
-
- wt = density * self.quadWts[iQuad] * jacobianDet[iQuad]
- cellMatrix[:] += wt * numpy.dot(n.transpose(), n)
- feutils.assembleMat(self.valsMatrix, cellMatrix, cell, self.fiberDim)
-
- cellVec = self._calculateMatrixLumped(cellMatrix)
- feutils.assembleVec(self.valsLumped, cellVec, cell, self.fiberDim)
- return
-
-
- def _calculateMatrixLumped(self, matrix):
- """
- Calculate lumped matrix associated with integration.
- """
-
- (nrows, ncols) = matrix.shape
- lumped = numpy.zeros( (ncols), dtype=numpy.float64)
-
- for iVertex in xrange(self.numBasis):
- i = numpy.asarray(range(self.numBasis))
- for iDim in xrange(self.fiberDim):
- iR = iVertex * self.fiberDim + iDim
- indices = self.fiberDim * i + iDim
- lumped[iR] = sum(matrix[iR, indices])
- return lumped
-
-
-# End of file
More information about the cig-commits
mailing list