[cig-commits] r13131 - in short/3D/PyLith/trunk: libsrc/feassemble libsrc/utils modulesrc/feassemble pylith/feassemble/quadrature unittests/libtests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Sat Oct 25 19:43:39 PDT 2008
Author: brad
Date: 2008-10-25 19:43:38 -0700 (Sat, 25 Oct 2008)
New Revision: 13131
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
short/3D/PyLith/trunk/libsrc/utils/Makefile.am
short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh
Log:
Added switch to turn ill-conditioning check on/off.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/libsrc/feassemble/ElasticityImplicit.cc 2008-10-26 02:43:38 UTC (rev 13131)
@@ -21,6 +21,7 @@
#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
#include "pylith/utils/array.hh" // USES double_array
#include "pylith/utils/macrodefs.h" // USES CALL_MEMBER_FN
+#include "pylith/utils/blas.h" // USES dgesvd
#include "petscmat.h" // USES PetscMat
#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
@@ -280,9 +281,6 @@
} // for
} // integrateResidual
-extern "C" {
- EXTERN void dgesvd_(const char*,const char*,PetscBLASInt *,PetscBLASInt*,PetscScalar *,PetscBLASInt*,PetscReal*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscScalar*,PetscBLASInt*,PetscBLASInt*);
-}
// ----------------------------------------------------------------------
// Compute stiffness matrix.
void
@@ -408,27 +406,35 @@
CALL_MEMBER_FN(*this, elasticityJacobianFn)(elasticConsts);
-#if 0
- int n = numBasis*spaceDim, lwork = 5*n, idummy, lierr;
- double *elemMat = new double[n*n];
- double *svalues = new double[n];
- double *work = new double[lwork];
- double minSV, maxSV, sdummy;
+ if (_quadrature->checkConditioning()) {
+ int n = numBasis*spaceDim;
+ int lwork = 5*n;
+ int idummy = 0;
+ int lierr = 0;
+ double *elemMat = new double[n*n];
+ double *svalues = new double[n];
+ double *work = new double[lwork];
+ double minSV = 0;
+ double maxSV = 0;
+ double sdummy = 0;
- for(int i = 0; i < n*n; ++i) {elemMat[i] = _cellMatrix[i];}
- dgesvd_("N", "N", &n, &n, elemMat, &n, svalues, &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr);
- if (lierr) {throw std::runtime_error("Lapack SVD failed");}
- minSV = svalues[n-1];
- maxSV = svalues[0];
- std::cout << "Element " << c_index << std::endl;
- for(int i = 0; i < n; ++i) {
- std::cout << " sV["<<i<<"] = " << svalues[i] << std::endl;
- }
- std::cout << " kappa(elemMat) = " << maxSV/minSV << std::endl;
- delete [] elemMat;
- delete [] svalues;
- delete [] work;
-#endif
+ const int n2 = n*n;
+ for (int i = 0; i < n2; ++i)
+ elemMat[i] = _cellMatrix[i];
+ dgesvd("N", "N", &n, &n, elemMat, &n, svalues,
+ &sdummy, &idummy, &sdummy, &idummy, work, &lwork, &lierr);
+ if (lierr)
+ throw std::runtime_error("Lapack SVD failed");
+ minSV = svalues[n-1];
+ maxSV = svalues[0];
+ std::cout << "Element " << c_index << std::endl;
+ for(int i = 0; i < n; ++i)
+ std::cout << " sV["<<i<<"] = " << svalues[i] << std::endl;
+ std::cout << " kappa(elemMat) = " << maxSV/minSV << std::endl;
+ delete [] elemMat;
+ delete [] svalues;
+ delete [] work;
+ } // if
// Assemble cell contribution into field. Not sure if this is correct for
// global stiffness matrix.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2008-10-26 02:43:38 UTC (rev 13131)
@@ -30,7 +30,8 @@
_numQuadPts(0),
_spaceDim(0),
_geometry(0),
- _precomputed(false)
+ _precomputed(false),
+ _checkConditioning(false)
{ // constructor
_quadPtsPre = new real_section_type(PETSC_COMM_WORLD);
_jacobianPre = new real_section_type(PETSC_COMM_WORLD);
@@ -64,7 +65,8 @@
_numQuadPts(q._numQuadPts),
_spaceDim(q._spaceDim),
_geometry(0),
- _precomputed(q._precomputed)
+ _precomputed(q._precomputed),
+ _checkConditioning(q._checkConditioning)
{ // copy constructor
if (0 != q._geometry)
_geometry = q._geometry->clone();
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2008-10-26 02:43:38 UTC (rev 13131)
@@ -127,8 +127,20 @@
*
* @returns Minimum allowable value for Jacobian
*/
- double minJacobian(void);
+ double minJacobian(void) const;
+ /** Set flag for checking ill-conditioning.
+ *
+ * @param flag True to check for ill-conditioning, false otherwise.
+ */
+ void checkConditioning(const bool flag);
+
+ /** Get flag for checking ill-conditioning.
+ *
+ * @returns True if checking for ill-conditioning, false otherwise.
+ */
+ bool checkConditioning(void) const;
+
/** Get coordinates of quadrature points in reference cell.
*
* @returns Array of coordinates of quadrature points in reference cell.
@@ -367,7 +379,6 @@
CellGeometry* _geometry; ///< Geometry of reference cell
- bool _precomputed;
/* Precomputation sections */
int _qTag, _jTag, _jDTag, _jITag, _bTag;
Obj<real_section_type> _quadPtsPre;
@@ -375,6 +386,9 @@
Obj<real_section_type> _jacobianDetPre;
Obj<real_section_type> _jacobianInvPre;
Obj<real_section_type> _basisDerivPre;
+
+ bool _precomputed; ///< True if we have computed geometry info
+ bool _checkConditioning; ///< True if checking for ill-conditioning
}; // Quadrature
#include "Quadrature.icc" // inline methods
Modified: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.icc 2008-10-26 02:43:38 UTC (rev 13131)
@@ -17,7 +17,7 @@
// Get minimum allowable Jacobian.
inline
double
-pylith::feassemble::Quadrature::minJacobian(void) {
+pylith::feassemble::Quadrature::minJacobian(void) const {
return _minJacobian;
}
@@ -28,6 +28,20 @@
_minJacobian = min;
}
+// Set flag for checking ill-conditioning.
+inline
+void
+pylith::feassemble::Quadrature::checkConditioning(const bool flag) {
+ _checkConditioning = flag;
+}
+
+// Get flag for checking ill-conditioning.
+inline
+bool
+pylith::feassemble::Quadrature::checkConditioning(void) const {
+ return _checkConditioning;
+}
+
// Get coordinates of quadrature points in reference cell.
inline
const pylith::double_array&
Modified: short/3D/PyLith/trunk/libsrc/utils/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/utils/Makefile.am 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/libsrc/utils/Makefile.am 2008-10-26 02:43:38 UTC (rev 13131)
@@ -16,6 +16,7 @@
subpkginclude_HEADERS = \
EventLogger.hh \
EventLogger.icc \
+ blas.h \
array.hh \
arrayfwd.hh \
constdefs.h \
Modified: short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src
===================================================================
--- short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/modulesrc/feassemble/feassemble.pyxe.src 2008-10-26 02:43:38 UTC (rev 13131)
@@ -651,6 +651,50 @@
return Quadrature_minJacobian_get(self.thisptr)
+ property checkConditioning:
+ def __set__(self, flag):
+ """Set minimum allowable Jacobian."""
+ # create shim for method 'checkConditioning'
+ #embed{ void Quadrature_checkConditioning_set(void* objVptr, int flag)
+ try {
+ assert(0 != objVptr);
+ ((pylith::feassemble::Quadrature*) objVptr)->checkConditioning(flag);
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ #}embed
+ Quadrature_checkConditioning_set(self.thisptr, flag)
+
+ def __get__(self):
+ """Get minimum allowable Jacobian."""
+ # create shim for method 'checkConditioning'
+ #embed{ int Quadrature_checkConditioning_get(void* objVptr)
+ int result = 0;
+ try {
+ assert(0 != objVptr);
+ result = ((pylith::feassemble::Quadrature*) objVptr)->checkConditioning();
+ } catch (const std::exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.what()));
+ } catch (const ALE::Exception& err) {
+ PyErr_SetString(PyExc_RuntimeError,
+ const_cast<char*>(err.msg().c_str()));
+ } catch (...) {
+ PyErr_SetString(PyExc_RuntimeError,
+ "Caught unknown C++ exception.");
+ } // try/catch
+ return result;
+ #}embed
+ return Quadrature_checkConditioning_get(self.thisptr)
+
+
property refGeometry:
def __set__(self, value):
"""
Modified: short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py
===================================================================
--- short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/pylith/feassemble/quadrature/Quadrature.py 2008-10-26 02:43:38 UTC (rev 13131)
@@ -39,6 +39,8 @@
##
## \b Properties
## @li \b min_jacobian Minimum allowable determinant of Jacobian.
+ ## @li \b check_conditoning Check element matrices for
+ ## ill-conditioning.
##
## \b Facilities
## @li \b cell Reference cell with basis functions and quadrature rules
@@ -48,6 +50,11 @@
minJacobian = pyre.inventory.float("min_jacobian", default=1.0e-06)
minJacobian.meta['tip'] = "Minimum allowable determinant of Jacobian."
+ checkConditioning = pyre.inventory.bool("check_conditioning",
+ default=False)
+ checkConditioning.meta['tip'] = \
+ "Check element matrices for ill-conditioning."
+
from pylith.feassemble.FIATSimplex import FIATSimplex
cell = pyre.inventory.facility("cell", family="reference_cell",
factory=FIATSimplex)
@@ -74,6 +81,7 @@
self._createCppHandle()
self.cppHandle.minJacobian = self.minJacobian
+ self.cppHandle.checkConditioning = self.checkConditioning
self._info.log("Initializing reference cell.")
cell = self.cell
@@ -109,6 +117,7 @@
"""
Component._configure(self)
self.minJacobian = self.inventory.minJacobian
+ self.checkConditioning = self.inventory.checkConditioning
self.cell = self.inventory.cell
return
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.cc 2008-10-26 02:43:38 UTC (rev 13131)
@@ -31,6 +31,7 @@
{ // testClone
// Semi-random values manually set to check cloning
const double minJacobianE = 1.0;
+ const bool checkConditioning = true;
const int cellDimE = 1;
const int numBasisE = 2;
const int numQuadPtsE = 1;
@@ -48,6 +49,7 @@
// Set values
Quadrature1D qOrig;
qOrig._minJacobian = minJacobianE;
+ qOrig._checkConditioning = checkConditioning;
qOrig._cellDim = cellDimE;
qOrig._numBasis = numBasisE;
qOrig._numQuadPts = numQuadPtsE;
@@ -94,6 +96,7 @@
CPPUNIT_ASSERT(0 != qCopy);
CPPUNIT_ASSERT_EQUAL(minJacobianE, qCopy->_minJacobian);
+ CPPUNIT_ASSERT_EQUAL(checkConditioning, qCopy->_checkConditioning);
CPPUNIT_ASSERT_EQUAL(cellDimE, qCopy->cellDim());
CPPUNIT_ASSERT_EQUAL(numBasisE, qCopy->numBasis());
CPPUNIT_ASSERT_EQUAL(numQuadPtsE, qCopy->numQuadPts());
@@ -167,6 +170,19 @@
} // testMinJacobian
// ----------------------------------------------------------------------
+// Test checkConditioning()
+void
+pylith::feassemble::TestQuadrature::testCheckConditioning(void)
+{ // testCheckConditioning
+ Quadrature1D q;
+ CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
+ q.checkConditioning(true);
+ CPPUNIT_ASSERT_EQUAL(true, q.checkConditioning());
+ q.checkConditioning(false);
+ CPPUNIT_ASSERT_EQUAL(false, q.checkConditioning());
+} // testCheckConditioning
+
+// ----------------------------------------------------------------------
// Test refGeometry()
void
pylith::feassemble::TestQuadrature::testRefGeometry(void)
Modified: short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh
===================================================================
--- short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh 2008-10-26 02:43:02 UTC (rev 13130)
+++ short/3D/PyLith/trunk/unittests/libtests/feassemble/TestQuadrature.hh 2008-10-26 02:43:38 UTC (rev 13131)
@@ -41,6 +41,7 @@
CPPUNIT_TEST( testClone );
CPPUNIT_TEST( testMinJacobian );
+ CPPUNIT_TEST( testCheckConditioning );
CPPUNIT_TEST( testRefGeometry );
CPPUNIT_TEST( testInitialize );
@@ -55,6 +56,9 @@
/// Test minJacobian()
void testMinJacobian(void);
+ void testCheckConditioning(void);
+ /// Test checkConditioning()
+
/// Test refGeometry()
void testRefGeometry(void);
More information about the CIG-COMMITS
mailing list