[cig-commits] r7701 - short/3D/PyLith/trunk/libsrc/bc
willic3 at geodynamics.org
willic3 at geodynamics.org
Wed Jul 18 14:13:21 PDT 2007
Author: willic3
Date: 2007-07-18 14:13:20 -0700 (Wed, 18 Jul 2007)
New Revision: 7701
Modified:
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
Log:
Put in orientation stuff stolen from faults.
Maybe we should have this in a separate package.
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-07-18 21:01:09 UTC (rev 7700)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.cc 2007-07-18 21:13:20 UTC (rev 7701)
@@ -30,5 +30,109 @@
_db = 0;
} // destructor
+// ----------------------------------------------------------------------
+// Compute weighted orientation of cell face for 1-D elements.
+void
+pylith::bc::BoundaryCondition::_orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient1D
+ assert(0 != orientation);
+ assert(1 == orientation->size());
+ (*orientation) = 1.0;
+} // _orient1D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of cell face for 2-D elements.
+void
+pylith::bc::BoundaryCondition::_orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient2D
+ const int orientSize = 4;
+ assert(0 != orientation);
+ assert(orientSize == orientation->size());
+ const int jacobianSize = 2;
+ assert(jacobianSize == jacobian.size());
+ // cellDim is 1
+ const int spaceDim = 2;
+
+ const double j1 = jacobian[0];
+ const double j2 = jacobian[1];
+ (*orientation)[0] = j1;
+ (*orientation)[1] = j2;
+ (*orientation)[2] = j2;
+ (*orientation)[3] = -j1;
+} // _orient2D
+
+// ----------------------------------------------------------------------
+// Compute weighted orientation of cell face for 3-D elements.
+void
+pylith::bc::BoundaryCondition::_orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir)
+{ // _orient3D
+ const int orientSize = 9;
+ assert(0 != orientation);
+ assert(orientSize == orientation->size());
+ const int jacobianSize = 6;
+ assert(jacobianSize == jacobian.size());
+ assert(3 == upDir.size());
+
+ const int cellDim = 2;
+ const int spaceDim = 3;
+
+ const double j00 = jacobian[0];
+ const double j01 = jacobian[1];
+ const double j10 = jacobian[2];
+ const double j11 = jacobian[3];
+ const double j20 = jacobian[4];
+ const double j21 = jacobian[5];
+
+ // Compute normal using Jacobian
+ double r0 = j10*j21 - j20*j11;
+ double r1 = -j00*j21 + j20*j01;
+ double r2 = j00*j11 - j10*j01;
+ // Make unit vector
+ double mag = sqrt(r0*r0 + r1*r1 + r2*r2);
+ assert(mag > 0.0);
+ r0 /= mag;
+ r1 /= mag;
+ r2 /= mag;
+
+ // Compute along-strike direction; cross product of "up" and normal
+ double p0 = upDir[1]*r2 - upDir[2]*r1;
+ double p1 = -upDir[0]*r2 + upDir[2]*r0;
+ double p2 = upDir[0]*r1 - upDir[1]*r0;
+ // Make unit vector
+ mag = sqrt(p0*p0 + p1*p1 + p2*p2);
+ assert(mag > 0.0);
+ p0 /= mag;
+ p1 /= mag;
+ p2 /= mag;
+
+ // Compute up-dip direction; cross product of normal and along-strike
+ const double q0 = r1*p2 - r2*p1;
+ const double q1 = -r0*p2 + r2*p0;
+ const double q2 = r0*p1 - r1*p0;
+ mag = sqrt(q0*q0 + q1*q1 + q2*q2);
+ assert(mag > 0.0);
+
+ const double wt = jacobianDet;
+ (*orientation)[0] = p0*wt;
+ (*orientation)[1] = p1*wt;
+ (*orientation)[2] = p2*wt;
+ (*orientation)[3] = q0*wt;
+ (*orientation)[4] = q1*wt;
+ (*orientation)[5] = q2*wt;
+ (*orientation)[6] = r0*wt;
+ (*orientation)[7] = r1*wt;
+ (*orientation)[8] = r2*wt;
+} // _orient3D
+
+
// End of file
Modified: short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-07-18 21:01:09 UTC (rev 7700)
+++ short/3D/PyLith/trunk/libsrc/bc/BoundaryCondition.hh 2007-07-18 21:13:20 UTC (rev 7701)
@@ -22,6 +22,7 @@
#include "pylith/utils/sievetypes.hh" // USES PETSc Mesh, real_section_type
#include "pylith/utils/petscfwd.h" // USES PETScMat
+#include "pylith/utils/array.hh" // USES double_array
#include <string> // HASA std::string
@@ -98,6 +99,84 @@
void initialize(const ALE::Obj<ALE::Mesh>& mesh,
const spatialdata::geocoords::CoordSys* cs) = 0;
+ // PROTECTED TYPEDEFS /////////////////////////////////////////////////
+protected :
+
+ /// Function type for orientation methods.
+ typedef void (*orient_fn_type)(double_array*,
+ const double_array&,
+ const double,
+ const double_array&);
+
+ // PROTECTED METHODS //////////////////////////////////////////////////
+protected :
+
+ /** Compute weighted orientation of cell face for
+ * 1-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to horizontal tangent that
+ * is not collinear with face normal (usually "up" direction).
+ */
+ static
+ void _orient1D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of cell face for
+ * 2-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to horizontal tangent that is
+ * not collinear with face normal (usually "up" direction).
+ */
+ static
+ void _orient2D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
+ /** Compute weighted orientation of cell face for
+ * 3-D elements. Orientation is either at vertices or quadrature
+ * points, depending on whether the arguments have been evaluated at
+ * the vertices or quadrature points.
+ *
+ * The orientation is returned as an array of direction cosines.
+ *
+ * size = spaceDim*spaceDim
+ * index = iDir*spaceDim + iComponent
+ *
+ * @param orientation Array of direction cosines.
+ * @param jacobian Jacobian matrix at point.
+ * @param jacobianDet Determinant of Jacobian matrix at point.
+ * @param upDir Direction perpendicular to horizontal tangent that is
+ * not collinear with face normal (usually "up" direction).
+ */
+ static
+ void _orient3D(double_array* orientation,
+ const double_array& jacobian,
+ const double jacobianDet,
+ const double_array& upDir);
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
More information about the cig-commits
mailing list