[cig-commits] r19627 - short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults
brad at geodynamics.org
brad at geodynamics.org
Tue Feb 14 14:44:57 PST 2012
Author: brad
Date: 2012-02-14 14:44:56 -0800 (Tue, 14 Feb 2012)
New Revision: 19627
Modified:
short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc
Log:
Removed unnecessary updates to dispRel. Minor cleanup of variable names.
Modified: short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc
===================================================================
--- short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-14 20:16:57 UTC (rev 19626)
+++ short/3D/PyLith/branches/pylith-scecdynrup/libsrc/pylith/faults/FaultCohesiveDyn.cc 2012-02-14 22:44:56 UTC (rev 19627)
@@ -530,7 +530,6 @@
// Allocate arrays for vertex values
scalar_array tractionTpdtVertex(spaceDim);
- scalar_array dTractionTpdtVertex(spaceDim);
scalar_array dDispRelVertex(spaceDim);
// Get sections
@@ -557,8 +556,8 @@
fields->get("dispIncr(t->t+dt)").section();
assert(!dispIncrSection.isNull());
+ scalar_array dTractionTpdtVertex(spaceDim);
scalar_array dLagrangeTpdtVertex(spaceDim);
- scalar_array dLagrangeTpdtVertexGlobal(spaceDim);
const ALE::Obj<RealSection>& dLagrangeTpdtSection =
_fields->get("sensitivity dLagrange").section();
assert(!dLagrangeTpdtSection.isNull());
@@ -694,7 +693,7 @@
} // if
// Step 2: Apply friction criterion to trial solution to get
- // change in Lagrange multiplier (dLagrangeTpdtVertex) in fault
+ // change in Lagrange multiplier (dTractionTpdtVertex) in fault
// coordinate system.
// Get friction properties and state variables.
@@ -702,25 +701,25 @@
// Use fault constitutive model to compute traction associated with
// friction.
- dLagrangeTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
const bool iterating = true; // Iterating to get friction
CALL_MEMBER_FN(*this,
- constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
+ constrainSolnSpaceFn)(&dTractionTpdtVertex,
t, slipTpdtVertex, slipRateVertex,
tractionTpdtVertex,
iterating);
// Rotate increment in traction back to global coordinate system.
- dLagrangeTpdtVertexGlobal = 0.0;
+ dLagrangeTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- dLagrangeTpdtVertexGlobal[iDim] +=
- orientationVertex[jDim*spaceDim+iDim] * dLagrangeTpdtVertex[jDim];
+ dLagrangeTpdtVertex[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] * dTractionTpdtVertex[jDim];
} // for
// Add in potential contribution from adjusting Lagrange
// multiplier for fault normal DOF of trial solution in Step 1.
- dLagrangeTpdtVertexGlobal[iDim] +=
+ dLagrangeTpdtVertex[iDim] +=
orientationVertex[indexN*spaceDim+iDim] * dTractionTpdtVertexNormal;
} // for
@@ -741,19 +740,19 @@
std::cout << ", lagrangeTIncrVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << lagrangeTIncrVertex[iDim];
+ std::cout << ", dTractionTpdtVertex: ";
+ for (int iDim=0; iDim < spaceDim; ++iDim)
+ std::cout << " " << dTractionTpdtVertex[iDim];
std::cout << ", dLagrangeTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << dLagrangeTpdtVertex[iDim];
- std::cout << ", dLagrangeTpdtVertexGlobal: ";
- for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertexGlobal[iDim];
std::cout << std::endl;
#endif
// Set change in Lagrange multiplier
- assert(dLagrangeTpdtVertexGlobal.size() ==
+ assert(dLagrangeTpdtVertex.size() ==
dLagrangeTpdtSection->getFiberDimension(v_fault));
- dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertexGlobal[0]);
+ dLagrangeTpdtSection->updatePoint(v_fault, &dLagrangeTpdtVertex[0]);
// Update displacement in trial solution (if necessary) so that it
// conforms to physical constraints.
@@ -1008,12 +1007,6 @@
std::cout << std::endl;
#endif
-#if 0 // TEST
- // Set change in relative displacement.
- assert(dispRelVertex.size() ==
- dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updatePoint(v_fault, &dispRelVertex[0]);
-#endif
// Update Lagrange multiplier increment.
assert(dLagrangeTpdtVertex.size() ==
dispIncrSection->getFiberDimension(v_lagrange));
@@ -1086,8 +1079,8 @@
// Allocate arrays for vertex values
scalar_array tractionTpdtVertex(spaceDim);
scalar_array lagrangeTpdtVertex(spaceDim);
+ scalar_array dTractionTpdtVertex(spaceDim);
scalar_array dLagrangeTpdtVertex(spaceDim);
- scalar_array dLagrangeTpdtVertexGlobal(spaceDim);
// Update time step in friction (can vary).
_friction->timeStep(_dt);
@@ -1272,20 +1265,20 @@
// Use fault constitutive model to compute traction associated with
// friction.
- dLagrangeTpdtVertex = 0.0;
+ dTractionTpdtVertex = 0.0;
const bool iterating = false; // No iteration for friction in lumped soln
CALL_MEMBER_FN(*this,
- constrainSolnSpaceFn)(&dLagrangeTpdtVertex,
+ constrainSolnSpaceFn)(&dTractionTpdtVertex,
t, slipVertex, slipRateVertex,
tractionTpdtVertex,
iterating);
// Rotate traction back to global coordinate system.
- dLagrangeTpdtVertexGlobal = 0.0;
+ dLagrangeTpdtVertex = 0.0;
for (int iDim=0; iDim < spaceDim; ++iDim) {
for (int jDim=0; jDim < spaceDim; ++jDim) {
- dLagrangeTpdtVertexGlobal[iDim] +=
- orientationVertex[jDim*spaceDim+iDim] * dLagrangeTpdtVertex[jDim];
+ dLagrangeTpdtVertex[iDim] +=
+ orientationVertex[jDim*spaceDim+iDim] * dTractionTpdtVertex[jDim];
} // for
} // for
@@ -1299,27 +1292,27 @@
std::cout << ", slipVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipVertex[iDim];
- std::cout << ", slipRateVertex: ";
+ std::cout << ", slipRateVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << slipRateVertex[iDim];
std::cout << ", orientationVertex: ";
for (int iDim=0; iDim < spaceDim*spaceDim; ++iDim)
std::cout << " " << orientationVertex[iDim];
- std::cout << ", tractionVertex: ";
+ std::cout << ", tractionVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << tractionTpdtVertex[iDim];
- std::cout << ", lagrangeTVertex: ";
+ std::cout << ", lagrangeTVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << lagrangeTVertex[iDim];
- std::cout << ", lagrangeTIncrVertex: ";
+ std::cout << ", lagrangeTIncrVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
std::cout << " " << lagrangeTIncrVertex[iDim];
- std::cout << ", dLagrangeTpdtVertex: ";
+ std::cout << ", dTractionTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertex[iDim];
- std::cout << ", dLagrangeTpdtVertexGlobal: ";
+ std::cout << " " << dTractionTpdtVertex[iDim];
+ std::cout << ", dLagrangeTpdtVertex: ";
for (int iDim=0; iDim < spaceDim; ++iDim)
- std::cout << " " << dLagrangeTpdtVertexGlobal[iDim];
+ std::cout << " " << dLagrangeTpdtVertex[iDim];
std::cout << std::endl;
#endif
@@ -1329,16 +1322,12 @@
assert(jacobianVertexN[iDim] > 0.0);
dispIncrVertexN[iDim] +=
- areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexN[iDim];
+ areaVertex * dLagrangeTpdtVertex[iDim] / jacobianVertexN[iDim];
dispIncrVertexP[iDim] -=
- areaVertex * dLagrangeTpdtVertexGlobal[iDim] / jacobianVertexP[iDim];
+ areaVertex * dLagrangeTpdtVertex[iDim] / jacobianVertexP[iDim];
- // Set increment in relative displacement.
- dispRelVertex[iDim] = -areaVertex * 2.0*dLagrangeTpdtVertexGlobal[iDim] /
- (jacobianVertexN[iDim] + jacobianVertexP[iDim]);
-
// Update increment in Lagrange multiplier.
- lagrangeTIncrVertex[iDim] += dLagrangeTpdtVertexGlobal[iDim];
+ lagrangeTIncrVertex[iDim] += dLagrangeTpdtVertex[iDim];
} // for
#if defined(DETAILED_EVENT_LOGGING)
@@ -1361,7 +1350,8 @@
} // if
// The Lagrange multiplier and relative displacement are NOT
- // assembled across processors.
+ // assembled across processors, so update even if Lagrange vertex
+ // is not local.
// Set Lagrange multiplier value. Value from preliminary solve is
// bogus due to artificial diagonal entry in Jacobian of 1.0.
@@ -1369,12 +1359,6 @@
dispIncrSection->getFiberDimension(v_lagrange));
dispIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
- // Update the relative displacement estimate based on adjustment
- // to the Lagrange multiplier values.
- assert(dispRelVertex.size() ==
- dispRelSection->getFiberDimension(v_fault));
- dispRelSection->updateAddPoint(v_fault, &dispRelVertex[0]);
-
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
#endif
More information about the CIG-COMMITS
mailing list