[cig-commits] r19033 - in short/3D/PyLith/branches/v1.6-revisedfault: libsrc/pylith/faults libsrc/pylith/problems pylith pylith/faults pylith/feassemble unittests/pytests/feassemble
brad at geodynamics.org
brad at geodynamics.org
Thu Oct 6 11:06:56 PDT 2011
Author: brad
Date: 2011-10-06 11:06:55 -0700 (Thu, 06 Oct 2011)
New Revision: 19033
Added:
short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATQuadrature.py
Modified:
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/problems/Solver.cc
short/3D/PyLith/branches/v1.6-revisedfault/pylith/Makefile.am
short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesive.py
short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesiveKin.py
short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATLagrange.py
short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATSimplex.py
short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/__init__.py
short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATLagrange.py
short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATSimplex.py
Log:
Merge from stable. Initial cleanup of fault implementation using collocated vertices/quadrature. Switch back to diagonal custom fault preconditioner.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-10-06 18:06:55 UTC (rev 19033)
@@ -787,7 +787,7 @@
*/
const double assemblyFactor = 5.0;
- //#define USE_DIAGONAL_PRECONDITIONER // Extract only the diagonal terms of P.
+#define USE_DIAGONAL_PRECONDITIONER // Extract only the diagonal terms of P.
const int setupEvent = _logger->eventId("FaPr setup");
const int computeEvent = _logger->eventId("FaPr compute");
@@ -813,12 +813,13 @@
int_array indicesLagrange(nrowsF);
#if defined(USE_DIAGONAL_PRECONDITIONER)
double_array preconditionerCell(nrowsF);
+ double_array basisProducts(numBasis);
#else
double_array preconditionerCell(matrixSizeF);
+ double_array basisProducts(numBasis*numBasis);
#endif
double_array jacobianCellP(matrixSizeF);
double_array jacobianCellN(matrixSizeF);
- double_array basisProducts(numBasis*numBasis);
// Get sections
const ALE::Obj<RealSection>& solutionSection = fields->solution().section();
@@ -965,6 +966,7 @@
const double_array& basis = _quadrature->basis();
const double_array& jacobianDet = _quadrature->jacobianDet();
+#if defined(USE_DIAGONAL_PRECONDITIONER)
// Compute product of basis functions.
// Want values summed over quadrature points
basisProducts = 0.0;
@@ -972,6 +974,17 @@
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
+ basisProducts[iBasis] += wt*basis[iQ+iBasis]*basis[iQ+iBasis];
+ } // for
+ } // for
+#else
+ // Compute product of basis functions.
+ // Want values summed over quadrature points
+ basisProducts = 0.0;
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+
+ for (int iBasis=0, iQ=iQuad*numBasis; iBasis < numBasis; ++iBasis) {
const double valI = wt*basis[iQ+iBasis];
for (int jBasis=0; jBasis < numBasis; ++jBasis) {
@@ -980,6 +993,7 @@
} // for
} // for
} // for
+#endif
preconditionerCell = 0.0;
@@ -988,24 +1002,12 @@
// First index Lagrange vertices
const int iB = iBasis*spaceDim;
- const int jBasis = iBasis; // Compute only diagonal entries
-
- // First index for Lagrange vertices
- const int jB = jBasis*spaceDim;
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ const int k = (iB+iDim)*nrowsF+(iB+iDim);
+ const double jacobianNPInv =
+ 1.0/jacobianCellN[k] + 1.0/jacobianCellP[k];
- for (int kBasis=0; kBasis < numBasis; ++kBasis) {
- const int kB = kBasis*spaceDim;
-
- const double valIK = basisProducts[iBasis*numBasis+kBasis];
- const double valJK = basisProducts[jBasis*numBasis+kBasis];
-
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- const int k = (kB+iDim)*nrowsF+(kB+iDim);
- const double jacobianNPInv =
- 1.0/jacobianCellN[k] + 1.0/jacobianCellP[k];
-
- preconditionerCell[iB+iDim] -= valIK*valJK*jacobianNPInv;
- } // for
+ preconditionerCell[iB+iDim] -= jacobianNPInv*basisProducts[iBasis];
} // for
} // for
#else
@@ -1039,8 +1041,8 @@
} // for
} // for
} // for
+ preconditionerCell *= assemblyFactor;
#endif
- preconditionerCell *= assemblyFactor;
#if 0 // DEBUGGING
#if defined(USE_DIAGONAL_PRECONDITIONER)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/problems/Solver.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/problems/Solver.cc 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/problems/Solver.cc 2011-10-06 18:06:55 UTC (rev 19033)
@@ -189,7 +189,7 @@
err = MatSetType(_jacobianPCFault, MATAIJ);
err = MatSetFromOptions(_jacobianPCFault); CHECK_PETSC_ERROR(err);
-#if 0
+#if 1
// Allocate just the diagonal.
err = MatSeqAIJSetPreallocation(_jacobianPCFault, 1,
PETSC_NULL); CHECK_PETSC_ERROR(err);
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/Makefile.am
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/Makefile.am 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/Makefile.am 2011-10-06 18:06:55 UTC (rev 19033)
@@ -54,6 +54,7 @@
feassemble/ElasticityExplicitLgDeform.py \
feassemble/ElasticityImplicit.py \
feassemble/ElasticityImplicitLgDeform.py \
+ feassemble/FIATQuadrature.py \
feassemble/FIATLagrange.py \
feassemble/FIATSimplex.py \
feassemble/Integrator.py \
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesive.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesive.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesive.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -92,6 +92,10 @@
# TEMPORARY
ModuleFaultCohesive.faultMeshFilename(self,
self.inventory.meshFilename)
+
+ # Hardwire collocated quadrature
+ self.faultQuadrature.cell.collocateQuad = True
+ self.faultQuadrature.cell.order = self.faultQuadrature.cell.degree
return
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesiveKin.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesiveKin.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/faults/FaultCohesiveKin.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -73,7 +73,6 @@
output = pyre.inventory.facility("output", family="output_manager",
factory=OutputFaultKin)
output.meta['tip'] = "Output manager associated with fault data."
-
# PUBLIC METHODS /////////////////////////////////////////////////////
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATLagrange.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATLagrange.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -35,6 +35,7 @@
raise ValueError("Dimension of Lagrange element must be 1, 2, or 3.")
return dim
+
# FIATLagrange class
class FIATLagrange(ReferenceCell):
"""
@@ -53,9 +54,10 @@
## Python object for managing FIATLagrange facilities and properties.
##
## \b Properties
- ## @li \b dimension Dimension of finite-element cell
- ## @li \b degree Degree of finite-element cell
- ## @li \b quad_order Order of quadrature rule
+ ## @li \b dimension Dimension of finite-element cell.
+ ## @li \b degree Degree of finite-element cell.
+ ## @li \b quad_order Order of quadrature rule.
+ ## @li \b collocate_quad Collocate quadrature points with vertices.
##
## \b Facilities
## @li None
@@ -72,7 +74,10 @@
order = pyre.inventory.int("quad_order", default=-1)
order.meta['tip'] = "Order of quadrature rule."
+ collocateQuad = pyre.inventory.bool("collocate_quad", default=False)
+ collocateQuad.meta['tip'] = "Collocate quadrature points with vertices."
+
# PUBLIC METHODS /////////////////////////////////////////////////////
def __init__(self, name="fiatlagrange"):
@@ -411,6 +416,7 @@
self.cellDim = self.inventory.dimension
self.degree = self.inventory.degree
self.order = self.inventory.order
+ self.collocateQuad = self.inventory.collocateQuad
if self.order == -1:
self.order = self.degree+1
@@ -460,11 +466,18 @@
"""
Setup quadrature rule for reference cell.
"""
+ from FIAT.reference_element import default_simplex
from FIAT.quadrature import make_quadrature
- from FIAT.reference_element import default_simplex
- return make_quadrature(default_simplex(1), self.order)
+ from FIATQuadrature import CollocatedQuadratureRule
+
+ if not self.collocateQuad:
+ q = make_quadrature(default_simplex(1), self.order)
+ else:
+ q = CollocatedQuadratureRule(default_simplex(1), self.order)
+ return q
+
def _setupElement(self):
"""
Setup the finite element for reference cell.
Copied: short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATQuadrature.py (from rev 19030, short/3D/PyLith/branches/v1.6-stable/pylith/feassemble/FIATQuadrature.py)
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATQuadrature.py (rev 0)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATQuadrature.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -0,0 +1,39 @@
+#!/usr/bin/env python
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard, U.S. Geological Survey
+# Charles A. Williams, GNS Science
+# Matthew G. Knepley, University of Chicago
+#
+# This code was developed as part of the Computational Infrastructure
+# for Geodynamics (http://geodynamics.org).
+#
+# Copyright (c) 2010-2011 University of California, Davis
+#
+# See COPYING for license information.
+#
+# ----------------------------------------------------------------------
+#
+
+## @file pylith/feassemble/FIATQuadrature.py
+##
+## @brief Python object for special FIAT quadrature schemes.
+
+from FIAT.quadrature import QuadratureRule
+class CollocatedQuadratureRule(QuadratureRule):
+ """
+ Quadrature points colocated with vertices.
+ """
+ def __init__(self, ref_el, m):
+ from FIAT.lagrange import Lagrange
+ vertices = Lagrange(ref_el, m).dual.get_nodes()
+ pts = [v.get_point_dict().keys()[0] for v in vertices]
+ npts = len(pts)
+ wts = (ref_el.volume()/npts,)*npts
+
+ QuadratureRule.__init__(self, ref_el, pts, wts)
+ return
+
+
+# End of file
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATSimplex.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATSimplex.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/FIATSimplex.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -39,6 +39,7 @@
name)
return name
+
# FIATSimplex class
class FIATSimplex(ReferenceCell):
"""
@@ -76,6 +77,9 @@
order = pyre.inventory.int("quad_order", default=-1)
order.meta['tip'] = "Order of quadrature rule [-1, order = degree]."
+ collocateQuad = pyre.inventory.bool("collocate_quad", default=False)
+ collocateQuad.meta['tip'] = "Collocate quadrature points with vertices."
+
# PUBLIC METHODS /////////////////////////////////////////////////////
@@ -163,6 +167,7 @@
self.shape = self.inventory.shape
self.degree = self.inventory.degree
self.order = self.inventory.order
+ self.collocateQuad = self.inventory.collocateQuad
if self.order == -1:
self.order = self.degree
except ValueError, err:
@@ -254,11 +259,18 @@
"""
Setup quadrature rule for reference cell.
"""
-
- import FIAT.quadrature
- return FIAT.quadrature.make_quadrature(self._getShape(), self.order)
+ from FIAT.reference_element import default_simplex
+ from FIAT.quadrature import make_quadrature
+ from FIATQuadrature import CollocatedQuadratureRule
+
+ if not self.collocateQuad:
+ q = make_quadrature(self._getShape(), self.order)
+ else:
+ q = CollocatedQuadratureRule(self._getShape(), self.order)
+ return q
+
def _setupBasisFns(self):
"""
Setup basis functions for reference cell.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/__init__.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/__init__.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/pylith/feassemble/__init__.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -26,6 +26,7 @@
'ElasticityExplicitLgDeform',
'ElasticityImplicit',
'ElasticityImplicitLgDeform',
+ 'FIATQuadrature',
'FIATLagrange',
'FIATSimplex',
'IntegratorElasticity',
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATLagrange.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATLagrange.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATLagrange.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -74,6 +74,40 @@
# ----------------------------------------------------------------------
+class Line2Collocated(Line2):
+
+ def __init__(self):
+ """
+ Setup line2 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0]])
+ quadPts = vertices[:]
+ quadWts = numpy.array( [1.0,1.0], dtype=numpy.float64 )
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (2, 2), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (2, 2, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q)],
+ dtype=numpy.float64).reshape( (2,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)]],
+ dtype=numpy.float64)
+ basisDeriv[iQuad] = deriv.reshape((2, 1))
+ iQuad += 1
+
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class Line3(object):
def __init__(self):
@@ -207,6 +241,49 @@
# ----------------------------------------------------------------------
+class Quad4Collocated(Quad4):
+
+ def __init__(self):
+ """
+ Setup quad4 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0],
+ [+1.0, -1.0],
+ [+1.0, +1.0],
+ [-1.0, +1.0]])
+ quadPts = numpy.array([[-1.0, -1.0],
+ [+1.0, -1.0],
+ [-1.0, +1.0],
+ [+1.0, +1.0]])
+ quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (4, 4), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (4, 4, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q)],
+ dtype=numpy.float64).reshape( (4,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)],
+ [self.N3p(q), self.N3q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((4, 2))
+ iQuad += 1
+
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class Quad9(object):
def __init__(self):
@@ -505,6 +582,63 @@
# ----------------------------------------------------------------------
+class Hex8Collocated(Hex8):
+
+ def __init__(self):
+ """
+ Setup hex8 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [+1.0, +1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0],
+ [+1.0, -1.0, +1.0],
+ [+1.0, +1.0, +1.0],
+ [-1.0, +1.0, +1.0]])
+ quadPts = numpy.array([[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [+1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0],
+ [+1.0, -1.0, +1.0],
+ [-1.0, +1.0, +1.0],
+ [+1.0, +1.0, +1.0]])
+ quadWts = numpy.array( [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (8, 8), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (8, 8, 3), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q),
+ self.N4(q), self.N5(q),
+ self.N6(q), self.N7(q)],
+ dtype=numpy.float64).reshape( (8,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q), self.N0r(q)],
+ [self.N1p(q), self.N1q(q), self.N1r(q)],
+ [self.N2p(q), self.N2q(q), self.N2r(q)],
+ [self.N3p(q), self.N3q(q), self.N3r(q)],
+ [self.N4p(q), self.N4q(q), self.N4r(q)],
+ [self.N5p(q), self.N5q(q), self.N5r(q)],
+ [self.N6p(q), self.N6q(q), self.N6r(q)],
+ [self.N7p(q), self.N7q(q), self.N7r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((8, 3))
+ iQuad += 1
+
+ self.cellDim = 3
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class Hex27(object):
def __init__(self):
@@ -1014,6 +1148,26 @@
return
+ def test_initialize_line2_collocated(self):
+ """
+ Test initialize() with line2 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 1
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=1)
+
+
+ cellE = Line2Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryLine1D
+ self.failUnless(isinstance(cell.geometry, GeometryLine1D))
+ return
+
+
def test_initialize_line3(self):
"""
Test initialize() with line3 cell.
@@ -1050,6 +1204,25 @@
return
+ def test_initialize_quad4_collocated(self):
+ """
+ Test initialize() with quad4 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 2
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Quad4Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryQuad2D
+ self.failUnless(isinstance(cell.geometry, GeometryQuad2D))
+ return
+
+
def test_initialize_quad9(self):
"""
Test initialize() with quad9 cell.
@@ -1086,6 +1259,25 @@
return
+ def test_initialize_hex8_collocated(self):
+ """
+ Test initialize() with hex8 cell.
+ """
+ cell = FIATLagrange()
+ cell.inventory.dimension = 3
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=3)
+
+ cellE = Hex8Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryHex3D
+ self.failUnless(isinstance(cell.geometry, GeometryHex3D))
+ return
+
+
def test_initialize_hex27(self):
"""
Test initialize() with hex27 cell.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATSimplex.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATSimplex.py 2011-10-06 18:05:23 UTC (rev 19032)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/pytests/feassemble/TestFIATSimplex.py 2011-10-06 18:06:55 UTC (rev 19033)
@@ -73,6 +73,40 @@
return 0.5
# ----------------------------------------------------------------------
+class Line2Collocated(Line2):
+
+ def __init__(self):
+ """
+ Setup line2 cell.
+ """
+ vertices = numpy.array([[-1.0], [1.0]])
+ quadPts = vertices[:]
+ quadWts = numpy.array( [1.0, 1.0], dtype=numpy.float64 )
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (2, 2), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (2, 2, 1), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q)],
+ dtype=numpy.float64).reshape( (2,) )
+ deriv = numpy.array([[self.N0p(q)], [self.N1p(q)]],
+ dtype=numpy.float64)
+ basisDeriv[iQuad] = deriv.reshape((2, 1))
+ iQuad += 1
+
+ self.cellDim = 1
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class Line3(object):
def __init__(self):
@@ -197,6 +231,43 @@
# ----------------------------------------------------------------------
+class Tri3Collocated(Tri3):
+
+ def __init__(self):
+ """
+ Setup tri33 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0],
+ [+1.0, -1.0],
+ [-1.0, +1.0]])
+ quadPts = vertices[:]
+ quadWts = numpy.array( [2.0/3.0, 2.0/3.0, 2.0/3.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (3, 3), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (3, 3, 2), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q), self.N2(q)],
+ dtype=numpy.float64).reshape( (3,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q)],
+ [self.N1p(q), self.N1q(q)],
+ [self.N2p(q), self.N2q(q)]])
+ basisDeriv[iQuad] = deriv.reshape((3, 2))
+ iQuad += 1
+
+ self.cellDim = 2
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class Tri6(object):
def __init__(self):
@@ -393,6 +464,46 @@
# ----------------------------------------------------------------------
+class Tet4Collocated(Tet4):
+
+ def __init__(self):
+ """
+ Setup tri33 cell.
+ """
+ vertices = numpy.array([[-1.0, -1.0, -1.0],
+ [+1.0, -1.0, -1.0],
+ [-1.0, +1.0, -1.0],
+ [-1.0, -1.0, +1.0]])
+ quadPts = vertices[:]
+ quadWts = numpy.array( [1.0/3.0, 1.0/3.0, 1.0/3.0, 1.0/3.0])
+
+ # Compute basis fns and derivatives at quadrature points
+ basis = numpy.zeros( (4, 4), dtype=numpy.float64)
+ basisDeriv = numpy.zeros( (4, 4, 3), dtype=numpy.float64)
+ iQuad = 0
+ for q in quadPts:
+ basis[iQuad] = numpy.array([self.N0(q), self.N1(q),
+ self.N2(q), self.N3(q)],
+ dtype=numpy.float64).reshape( (4,) )
+ deriv = numpy.array([[self.N0p(q), self.N0q(q), self.N0r(q)],
+ [self.N1p(q), self.N1q(q), self.N1r(q)],
+ [self.N2p(q), self.N2q(q), self.N2r(q)],
+ [self.N3p(q), self.N3q(q), self.N3r(q)]])
+ basisDeriv[iQuad] = deriv.reshape((4, 3))
+ iQuad += 1
+
+ self.cellDim = 3
+ self.numCorners = len(vertices)
+ self.numQuadPts = len(quadPts)
+ self.vertices = vertices
+ self.quadPts = quadPts
+ self.quadWts = quadWts
+ self.basis = basis
+ self.basisDeriv = basisDeriv
+ return
+
+
+# ----------------------------------------------------------------------
class TestFIATSimplex(unittest.TestCase):
"""
Unit testing of FIATSimplex object.
@@ -426,9 +537,10 @@
Test initialize() with line2 cell.
"""
cell = FIATSimplex()
- cell.shape = "line"
- cell.degree = 1
- cell.order = 1
+ cell.inventory.shape = "line"
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell._configure()
cell.initialize(spaceDim=1)
cellE = Line2()
@@ -438,14 +550,34 @@
return
+ def test_initialize_line2_collodated(self):
+ """
+ Test initialize() with line2 cell.
+ """
+ cell = FIATSimplex()
+ cell.inventory.shape = "line"
+ cell.inventory.degree = 1
+ cell.inventory.order = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=1)
+
+ cellE = Line2Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryLine1D
+ self.failUnless(isinstance(cell.geometry, GeometryLine1D))
+ return
+
+
def test_initialize_line3(self):
"""
Test initialize() with line3 cell.
"""
cell = FIATSimplex()
- cell.shape = "line"
- cell.degree = 2
- cell.order = 2
+ cell.inventory.shape = "line"
+ cell.inventory.degree = 2
+ cell.inventory.order = 2
+ cell._configure()
cell.initialize(spaceDim=2)
cellE = Line3()
@@ -472,6 +604,24 @@
return
+ def test_initialize_tri3_collocated(self):
+ """
+ Test initialize() with tri3 cell.
+ """
+ cell = FIATSimplex()
+ cell.inventory.shape = "triangle"
+ cell.inventory.degree = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=2)
+
+ cellE = Tri3Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryTri2D
+ self.failUnless(isinstance(cell.geometry, GeometryTri2D))
+ return
+
+
def test_initialize_tri6(self):
"""
Test initialize() with tri6 cell.
@@ -506,6 +656,24 @@
return
+ def test_initialize_tet4_collocated(self):
+ """
+ Test initialize() with tet4 cell.
+ """
+ cell = FIATSimplex()
+ cell.inventory.shape = "tetrahedron"
+ cell.inventory.degree = 1
+ cell.inventory.collocateQuad = True
+ cell._configure()
+ cell.initialize(spaceDim=3)
+
+ cellE = Tet4Collocated()
+ self._checkVals(cellE, cell)
+ from pylith.feassemble.CellGeometry import GeometryTet3D
+ self.failUnless(isinstance(cell.geometry, GeometryTet3D))
+ return
+
+
def test_factory(self):
"""
Test factory method.
More information about the CIG-COMMITS
mailing list