[cig-commits] r15560 - short/3D/PyLith/branches/pylith-friction/libsrc/faults
surendra at geodynamics.org
surendra at geodynamics.org
Wed Aug 19 13:16:27 PDT 2009
Author: surendra
Date: 2009-08-19 13:16:27 -0700 (Wed, 19 Aug 2009)
New Revision: 15560
Modified:
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh
Log:
Worked on integrate residual
Modified: short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-19 17:21:51 UTC (rev 15559)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.cc 2009-08-19 20:16:27 UTC (rev 15560)
@@ -52,7 +52,22 @@
const double upDir[3],
const double normalDir[3])
{ // initialize
- throw std::logic_error("FaultCohesiveDyn::initialize() not implemented.");
+ assert(0 != upDir);
+ assert(0 != normalDir);
+ assert(0 != _quadrature);
+ assert(0 != _normalizer);
+
+ // _calcOrientation()
+ // Compute orientation at vertices in fault mesh.
+ _calcOrientation(upDir, normalDir);
+
+ // Compute tributary area for each vertex in fault mesh.
+ _calcArea();
+
+ //_getInitialTractions()
+
+
+ // _initializeConstitutiveModel()
} // initialize
// ----------------------------------------------------------------------
@@ -240,10 +255,286 @@
// ----------------------------------------------------------------------
// Calculate orientation at fault vertices.
void
-pylith::faults::FaultCohesiveDyn::_calcOrientation(const double upDir[3],
+pylith::faults::FaultCohesive::_calcOrientation(const double upDir[3],
const double normalDir[3])
{ // _calcOrientation
- throw std::logic_error("FaultCohesiveDyn::_calcOrientation() not implemented.");
+ 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/FaultCohesiveDyn.hh
===================================================================
--- short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh 2009-08-19 17:21:51 UTC (rev 15559)
+++ short/3D/PyLith/branches/pylith-friction/libsrc/faults/FaultCohesiveDyn.hh 2009-08-19 20:16:27 UTC (rev 15560)
@@ -131,6 +131,9 @@
void _calcOrientation(const double upDir[3],
const double normalDir[3]);
+ /// Calculate fault area field.
+ void _calcArea(void);
+
// NOT IMPLEMENTED ////////////////////////////////////////////////////
private :
More information about the CIG-COMMITS
mailing list