[cig-commits] r11169 - short/3D/PyLith/trunk/libsrc/faults

[brad at geodynamics.org]

Author: brad
Date: 2008-02-15 14:11:44 -0800 (Fri, 15 Feb 2008)
New Revision: 11169

Modified:
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
   short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
Log:
Added computation of change in tractions on fault surface.

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2008-02-15 22:09:53 UTC (rev 11168)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.cc	2008-02-15 22:11:44 UTC (rev 11169)
@@ -470,13 +470,16 @@
 { // cellField
   assert(!_faultMesh.isNull());
   assert(!_orientation.isNull());
+  assert(0 != fields);
   assert(0 != _eqsrc);
 
   const int cohesiveDim = _faultMesh->getDimension();
 
   if (0 == strcasecmp("traction_change", name)) {
     _allocateBufferCellVector();
-    // ADD STUFF HERE
+    *fieldType = VECTOR_FIELD;
+    const ALE::Obj<real_section_type>& solution = fields->getSolution();
+    _calcTractionsChange(&_bufferCellVector, solution);
     return _bufferCellVector;
   } // if
 
@@ -776,6 +779,53 @@
 } // _calcConditioning
 
 // ----------------------------------------------------------------------
+// Compute change in tractions on fault surface using solution.
+void
+pylith::faults::FaultCohesiveKin::_calcTractionsChange(
+				 ALE::Obj<real_section_type>* tractions,
+				 const ALE::Obj<real_section_type>& solution)
+{ // _calcTractionsChange
+  assert(0 != tractions);
+  assert(!tractions->isNull());
+  assert(!solution.isNull());
+  assert(!_faultMesh.isNull());
+  assert(!_pseudoStiffness.isNull());
+  assert(0 != _quadrature);
+
+  const ALE::Obj<Mesh::label_sequence>& cells = 
+    _faultMesh->heightStratum(0);
+  const Mesh::label_sequence::iterator cellsEnd = cells->end();
+  
+  const int numBasis = _quadrature->numBasis();
+  const int spaceDim = _quadrature->spaceDim();
+  const int numQuadPts = _quadrature->numQuadPts();
+  const double_array& basis = _quadrature->basis();
+
+  double_array solutionCell(numBasis*spaceDim);
+  double_array stiffnessCell(numBasis);
+  double_array tractionsCell(numQuadPts*spaceDim);
+
+  int count = 0;
+  for (Mesh::label_sequence::iterator c_iter=cells->begin();
+       c_iter != cellsEnd;
+       ++c_iter, ++count) {
+    _faultMesh->restrict(solution, *c_iter, 
+			 &solutionCell[0], solutionCell.size());
+    _faultMesh->restrict(_pseudoStiffness, *c_iter, 
+			 &stiffnessCell[0], stiffnessCell.size());
+    for (int iQuad=0; iQuad < numQuadPts; ++iQuad)
+      for (int iBasis=0; iBasis < numBasis; ++iBasis)
+	for (int iDim=0; iDim < spaceDim; ++iDim)
+	  tractionsCell[iQuad*spaceDim+iDim] = basis[iBasis] *
+	    solutionCell[iBasis*spaceDim+iDim] * stiffnessCell[iBasis];
+    (*tractions)->updatePoint(*c_iter, &tractionsCell[0]);
+  } // for
+
+  //solution->view("SOLUTION");
+  //(*tractions)->view("TRACTIONS");
+} // _calcTractionsChange
+
+// ----------------------------------------------------------------------
 // Allocate scalar field for output of vertex information.
 void
 pylith::faults::FaultCohesiveKin::_allocateBufferVertexScalar(void)
@@ -814,7 +864,9 @@
 pylith::faults::FaultCohesiveKin::_allocateBufferCellVector(void)
 { // _allocateBufferCellVector
   assert(0 != _quadrature);
-  const int fiberDim = _quadrature->spaceDim();
+  const int numQuadPts = _quadrature->numQuadPts();
+  const int spaceDim = _quadrature->spaceDim();
+  const int fiberDim = numQuadPts * spaceDim;
   if (_bufferCellVector.isNull()) {
     _bufferCellVector = new real_section_type(_faultMesh->comm(), 
 					      _faultMesh->debug());

Modified: short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2008-02-15 22:09:53 UTC (rev 11168)
+++ short/3D/PyLith/trunk/libsrc/faults/FaultCohesiveKin.hh	2008-02-15 22:11:44 UTC (rev 11169)
@@ -193,6 +193,14 @@
   void _calcConditioning(const spatialdata::geocoords::CoordSys* cs,
 			 spatialdata::spatialdb::SpatialDB* matDB);
 
+  /** Compute change in tractions on fault surface using solution.
+   *
+   * @param tractions Field for tractions.
+   * @param solution Solution field.
+   */
+  void _calcTractionsChange(ALE::Obj<real_section_type>* tractions,
+			    const ALE::Obj<real_section_type>& solution);
+
   /// Allocate scalar field for output of vertex information.
   void _allocateBufferVertexScalar(void);