[cig-commits] r4558 - in short/3D/PyLith/trunk/libsrc: . feassemble
baagaard at geodynamics.org
baagaard at geodynamics.org
Fri Sep 15 15:50:47 PDT 2006
Author: baagaard
Date: 2006-09-15 15:50:47 -0700 (Fri, 15 Sep 2006)
New Revision: 4558
Added:
short/3D/PyLith/trunk/libsrc/feassemble/
short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
Log:
Started work on quadrature interface.
Added: short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2006-09-15 18:25:07 UTC (rev 4557)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Makefile.am 2006-09-15 22:50:47 UTC (rev 4558)
@@ -0,0 +1,36 @@
+# -*- Makefile -*-
+#
+# ----------------------------------------------------------------------
+#
+# Brad T. Aagaard
+# U.S. Geological Survey
+#
+# <LicenseText>
+#
+# ----------------------------------------------------------------------
+#
+
+subpackage = feassemble
+include $(top_srcdir)/subpackage.am
+
+lib_LTLIBRARIES = libpylithfeassemble.la
+
+libpylithmeshio_la_SOURCES = \
+ Quadrature.cc
+
+subpkginclude_HEADERS = \
+ Quadrature.hh
+
+noinst_HEADERS =
+
+INCLUDES += $(PETSC_INCLUDE)
+
+# export
+clean-local: clean-subpkgincludeHEADERS
+BUILT_SOURCES = export-subpkgincludeHEADERS
+CLEANFILES = export-subpkgincludeHEADERS
+
+# version
+# $Id$
+
+# End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2006-09-15 18:25:07 UTC (rev 4557)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.cc 2006-09-15 22:50:47 UTC (rev 4558)
@@ -0,0 +1,74 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include "Quadrature.hh" // implementation of class methods
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::Quadrature::Quadrature(void) :
+ _pBasisFns(0),
+ _pBasisFnsDeriv(0),
+ _pQuadPts(0),
+ _pQuadWts(0),
+ _numDims(0),
+ _numCorners(0)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::Quadrature::~Quadrature(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+// Set basis functions and their derivatives and coordinates and
+// weights of the quadrature points.
+void
+pylith::feassemble::Quadrature::initialize(const double* pBasisFns,
+ const double* pBasisFnsDeriv,
+ const double* pQuadPts,
+ const double* pQuadWts,
+ const int numDims,
+ const int numCorners,
+ const int numQuadPts)
+{ // initialize
+ if (0 == pBasisFns ||
+ 0 == pBasisFnsDeriv ||
+ 0 == pQuadPts ||
+ 0 == pQuadWts ||
+ numDims < 1 || numDims > 3)
+ throw std::runtime_error("Basis functions and their derivatives "
+ "and quadrature points must all be specified.");
+
+ int size = numCorners * numQuadPts;
+ delete[] _pBasisFns;
+ _pBasisFns = (size > 0) ? new double[size] : 0;
+
+ size = numCorners * numQuadPts * numDims;
+ delete[] _pBasisFnsDeriv;
+ _pBasisFnsDeriv = (size > 0) ? new double[size] : 0;
+
+ size = numQuadPts * numDims;
+ delete[] _pQuadPts;
+ _pQuadPts = (size > 0) ? new double[size] : 0;
+
+ size = numQuadPts;
+ delete[] _pQuadWts;
+ _pQuadWts = (size > 0) ? new double[size] : 0;
+
+ _numDims = numDims;
+ _numCorners = numCorners;
+ _numQuadPts = numQuadPts;
+} // initialize
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-15 18:25:07 UTC (rev 4557)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature.hh 2006-09-15 22:50:47 UTC (rev 4558)
@@ -0,0 +1,128 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+// @file pylith/feassemble/Quadrature.hh
+
+// @brief Abstract base class for integrating over finite-elements
+// using quadratures.
+//
+// This object holds the basis functions and their derivatives
+// evaluated at the quadrature points, and the coordinates and weights
+// of the quadrature points.
+
+#if !defined(pylith_feassemble_quadrature_hh)
+#define pylith_feassemble_quadrature_hh
+
+namespace pylith {
+ namespace feassemble {
+ class Quadrature;
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::Quadrature
+{ // Quadrature
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ Quadrature(void);
+
+ /// Destructor
+ virtual ~Quadrature(void);
+
+ /** Compute geometric quantities for a cell.
+ *
+ */
+ virtual void compute(const Obj<section_type>& coordinates,
+ const point_type& cell,
+ double* pV,
+ double* pJacobian,
+ double* pJacobianInv,
+ double* pJacobianDet) = 0;
+
+ /** Set basis functions and their derivatives and coordinates and
+ * weights of the quadrature points.
+ *
+ * @param pBasisFns Array of basis functions evaluated at quadrature pts
+ * index = iVertex*numDims + iDimension
+ *
+ * @param pBasisFnsDeriv Array of basis function derivaties evaluated
+ * at quadrature pts
+ * index = iVertex*numDims + iDimension??
+ *
+ * @param pQuadPts Array of coordinates of quadrature points in
+ * reference element
+ * index = iQuadPt*numDims + iDimension
+ *
+ * @param pQuadWts Array of weights of quadrature points
+ * index = iQuadPt
+ *
+ * @param numDims Number of dimensions
+ * @param numCorners Number of vertices in a cell
+ * @param numQuadPts Number of quadrature points
+ */
+ void initialize(const double* pBasisFns,
+ const double* pBasisFnsDeriv,
+ const double* pQuadPts,
+ const double* pQuadWts,
+ const int numDims,
+ const int numCorners,
+ const int numQuadPts);
+
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+ double _jacobianTol; ///< Tolernace for small Jacobian determinant
+
+ /** Array of basis functions evaluated at the quadrature points.
+ *
+ * N1Qp1, N1Qp2, ...
+ * N2Qp1, N2Qp2, ...
+ *
+ * size = numCorners * numQuadPts
+ * index = iBasis*numQuadPts + iQuadPt
+ */
+ double* _pBasisFns; ///< Array of basis fns evaluated at quad pts
+
+ /** Array of basis functions evaluated at the quadrature points.
+ *
+ * N1xQp1, N1yQp1, N1zQp1, N1xQp2, N1yQp2, N1zQp2, ...
+ * N2xQp1, N2yQp1, N2zQp1, N2xQp2, N2yQp2, N2zQp2, ...
+ *
+ * size = numCorners * numQuadPts * numDeriv
+ * index = iBasis*numQuadPts*numDeriv + iQuadPt*numDeriv + iDeriv
+ */
+ double* _pBasisFnsDeriv;
+
+ /** Array of coordinates of quadrature points.
+ *
+ * size = numQuadPts * numDims
+ * index = iQuadPts*numDims + iDim
+ */
+ double* _pQuadPts; ///< Array of coordinates of quadrature points
+
+ /** Array of weights of quadrature points.
+ *
+ * WtQp1, WtQp2, ...
+ */
+ double* _pQuadWts; ///< Array of weights of quadrature points
+
+ int _numDims; ///< Number of dimensions
+ int _numCorners; ///< Number of vertices in cell
+ int _numQuadPts; ///< Number of quadrature points
+
+}; // Quadrature
+
+#endif // pylith_feassemble_quadrature_hh
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-15 18:25:07 UTC (rev 4557)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.cc 2006-09-15 22:50:47 UTC (rev 4558)
@@ -0,0 +1,97 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#include "ElemGeometry1D.hh" // implementation of class methods
+
+// ----------------------------------------------------------------------
+// Constructor
+pylith::feassemble::ElemGeometry1D::ElemGeometry1D(void)
+{ // constructor
+} // constructor
+
+// ----------------------------------------------------------------------
+// Destructor
+pylith::feassemble::ElemGeometry1D::~ElemGeometry1D(void)
+{ // destructor
+} // destructor
+
+// ----------------------------------------------------------------------
+/** Compute geometry of element object.
+ *
+ */
+void compute(const Obj<section_type>& coordinates,
+ const point_type& cell,
+ value_type* pV,
+ value_type* pJacobian,
+ valye_type* pJacobianInv,
+ value_type& jacobianDet)
+{ // compute
+ assert(0 != pV);
+ assert(0 != pJacobian);
+ assert(0 != pJacobianInv);
+
+ const int ndims = 1;
+ const topology_type::patch_type patch = 0;
+ const value_type* coords = coordinates->restrict(patch, cell);
+
+ // Matt: What is pV???
+ // Are they the coordinates of the vertices of the cell?
+ pV[0] = coords[0];
+
+
+
+
+ value_type* invDet;
+
+
+
+
+
+
+ for(int d = 0; d < _dim; d++) {
+ for(int f = 0; f < _dim; f++) {
+ J[d*_dim+f] = 0.5*(coords[(f+1)*_dim+d] - coords[0*_dim+d]);
+ }
+ }
+ if (_dim == 1) {
+ detJ = J[0];
+ } else if (_dim == 2) {
+ detJ = J[0]*J[3] - J[1]*J[2];
+ } else if (_dim == 3) {
+ detJ = J[0*3+0]*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]) +
+ J[0*3+1]*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]) +
+ J[0*3+2]*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
+ }
+ if (invJ) {
+ invDet = 1.0/detJ;
+ if (_dim == 2) {
+ invJ[0] = invDet*J[3];
+ invJ[1] = -invDet*J[1];
+ invJ[2] = -invDet*J[2];
+ invJ[3] = invDet*J[0];
+ } else if (_dim == 3) {
+ // FIX: This may be wrong
+ invJ[0*3+0] = invDet*(J[1*3+1]*J[2*3+2] - J[1*3+2]*J[2*3+1]);
+ invJ[0*3+1] = invDet*(J[1*3+2]*J[2*3+0] - J[1*3+0]*J[2*3+2]);
+ invJ[0*3+2] = invDet*(J[1*3+0]*J[2*3+1] - J[1*3+1]*J[2*3+0]);
+ invJ[1*3+0] = invDet*(J[0*3+1]*J[2*3+2] - J[0*3+2]*J[2*3+1]);
+ invJ[1*3+1] = invDet*(J[0*3+2]*J[2*3+0] - J[0*3+0]*J[2*3+2]);
+ invJ[1*3+2] = invDet*(J[0*3+0]*J[2*3+1] - J[0*3+1]*J[2*3+0]);
+ invJ[2*3+0] = invDet*(J[0*3+1]*J[1*3+2] - J[0*3+2]*J[1*3+1]);
+ invJ[2*3+1] = invDet*(J[0*3+2]*J[1*3+0] - J[0*3+0]*J[1*3+2]);
+ invJ[2*3+2] = invDet*(J[0*3+0]*J[1*3+1] - J[0*3+1]*J[1*3+0]);
+ }
+ }
+
+} // compute
+
+// End of file
Added: short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh 2006-09-15 18:25:07 UTC (rev 4557)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Quadrature1D.hh 2006-09-15 22:50:47 UTC (rev 4558)
@@ -0,0 +1,50 @@
+// -*- C++ -*-
+//
+// ======================================================================
+//
+// Brad T. Aagaard
+// U.S. Geological Survey
+//
+// {LicenseText}
+//
+// ======================================================================
+//
+
+#if !defined(pylith_feassemble_elemgeometry1d_hh)
+#define pylith_feassemble_elemgeometry1d_hh
+
+#include "ElemGeometry.hh"
+
+namespace pylith {
+ namespace feassemble {
+ class ElemGeometry1D;
+ } // feassemble
+} // pylith
+
+class pylith::feassemble::ElemGeometry1D
+{ // ElemGeometry1D
+
+// PUBLIC MEMBERS ///////////////////////////////////////////////////////
+public :
+
+ /// Constructor
+ ElemGeometry1D(void);
+
+ /// Destructor
+ ~ElemGeometry1D(void);
+
+ /** Compute geometry of element object.
+ *
+ */
+ void compute(const Obj<section_type>& coordinates,
+ const point_type& cell,
+ value_type* pV,
+ value_type* pJacobian,
+ valye_type* pJacobianInv,
+ value_type* pJacobianDet);
+
+}; // ElemGeometry1D
+
+#endif // pylith_feassemble_elemgeometry1d_hh
+
+// End of file
More information about the cig-commits
mailing list