[cig-commits] r15556 - short/3D/PyLith/branches/pylith-friction/libsrc/faults
brad at geodynamics.org
brad at geodynamics.org
Tue Aug 18 16:00:19 PDT 2009
Author: brad
Date: 2009-08-18 16:00:18 -0700 (Tue, 18 Aug 2009)
New Revision: 15556
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.hh
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh
Log:
Moved orientation and area routines back to FaultCohesiveKin from FaultCohesive because they are vertex-based and won't be shared with FaultCohesiveDyn.
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.cc 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.cc 2009-08-18 23:00:18 UTC (rev 15556)
@@ -157,289 +157,5 @@
} // if/else
} // adjustTopology
-// ----------------------------------------------------------------------
-// Calculate orientation at fault vertices.
-void
-pylith::faults::FaultCohesive::_calcOrientation(const double upDir[3],
- const double normalDir[3])
-{ // _calcOrientation
- assert(0 != upDir);
- assert(0 != normalDir);
- assert(0 != _faultMesh);
- assert(0 != _fields);
- double_array upDirArray(upDir, 3);
-
- // Get vertices in fault mesh.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- // Containers for orientation information.
- const int cohesiveDim = _faultMesh->dimension();
- const int numBasis = _quadrature->numBasis();
- const int spaceDim = _quadrature->spaceDim();
- const int orientationSize = spaceDim*spaceDim;
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
- const double_array& verticesRef = cellGeometry.vertices();
- const int jacobianSize = (cohesiveDim > 0) ? spaceDim * cohesiveDim : 1;
- const double_array& quadWts = _quadrature->quadWts();
- double_array jacobian(jacobianSize);
- double jacobianDet = 0;
- double_array orientationVertex(orientationSize);
- double_array coordinatesCell(numBasis*spaceDim);
- double_array refCoordsVertex(cohesiveDim);
-
- // Allocate orientation field.
- _fields->add("orientation", "orientation");
- topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- orientation.newSection(slip, orientationSize);
- const ALE::Obj<RealSection>& orientationSection = orientation.section();
- assert(!orientationSection.isNull());
- // Create subspaces for along-strike, up-dip, and normal directions
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
- orientationSection->addSpace();
- for (int iDim=0; iDim <= cohesiveDim; ++iDim)
- orientationSection->setFiberDimension(vertices, spaceDim, iDim);
- orientation.allocate();
- orientation.zero();
-
- // Get fault cells.
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
- // Compute orientation of fault at constraint vertices
-
- // Get section containing coordinates of vertices
- const ALE::Obj<RealSection>& coordinatesSection =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinatesSection.isNull());
- topology::Mesh::RestrictVisitor coordinatesVisitor(*coordinatesSection,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- // Set orientation function
- assert(cohesiveDim == _quadrature->cellDim());
- assert(spaceDim == _quadrature->spaceDim());
-
- // Loop over cohesive cells, computing orientation weighted by
- // jacobian at constraint vertices
-
- const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
- assert(!sieve.isNull());
- typedef ALE::SieveAlg<SieveSubMesh> SieveAlg;
-
- ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0, faultSieveMesh->depth())));
-
- for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- // Get orientations at fault cell's vertices.
- coordinatesVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordinatesVisitor);
-
- ncV.clear();
- ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
- const int coneSize = ncV.getSize();
- const Mesh::point_type *cone = ncV.getPoints();
-
- for (int v=0; v < coneSize; ++v) {
- // Compute Jacobian and determinant of Jacobian at vertex
- memcpy(&refCoordsVertex[0], &verticesRef[v*cohesiveDim],
- cohesiveDim*sizeof(double));
- cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell,
- refCoordsVertex);
-
- // Compute orientation
- cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
- upDirArray);
-
- // Update orientation
- orientationSection->updateAddPoint(cone[v], &orientationVertex[0]);
- } // for
- } // for
-
- //orientation.view("ORIENTATION BEFORE COMPLETE");
-
- // Assemble orientation information
- orientation.complete();
-
- // Loop over vertices, make orientation information unit magnitude
- double_array vertexDir(orientationSize);
- int count = 0;
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter, ++count) {
- orientationVertex = 0.0;
- orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
- orientationVertex.size());
- for (int iDim=0; iDim < spaceDim; ++iDim) {
- double mag = 0;
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- mag += pow(orientationVertex[index+jDim],2);
- mag = sqrt(mag);
- assert(mag > 0.0);
- for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
- orientationVertex[index+jDim] /= mag;
- } // for
-
- orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
- } // for
- PetscLogFlops(count * orientationSize * 4);
-
- if (2 == cohesiveDim && vertices->size() > 0) {
- // Check orientation of first vertex, if dot product of fault
- // normal with preferred normal is negative, flip up/down dip
- // direction.
- //
- // If the user gives the correct normal direction (points from
- // footwall to ahanging wall), we should end up with
- // left-lateral-slip, reverse-slip, and fault-opening for positive
- // slip values.
- //
- // When we flip the up/down dip direction, we create a left-handed
- // strike/dip/normal coordinate system, but it gives the correct
- // sense of slip. In reality the strike/dip/normal directions that
- // are used are the opposite of what we would want, but we cannot
- // flip the fault normal direction because it is tied to how the
- // cohesive cells are created.
-
- assert(vertices->size() > 0);
- orientationSection->restrictPoint(*vertices->begin(), &orientationVertex[0],
- orientationVertex.size());
-
- assert(3 == spaceDim);
- double_array normalDirVertex(&orientationVertex[6], 3);
- const double normalDot =
- normalDir[0]*normalDirVertex[0] +
- normalDir[1]*normalDirVertex[1] +
- normalDir[2]*normalDirVertex[2];
-
- const int istrike = 0;
- const int idip = 3;
- const int inormal = 6;
- if (normalDot < 0.0) {
- // Flip dip direction
- for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
- v_iter != verticesEnd;
- ++v_iter) {
- orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
- orientationVertex.size());
- assert(9 == orientationSection->getFiberDimension(*v_iter));
- for (int iDim=0; iDim < 3; ++iDim) // flip dip
- orientationVertex[idip+iDim] *= -1.0;
-
- // Update direction
- orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
- } // for
- PetscLogFlops(5 + count * 3);
- } // if
- } // if
-
- //orientation.view("ORIENTATION");
-} // _calcOrientation
-
-// ----------------------------------------------------------------------
-void
-pylith::faults::FaultCohesive::_calcArea(void)
-{ // _calcArea
- assert(0 != _faultMesh);
- assert(0 != _fields);
-
- // Containers for area information
- const int cellDim = _quadrature->cellDim();
- const int numBasis = _quadrature->numBasis();
- const int numQuadPts = _quadrature->numQuadPts();
- const int spaceDim = _quadrature->spaceDim();
- const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
- const double_array& quadWts = _quadrature->quadWts();
- assert(quadWts.size() == numQuadPts);
- double jacobianDet = 0;
- double_array areaCell(numBasis);
-
- // Get vertices in fault mesh.
- const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
- assert(!faultSieveMesh.isNull());
- const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
- faultSieveMesh->depthStratum(0);
- const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
- const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
-
- // Allocate area field.
- _fields->add("area", "area");
-
- topology::Field<topology::SubMesh>& area = _fields->get("area");
- const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
- area.newSection(slip, 1);
- area.allocate();
- area.zero();
- const ALE::Obj<RealSection>& areaSection = area.section();
- assert(!areaSection.isNull());
- topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
-
- double_array coordinatesCell(numBasis*spaceDim);
- const ALE::Obj<RealSection>& coordinates =
- faultSieveMesh->getRealSection("coordinates");
- assert(!coordinates.isNull());
- topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
- coordinatesCell.size(),
- &coordinatesCell[0]);
-
- const ALE::Obj<SieveSubMesh::label_sequence>& cells =
- faultSieveMesh->heightStratum(0);
- assert(!cells.isNull());
- const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
- const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
-
- // Loop over cells in fault mesh, compute area
- for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
- c_iter != cellsEnd;
- ++c_iter) {
- areaCell = 0.0;
-
- // Compute geometry information for current cell
-#if defined(PRECOMPUTE_GEOMETRY)
- _quadrature->retrieveGeometry(*c_iter);
-#else
- coordsVisitor.clear();
- faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
- _quadrature->computeGeometry(coordinatesCell, *c_iter);
-#endif
-
- // Get cell geometry information that depends on cell
- const double_array& basis = _quadrature->basis();
- const double_array& jacobianDet = _quadrature->jacobianDet();
-
- // Compute area
- for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
- const double wt = quadWts[iQuad] * jacobianDet[iQuad];
- for (int iBasis=0; iBasis < numBasis; ++iBasis) {
- const double dArea = wt*basis[iQuad*numBasis+iBasis];
- areaCell[iBasis] += dArea;
- } // for
- } // for
- areaVisitor.clear();
- faultSieveMesh->updateAdd(*c_iter, areaVisitor);
-
- PetscLogFlops( numQuadPts*(1+numBasis*2) );
- } // for
-
- // Assemble area information
- area.complete();
-
-#if 0 // DEBUGGING
- area.view("AREA");
- //_faultMesh->getSendOverlap()->view("Send fault overlap");
- //_faultMesh->getRecvOverlap()->view("Receive fault overlap");
-#endif
-} // _calcArea
-
-
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.hh 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesive.hh 2009-08-18 23:00:18 UTC (rev 15556)
@@ -92,24 +92,6 @@
virtual
bool useLagrangeConstraints(void) const = 0;
- // PROTECTED METHODS //////////////////////////////////////////////////
-protected :
-
- /** Calculate orientation at fault vertices.
- *
- * @param upDir Direction perpendicular to along-strike direction that is
- * not collinear with fault normal (usually "up" direction but could
- * be up-dip direction; only applies to fault surfaces in a 3-D domain).
- * @param normalDir General preferred direction for fault normal
- * (used to pick which of two possible normal directions for
- * interface; only applies to fault surfaces in a 3-D domain).
- */
- void _calcOrientation(const double upDir[3],
- const double normalDir[3]);
-
- /// Calculate fault area field.
- void _calcArea(void);
-
// PROTECTED MEMBERS //////////////////////////////////////////////////
protected :
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-18 23:00:18 UTC (rev 15556)
@@ -144,5 +144,13 @@
throw std::logic_error("FaultCohesiveDyn::vertexField() not implemented.");
} // cellField
+// ----------------------------------------------------------------------
+// Calculate orientation at fault vertices.
+void
+pylith::faults::FaultCohesiveDyn::_calcOrientation(const double upDir[3],
+ const double normalDir[3])
+{ // _calcOrientation
+ throw std::logic_error("FaultCohesiveDyn::_calcOrientation() not implemented.");
+} // _calcOrientation
// End of file
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh 2009-08-18 23:00:18 UTC (rev 15556)
@@ -116,6 +116,21 @@
*/
bool useLagrangeConstraints(void) const;
+ // PRIVATE METHODS ////////////////////////////////////////////////////
+private :
+
+ /** Calculate orientation at quadrature points.
+ *
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void _calcOrientation(const double upDir[3],
+ const double normalDir[3]);
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.cc 2009-08-18 23:00:18 UTC (rev 15556)
@@ -835,6 +835,290 @@
} // cellField
// ----------------------------------------------------------------------
+// Calculate orientation at fault vertices.
+void
+pylith::faults::FaultCohesiveKin::_calcOrientation(const double upDir[3],
+ const double normalDir[3])
+{ // _calcOrientation
+ assert(0 != upDir);
+ assert(0 != normalDir);
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+
+ double_array upDirArray(upDir, 3);
+
+ // Get vertices in fault mesh.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Containers for orientation information.
+ const int cohesiveDim = _faultMesh->dimension();
+ const int numBasis = _quadrature->numBasis();
+ const int spaceDim = _quadrature->spaceDim();
+ const int orientationSize = spaceDim*spaceDim;
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& verticesRef = cellGeometry.vertices();
+ const int jacobianSize = (cohesiveDim > 0) ? spaceDim * cohesiveDim : 1;
+ const double_array& quadWts = _quadrature->quadWts();
+ double_array jacobian(jacobianSize);
+ double jacobianDet = 0;
+ double_array orientationVertex(orientationSize);
+ double_array coordinatesCell(numBasis*spaceDim);
+ double_array refCoordsVertex(cohesiveDim);
+
+ // Allocate orientation field.
+ _fields->add("orientation", "orientation");
+ topology::Field<topology::SubMesh>& orientation = _fields->get("orientation");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
+ orientation.newSection(slip, orientationSize);
+ const ALE::Obj<RealSection>& orientationSection = orientation.section();
+ assert(!orientationSection.isNull());
+ // Create subspaces for along-strike, up-dip, and normal directions
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ orientationSection->addSpace();
+ for (int iDim=0; iDim <= cohesiveDim; ++iDim)
+ orientationSection->setFiberDimension(vertices, spaceDim, iDim);
+ orientation.allocate();
+ orientation.zero();
+
+ // Get fault cells.
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Compute orientation of fault at constraint vertices
+
+ // Get section containing coordinates of vertices
+ const ALE::Obj<RealSection>& coordinatesSection =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinatesSection.isNull());
+ topology::Mesh::RestrictVisitor coordinatesVisitor(*coordinatesSection,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ // Set orientation function
+ assert(cohesiveDim == _quadrature->cellDim());
+ assert(spaceDim == _quadrature->spaceDim());
+
+ // Loop over cohesive cells, computing orientation weighted by
+ // jacobian at constraint vertices
+
+ const ALE::Obj<SieveSubMesh::sieve_type>& sieve = faultSieveMesh->getSieve();
+ assert(!sieve.isNull());
+ typedef ALE::SieveAlg<SieveSubMesh> SieveAlg;
+
+ ALE::ISieveVisitor::NConeRetriever<SieveMesh::sieve_type> ncV(*sieve, (size_t) pow(sieve->getMaxConeSize(), std::max(0, faultSieveMesh->depth())));
+
+ for (SieveSubMesh::label_sequence::iterator c_iter=cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ // Get orientations at fault cell's vertices.
+ coordinatesVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordinatesVisitor);
+
+ ncV.clear();
+ ALE::ISieveTraversal<SieveSubMesh::sieve_type>::orientedClosure(*sieve, *c_iter, ncV);
+ const int coneSize = ncV.getSize();
+ const Mesh::point_type *cone = ncV.getPoints();
+
+ for (int v=0; v < coneSize; ++v) {
+ // Compute Jacobian and determinant of Jacobian at vertex
+ memcpy(&refCoordsVertex[0], &verticesRef[v*cohesiveDim],
+ cohesiveDim*sizeof(double));
+ cellGeometry.jacobian(&jacobian, &jacobianDet, coordinatesCell,
+ refCoordsVertex);
+
+ // Compute orientation
+ cellGeometry.orientation(&orientationVertex, jacobian, jacobianDet,
+ upDirArray);
+
+ // Update orientation
+ orientationSection->updateAddPoint(cone[v], &orientationVertex[0]);
+ } // for
+ } // for
+
+ //orientation.view("ORIENTATION BEFORE COMPLETE");
+
+ // Assemble orientation information
+ orientation.complete();
+
+ // Loop over vertices, make orientation information unit magnitude
+ double_array vertexDir(orientationSize);
+ int count = 0;
+ for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter, ++count) {
+ orientationVertex = 0.0;
+ orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
+ orientationVertex.size());
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ double mag = 0;
+ for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
+ mag += pow(orientationVertex[index+jDim],2);
+ mag = sqrt(mag);
+ assert(mag > 0.0);
+ for (int jDim=0, index=iDim*spaceDim; jDim < spaceDim; ++jDim)
+ orientationVertex[index+jDim] /= mag;
+ } // for
+
+ orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+ } // for
+ PetscLogFlops(count * orientationSize * 4);
+
+ if (2 == cohesiveDim && vertices->size() > 0) {
+ // Check orientation of first vertex, if dot product of fault
+ // normal with preferred normal is negative, flip up/down dip
+ // direction.
+ //
+ // If the user gives the correct normal direction (points from
+ // footwall to ahanging wall), we should end up with
+ // left-lateral-slip, reverse-slip, and fault-opening for positive
+ // slip values.
+ //
+ // When we flip the up/down dip direction, we create a left-handed
+ // strike/dip/normal coordinate system, but it gives the correct
+ // sense of slip. In reality the strike/dip/normal directions that
+ // are used are the opposite of what we would want, but we cannot
+ // flip the fault normal direction because it is tied to how the
+ // cohesive cells are created.
+
+ assert(vertices->size() > 0);
+ orientationSection->restrictPoint(*vertices->begin(), &orientationVertex[0],
+ orientationVertex.size());
+
+ assert(3 == spaceDim);
+ double_array normalDirVertex(&orientationVertex[6], 3);
+ const double normalDot =
+ normalDir[0]*normalDirVertex[0] +
+ normalDir[1]*normalDirVertex[1] +
+ normalDir[2]*normalDirVertex[2];
+
+ const int istrike = 0;
+ const int idip = 3;
+ const int inormal = 6;
+ if (normalDot < 0.0) {
+ // Flip dip direction
+ for (SieveSubMesh::label_sequence::iterator v_iter=verticesBegin;
+ v_iter != verticesEnd;
+ ++v_iter) {
+ orientationSection->restrictPoint(*v_iter, &orientationVertex[0],
+ orientationVertex.size());
+ assert(9 == orientationSection->getFiberDimension(*v_iter));
+ for (int iDim=0; iDim < 3; ++iDim) // flip dip
+ orientationVertex[idip+iDim] *= -1.0;
+
+ // Update direction
+ orientationSection->updatePoint(*v_iter, &orientationVertex[0]);
+ } // for
+ PetscLogFlops(5 + count * 3);
+ } // if
+ } // if
+
+ //orientation.view("ORIENTATION");
+} // _calcOrientation
+
+// ----------------------------------------------------------------------
+void
+pylith::faults::FaultCohesiveKin::_calcArea(void)
+{ // _calcArea
+ assert(0 != _faultMesh);
+ assert(0 != _fields);
+
+ // Containers for area information
+ const int cellDim = _quadrature->cellDim();
+ const int numBasis = _quadrature->numBasis();
+ const int numQuadPts = _quadrature->numQuadPts();
+ const int spaceDim = _quadrature->spaceDim();
+ const feassemble::CellGeometry& cellGeometry = _quadrature->refGeometry();
+ const double_array& quadWts = _quadrature->quadWts();
+ assert(quadWts.size() == numQuadPts);
+ double jacobianDet = 0;
+ double_array areaCell(numBasis);
+
+ // Get vertices in fault mesh.
+ const ALE::Obj<SieveSubMesh>& faultSieveMesh = _faultMesh->sieveMesh();
+ assert(!faultSieveMesh.isNull());
+ const ALE::Obj<SieveSubMesh::label_sequence>& vertices =
+ faultSieveMesh->depthStratum(0);
+ const SieveSubMesh::label_sequence::iterator verticesBegin = vertices->begin();
+ const SieveSubMesh::label_sequence::iterator verticesEnd = vertices->end();
+
+ // Allocate area field.
+ _fields->add("area", "area");
+
+ topology::Field<topology::SubMesh>& area = _fields->get("area");
+ const topology::Field<topology::SubMesh>& slip = _fields->get("slip");
+ area.newSection(slip, 1);
+ area.allocate();
+ area.zero();
+ const ALE::Obj<RealSection>& areaSection = area.section();
+ assert(!areaSection.isNull());
+ topology::Mesh::UpdateAddVisitor areaVisitor(*areaSection, &areaCell[0]);
+
+ double_array coordinatesCell(numBasis*spaceDim);
+ const ALE::Obj<RealSection>& coordinates =
+ faultSieveMesh->getRealSection("coordinates");
+ assert(!coordinates.isNull());
+ topology::Mesh::RestrictVisitor coordsVisitor(*coordinates,
+ coordinatesCell.size(),
+ &coordinatesCell[0]);
+
+ const ALE::Obj<SieveSubMesh::label_sequence>& cells =
+ faultSieveMesh->heightStratum(0);
+ assert(!cells.isNull());
+ const SieveSubMesh::label_sequence::iterator cellsBegin = cells->begin();
+ const SieveSubMesh::label_sequence::iterator cellsEnd = cells->end();
+
+ // Loop over cells in fault mesh, compute area
+ for (SieveSubMesh::label_sequence::iterator c_iter = cellsBegin;
+ c_iter != cellsEnd;
+ ++c_iter) {
+ areaCell = 0.0;
+
+ // Compute geometry information for current cell
+#if defined(PRECOMPUTE_GEOMETRY)
+ _quadrature->retrieveGeometry(*c_iter);
+#else
+ coordsVisitor.clear();
+ faultSieveMesh->restrictClosure(*c_iter, coordsVisitor);
+ _quadrature->computeGeometry(coordinatesCell, *c_iter);
+#endif
+
+ // Get cell geometry information that depends on cell
+ const double_array& basis = _quadrature->basis();
+ const double_array& jacobianDet = _quadrature->jacobianDet();
+
+ // Compute area
+ for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
+ const double wt = quadWts[iQuad] * jacobianDet[iQuad];
+ for (int iBasis=0; iBasis < numBasis; ++iBasis) {
+ const double dArea = wt*basis[iQuad*numBasis+iBasis];
+ areaCell[iBasis] += dArea;
+ } // for
+ } // for
+ areaVisitor.clear();
+ faultSieveMesh->updateAdd(*c_iter, areaVisitor);
+
+ PetscLogFlops( numQuadPts*(1+numBasis*2) );
+ } // for
+
+ // Assemble area information
+ area.complete();
+
+#if 0 // DEBUGGING
+ area.view("AREA");
+ //_faultMesh->getSendOverlap()->view("Send fault overlap");
+ //_faultMesh->getRecvOverlap()->view("Receive fault overlap");
+#endif
+} // _calcArea
+
+// ----------------------------------------------------------------------
// Compute change in tractions on fault surface using solution.
// NOTE: We must convert vertex labels to fault vertex labels
void
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh 2009-08-18 22:47:57 UTC (rev 15555)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveKin.hh 2009-08-18 23:00:18 UTC (rev 15556)
@@ -206,6 +206,21 @@
// PRIVATE METHODS ////////////////////////////////////////////////////
private :
+ /** Calculate orientation at fault vertices.
+ *
+ * @param upDir Direction perpendicular to along-strike direction that is
+ * not collinear with fault normal (usually "up" direction but could
+ * be up-dip direction; only applies to fault surfaces in a 3-D domain).
+ * @param normalDir General preferred direction for fault normal
+ * (used to pick which of two possible normal directions for
+ * interface; only applies to fault surfaces in a 3-D domain).
+ */
+ void _calcOrientation(const double upDir[3],
+ const double normalDir[3]);
+
+ /// Calculate fault area field.
+ void _calcArea(void);
+
/** Compute change in tractions on fault surface using solution.
*
* @param tractions Field for tractions.
More information about the CIG-COMMITS
mailing list