[cig-commits] r18971 - in short/3D/PyLith/branches/v1.6-revisedfault: libsrc/pylith/faults unittests/libtests/faults unittests/libtests/faults/data
brad at geodynamics.org
brad at geodynamics.org
Mon Sep 26 13:19:26 PDT 2011
Author: brad
Date: 2011-09-26 13:19:26 -0700 (Mon, 26 Sep 2011)
New Revision: 18971
Modified:
short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/adjustsoln.py
Log:
Updated FaultCohesiveKin::adjustSolnLumped() and corresponding test data.
Modified: short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/libsrc/pylith/faults/FaultCohesiveLagrange.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -1295,31 +1295,28 @@
const topology::Field<
topology::Mesh>& jacobian)
{ // adjustSolnLumped
- /// Member prototype for _adjustSolnLumpedXD()
- typedef void (pylith::faults::FaultCohesiveLagrange::*adjustSolnLumped_fn_type)
- (double_array*, double_array*, double_array*,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
- const double_array&, const double_array&,
- const double_array&, const double_array&);
-
assert(0 != fields);
assert(0 != _quadrature);
- // Cohesive cells with conventional vertices i and j, and constraint
- // vertex k require 2 adjustments to the solution:
+ // Cohesive cells with conventional vertices N and P, and constraint
+ // vertex L require 2 adjustments to the solution:
//
- // * DOF k: Compute increment in Lagrange multipliers
- // dl_k = S^{-1} (-C_ki (A_i^{-1} r_i - C_kj A_j^{-1} r_j + u_i - u_j) - d_k)
- // S = C_ki (A_i^{-1} + A_j^{-1}) C_ki^T
+ // * DOF L: Compute increment in Lagrange multipliers
//
+ // d\vec{l}_p = \tensor{S}^{-1} \cdot (-\vec{r}_p^{*} + \tensor{L}_p
+ // \cdot (d\vec{u}_{n+}^* - d\vec{u}_{n-}^*))
+ // \tensor{S} = \tensor{L}_p \cdot
+ // (\tensor{K}_{n+n+}{-1} + \tensor{K}_{n-n-})
+ // \cdot \tensor{L}_p^T
+ //
// * DOF i and j: Adjust displacement increment (solution) to create slip
// consistent with Lagrange multiplier constraints
- // du_i = +A_i^-1 C_ki^T dlk
- // du_j = -A_j^-1 C_kj^T dlk
+ //
+ // d\vec{u}_{n+} = d\vec{u}_{n+}^{*} - \tensor{K}_{n+n+}^{-1} \cdot
+ // \tensor{L}_p^T \cdot d\vec{l}_p
+ // d\vec{u}_{n-} = d\vec{u}_{n-}^{*} + \tensor{K}_{n-n-}^{-1} \cdot
+ // \tensor{L}_p^T \cdot d\vec{l}_p
- // :TODO: FIX THIS
-
const int setupEvent = _logger->eventId("FaAS setup");
const int geometryEvent = _logger->eventId("FaAS geometry");
const int computeEvent = _logger->eventId("FaAS compute");
@@ -1332,29 +1329,16 @@
const int spaceDim = _quadrature->spaceDim();
// Get section information
- double_array slipVertex(spaceDim);
- const ALE::Obj<RealSection>& slipSection = _fields->get("slip").section();
- assert(!slipSection.isNull());
-
- double_array jacobianVertexN(spaceDim);
- double_array jacobianVertexP(spaceDim);
const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
assert(!jacobianSection.isNull());
- double_array residualVertexN(spaceDim);
- double_array residualVertexP(spaceDim);
const ALE::Obj<RealSection>& residualSection =
fields->get("residual").section();
assert(!residualSection.isNull());
- double_array dispTVertexN(spaceDim);
- double_array dispTVertexP(spaceDim);
- const ALE::Obj<RealSection>& dispTSection = fields->get("disp(t)").section();
- assert(!dispTSection.isNull());
-
double_array dispTIncrVertexN(spaceDim);
double_array dispTIncrVertexP(spaceDim);
- double_array lagrangeTIncrVertex(spaceDim);
+ double_array dispTIncrVertexL(spaceDim); // Lagrange multiplier
const ALE::Obj<RealSection>& dispTIncrSection = fields->get(
"dispIncr(t->t+dt)").section();
assert(!dispTIncrSection.isNull());
@@ -1366,29 +1350,10 @@
const ALE::Obj<SieveMesh>& sieveMesh = fields->mesh().sieveMesh();
assert(!sieveMesh.isNull());
const ALE::Obj<SieveMesh::order_type>& globalOrder =
- sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default", jacobianSection);
+ sieveMesh->getFactory()->getGlobalOrder(sieveMesh, "default",
+ jacobianSection);
assert(!globalOrder.isNull());
- adjustSolnLumped_fn_type adjustSolnLumpedFn;
- switch (spaceDim) { // switch
- case 1:
- adjustSolnLumpedFn =
- &pylith::faults::FaultCohesiveLagrange::_adjustSolnLumped1D;
- break;
- case 2:
- adjustSolnLumpedFn =
- &pylith::faults::FaultCohesiveLagrange::_adjustSolnLumped2D;
- break;
- case 3:
- adjustSolnLumpedFn =
- &pylith::faults::FaultCohesiveLagrange::_adjustSolnLumped3D;
- break;
- default :
- assert(0);
- throw std::logic_error("Unknown spatial dimension in "
- "FaultCohesiveLagrange::adjustSolnLumped.");
- } // switch
-
_logger->eventEnd(setupEvent);
#if !defined(DETAILED_EVENT_LOGGING)
@@ -1410,42 +1375,66 @@
_logger->eventBegin(restrictEvent);
#endif
- // Get slip at fault cell's vertices.
- slipSection->restrictPoint(v_fault, &slipVertex[0], slipVertex.size());
-
// Get residual at cohesive cell's vertices.
- residualSection->restrictPoint(v_negative, &residualVertexN[0],
- residualVertexN.size());
- residualSection->restrictPoint(v_positive, &residualVertexP[0],
- residualVertexP.size());
-
+ assert(spaceDim == residualSection->getFiberDimension(v_lagrange));
+ const double* residualVertexL = residualSection->restrictPoint(v_lagrange);
+ assert(residualVertexL);
+
// Get jacobian at cohesive cell's vertices.
- jacobianSection->restrictPoint(v_negative, &jacobianVertexN[0],
- jacobianVertexN.size());
- jacobianSection->restrictPoint(v_positive, &jacobianVertexP[0],
- jacobianVertexP.size());
+ assert(spaceDim == jacobianSection->getFiberDimension(v_negative));
+ const double* jacobianVertexN = jacobianSection->restrictPoint(v_negative);
+ assert(jacobianVertexN);
- // Get disp(t) at cohesive cell's vertices.
- dispTSection->restrictPoint(v_negative, &dispTVertexN[0],
- dispTVertexN.size());
- dispTSection->restrictPoint(v_positive, &dispTVertexP[0],
- dispTVertexP.size());
+ assert(spaceDim == jacobianSection->getFiberDimension(v_positive));
+ const double* jacobianVertexP = jacobianSection->restrictPoint(v_positive);
+ assert(jacobianVertexP);
+ assert(spaceDim == jacobianSection->getFiberDimension(v_lagrange));
+ const double* jacobianVertexL = jacobianSection->restrictPoint(v_lagrange);
+ assert(jacobianVertexL);
+
+ // Get dispIncr(t) at cohesive cell's vertices.
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_negative));
+ dispTIncrSection->restrictPoint(v_negative, &dispTIncrVertexN[0],
+ dispTIncrVertexN.size());
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_positive));
+ dispTIncrSection->restrictPoint(v_positive, &dispTIncrVertexP[0],
+ dispTIncrVertexP.size());
+
+ assert(spaceDim == dispTIncrSection->getFiberDimension(v_lagrange));
+ dispTIncrSection->restrictPoint(v_lagrange, &dispTIncrVertexL[0],
+ dispTIncrVertexL.size());
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(restrictEvent);
_logger->eventBegin(computeEvent);
#endif
-#if 0
- CALL_MEMBER_FN(*this,
- adjustSolnLumpedFn)(&lagrangeTIncrVertex,
- &dispTIncrVertexN, &dispTIncrVertexP,
- slipVertex, orientationVertex,
- dispTVertexN, dispTVertexP,
- residualVertexN, residualVertexP,
- jacobianVertexN, jacobianVertexP);
-#endif
+ for (int iDim=0; iDim < spaceDim; ++iDim) {
+ assert(jacobianVertexL[iDim] > 0.0);
+ const double Sinv = (jacobianVertexP[iDim] + jacobianVertexN[iDim]) /
+ (jacobianVertexL[iDim]*jacobianVertexL[iDim]);
+ dispTIncrVertexL[iDim] = Sinv *
+ (-residualVertexL[iDim] +
+ jacobianVertexL[iDim] *
+ (dispTIncrVertexP[iDim] - dispTIncrVertexN[iDim]));
+ assert(jacobianVertexN[iDim] > 0.0);
+ dispTIncrVertexN[iDim] =
+ +jacobianVertexL[iDim]/jacobianVertexN[iDim]*dispTIncrVertexL[iDim];
+
+ assert(jacobianVertexP[iDim] > 0.0);
+ dispTIncrVertexP[iDim] =
+ -jacobianVertexL[iDim]/jacobianVertexP[iDim]*dispTIncrVertexL[iDim];
+
+ std::cout << "iDim: " << iDim
+ << ", jacobianL: " << jacobianVertexL[iDim]
+ << ", Sinv: " << Sinv
+ << std::endl;
+
+ } // for
+
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(computeEvent);
_logger->eventBegin(updateEvent);
@@ -1462,10 +1451,10 @@
dispTIncrAdjSection->updateAddPoint(v_positive, &dispTIncrVertexP[0]);
// Set Lagrange multiplier value. Value from preliminary solve is
- // bogus due to artificial diagonal entry of 1.0.
- assert(lagrangeTIncrVertex.size() ==
+ // bogus due to artificial diagonal entry.
+ assert(dispTIncrVertexL.size() ==
dispTIncrSection->getFiberDimension(v_lagrange));
- dispTIncrSection->updatePoint(v_lagrange, &lagrangeTIncrVertex[0]);
+ dispTIncrSection->updatePoint(v_lagrange, &dispTIncrVertexL[0]);
#if defined(DETAILED_EVENT_LOGGING)
_logger->eventEnd(updateEvent);
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/TestFaultCohesiveKin.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -586,13 +586,13 @@
fields.get("dispIncr adjust");
solution += dispIncrAdj;
- //solution.view("SOLUTION AFTER ADJUSTMENT"); // DEBUGGING
+ solution.view("SOLUTION AFTER ADJUSTMENT"); // DEBUGGING
const ALE::Obj<RealSection>& solutionSection = solution.section();
CPPUNIT_ASSERT(!solutionSection.isNull());
int i = 0;
- const double tolerance = 1.0e-06;
+ const double tolerance = 2.0e-06;
const double* solutionE = _data->fieldIncrAdjusted;
for (SieveMesh::label_sequence::iterator v_iter=verticesBegin;
v_iter != verticesEnd;
@@ -601,13 +601,15 @@
CPPUNIT_ASSERT_EQUAL(spaceDim, fiberDim);
const double* solutionVertex = solutionSection->restrictPoint(*v_iter);
CPPUNIT_ASSERT(0 != solutionVertex);
- for (int iDim=0; iDim < spaceDim; ++iDim, ++i)
+ for (int iDim=0; iDim < spaceDim; ++iDim, ++i) {
+ std::cout << "valE: " << solutionE[i] << ", val: " << solutionVertex[iDim] << std::endl;
if (0.0 != solutionE[i])
CPPUNIT_ASSERT_DOUBLES_EQUAL(1.0, solutionVertex[iDim]/solutionE[i],
tolerance);
else
CPPUNIT_ASSERT_DOUBLES_EQUAL(solutionE[i], solutionVertex[iDim],
tolerance);
+ } // for
} // for
} // testAdjustSolnLumped
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataHex8.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataHex8.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -1827,22 +1827,22 @@
3.2, 4.2, 5.2,
3.3, 4.3, 5.3,
3.4, 4.4, 5.4,
- 10.1749505588, 9.31910094671, 1.19608248275, // 6
- 10.1039043952, 9.3565301613, 1.55916596774, // 7
- 10.0557317781, 9.37031383702, 1.88484239309, // 8
- 9.46707670425, 8.80701959867, 2.73150557445, // 9
+ 1.205353233, -0.669560195189, 2.4212859344, // 6
+ 1.31708438297, -0.532975672131, 2.53979527703, // 7
+ 1.44967569358, -0.34540170722, 2.68738185699, // 8
+ 1.6022763433, -0.117061939422, 2.86063862983, // 9
3.9, 4.9, 5.9,
3.0, 4.0, 5.0,
3.1, 4.1, 5.1,
3.2, 4.2, 5.2,
- -4.40186602934, -1.26050109236, 10.2660586738, // 14
- -3.43749802151, -0.573632172048, 9.81022296774, // 15
- -2.65573177809, 0.0296861629776, 9.51515760691, // 16
- -1.46880950929, 1.10387616968, 8.80699471894, // 17
- -7.22865142007, -6.45587627588, -10.0124258381, // 18
- -7.61044825807, -6.46533445161, -10.4062470323, // 19
- -7.93953352294, -6.48576793174, -10.8047440227, // 20
- -7.21263527761, -5.52328996599, -10.2007380677, // 21
+ 5.94766934654, 10.2648771483, 8.85236238339, // 14
+ 5.9351099915, 9.97517405027, 8.76421837117, // 15
+ 5.95032430642, 9.74540170722, 8.71261814301, // 16
+ 5.98205399056, 9.55826920577, 8.68465814016, // 17
+ -3.4419701505, -7.75434029278, -4.6180710984, // 18
+ -3.65266498725, -8.21276107541, -4.89632755675, // 19
+ -3.82555132091, -8.57718290227, -5.12145084311, // 20
+ -3.95590258206, -8.85071149096, -5.29085046631, // 21
};
pylith::faults::CohesiveKinDataHex8::CohesiveKinDataHex8(void)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataLine2.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -154,10 +154,10 @@
const double pylith::faults::CohesiveKinDataLine2::_fieldIncrAdjusted[] = {
1.1,
- -3.09368089375, // 3
+ 5.16324319611, // 3
1.3,
- 5.33587415261, // 5
- -9.44609796626, // 6
+ -2.23297292977, // 5
+ 8.719135031441999, // 6
};
pylith::faults::CohesiveKinDataLine2::CohesiveKinDataLine2(void)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataQuad4.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -404,14 +404,14 @@
const double pylith::faults::CohesiveKinDataQuad4::_fieldIncrAdjusted[] = {
3.1, 4.1,
3.2, 4.2,
- -4.09563227733, -3.248612969, // 4
- -2.78814167707, -2.11773476302, // 5
+ 3.45041520801, 8.42020518055, // 4
+ 3.78935267675, 8.84420957728, // 5
3.5, 4.5,
3.6, 4.6,
- 9.35548350619, 10.472468741, // 8
- 8.459683341, 9.70254140433, // 9
- -9.8131968597, -9.61432196053, // 10
- -9.12482866822, -8.66339834789, // 11
+ 3.58497660564, 1.54925486193, // 8
+ 3.61310855397, 1.62531925885, // 9
+ 0.19553977041, 5.35626673472, // 10
+ 0.545093747451, 6.22189340819, // 11
};
pylith::faults::CohesiveKinDataQuad4::CohesiveKinDataQuad4(void)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTet4.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTet4.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -587,16 +587,16 @@
const double pylith::faults::CohesiveKinDataTet4::_fieldIncrAdjusted[] = {
3.1, 4.1, 5.1,
- -4.70012229942, -3.03138431537, -1.50390580878, // 3
- -4.10716666377, -2.47185604467, -0.927412749378, // 4
- -3.34252157978, -1.47021951405, -0.324816213039, // 5
+ 4.05776047684, 8.66419800363, 4.4475545411, // 3
+ 4.41507389392, 9.10077378393, 4.76947881673, // 4
+ 2.65876279931, 6.30839037001, 3.47575435673, // 5
3.5, 4.5, 5.5,
- 9.52509172457, 10.0235382365, 10.6279293566, // 7
- 9.14962036828, 9.69078492115, 10.2975758746, // 8
- 12.4395302117, 12.2183073197, 13.0147426983, // 9
- -8.67766117844, -8.04468697053, -9.48014675931, // 10
- -8.80341285807, -8.09563657419, -9.62931666291, // 11
- -8.21830731966, -8.01474269826, -9.43953021169, // 12
+ 2.95667964237, 1.25185149728, 6.16433409418, // 7
+ 2.99466885439, 1.33277448939, 6.18315418792, // 8
+ 4.03773208096, 1.32825348198, 7.69394390058, // 9
+ 3.08793771664, 16.0711128131, -2.70880365204, // 10
+ 4.34878818628, 18.7230177573, -2.06903261476, // 11
+ -3.11319624288, 8.01523955405, -8.08183170175, // 12
};
pylith::faults::CohesiveKinDataTet4::CohesiveKinDataTet4(void)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/CohesiveKinDataTri3.cc 2011-09-26 20:19:26 UTC (rev 18971)
@@ -325,13 +325,13 @@
const double pylith::faults::CohesiveKinDataTri3::_fieldIncrAdjusted[] = {
3.1, 4.1,
- 11.4124508246, 12.5863689652, // 3
- 11.1489656107, 12.3019463023, // 4
+ 2.82834517219, -0.0422000510363, // 3
+ 3.14958479199, 0.17979481945, // 4
3.4, 4.4,
- -3.06996065966, -2.20909517212, // 6
- -2.30215017286, -1.41913540767, // 7
- -10.0636427582, -9.85494098949, // 8
- -10.402530193, -10.2036552939, // 9
+ 3.79732386225, 7.89376004083, // 6
+ 3.81502339436, 7.85074513807, // 7
+ -0.445985793369, -5.09064006124, // 8
+ -0.19553977041, -5.35626673471, // 9
};
pylith::faults::CohesiveKinDataTri3::CohesiveKinDataTri3(void)
Modified: short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/adjustsoln.py
===================================================================
--- short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/adjustsoln.py 2011-09-23 20:29:51 UTC (rev 18970)
+++ short/3D/PyLith/branches/v1.6-revisedfault/unittests/libtests/faults/data/adjustsoln.py 2011-09-26 20:19:26 UTC (rev 18971)
@@ -3,358 +3,205 @@
# ----------------------------------------------------------------------
if test == "line2":
- C = 1.0
- dk = 1.89546413727; l = 1.5
- Ai = 2.2; ri = 7.5; ui = 7.2; dui = 1.2
- Aj = 2.4; rj = -7.5; uj = 7.4; duj = 1.4
+ jL = 1.0; rL = -1.69546413727;
+ jN = 2.2; duN = 1.2
+ jP = 2.4; duP = 1.4
- Si = (Ai * Aj) / (Ai + Aj)
- Aru = ri/Ai - rj/Aj + ui - uj
- Aruslip = -C*Aru - dk
- dlp = Si * Aruslip
+ Sinv = (jN + jP) /jL**2
+ duL = Sinv * (-rL + jL*(duP-duN));
- ddui = +C / Ai * dlp
- dduj = -C / Aj * dlp
+ dduN = jL / jN * duL
+ dduP = -jL / jP * duL
#print "Aru",Aru
#print "Aruslip",Aruslip
#print "Si",Si
#print "dlip",dlp
- print dui+ddui,duj+dduj,dlp
+ print duN+dduN,duP+dduP,duL
# ----------------------------------------------------------------------
elif test == "tri3" or test == "quad4":
if test == "tri3":
- Cpx = 0.0
- Cpy = -1.0
- Cqx = -1.0
- Cqy = 0.0
-
if vertex == 0:
- # Lagrange vertex 8, vertex i: 3, vertex j: 6
- dkp = 1.89546413727; lp = 3.6
- dkq = 0.08241148423; lq = 4.6
-
- # vertex i
- Aix = 1.2; Aiy = 1.2
- rix = -9.6; riy = -8.6
- uix = 8.2; uiy = 9.2
- duix = 3.2; duiy = 4.2
-
- # vertex j
- Ajx = 1.5; Ajy = 1.5
- rjx = +9.6; rjy = +8.6
- ujx = 8.5; ujy = 9.5
- dujx = 3.5; dujy = 4.5
-
+ # Lagrange vertex 8, vertex N: 3, vertex P: 6
+ jL = 1.0
+ rLx = 0.5*(8.5-8.2 + 8.7-8.3) + 0.5*(0.08241148423+0.14794836271)
+ rLy = 0.5*(9.5-9.2 + 9.7-9.3) + 0.5*(1.89546413727+1.77538035254)
+ jN = 1.2; duNx = 3.2; duNy = 4.2;
+ jP = 1.5; duPx = 3.5; duPy = 4.5;
+
elif vertex == 1:
- # Lagrange vertex 9, vertex i: 4, vertex j: 7
- dkp = 1.77538035254; lp = 3.8
- dkq = 0.14794836271; lq = 4.8
+ # Lagrange vertex 9, vertex N: 4, vertex P: 7
+ jL = 1.0
+ rLx = 0.5*(8.5-8.2 + 8.7-8.3) + 0.5*(0.08241148423+0.14794836271)
+ rLy = 0.5*(9.5-9.2 + 9.7-9.3) + 0.5*(1.89546413727+1.77538035254)
+ jN = 1.3; duNx = 3.3; duNy = 4.3;
+ jP = 1.7; duPx = 3.7; duPy = 4.7;
- # vertex i
- Aix = 1.3; Aiy = 1.3
- rix = -9.8; riy = -8.8
- uix = 8.3; uiy = 9.3
- duix = 3.3; duiy = 4.3
-
- # vertex j
- Ajx = 1.7; Ajy = 1.7
- rjx = +9.8; rjy = +8.8
- ujx = 8.7; ujy = 9.7
- dujx = 3.7; dujy = 4.7
-
elif test == "quad4":
- Cpx = 0.0
- Cpy = 1.0
- Cqx = 1.0
- Cqy = 0.0
-
if vertex == 0:
- # Lagrange vertex 10, vertex i: 4, vertex j: 8
- dkp = 1.77538035254; lp = 3.8
- dkq = 0.14794836271; lq = 4.8
+ # Lagrange vertex 10, vertex N: 4, vertex P: 8
+ jL = 1.0
+ rLx = 0.5*(8.7-8.3 + 8.9-8.4) - 0.5*(0.14794836271+0.08241148423)
+ rLy = 0.5*(9.7-9.3 + 9.9-9.4) - 0.5*(1.77538035254+1.89546413727)
+ jN = 1.3; duNx = 3.3; duNy = 4.3;
+ jP = 1.7; duPx = 3.7; duPy = 4.7;
- # vertex i: 4
- Aix = 1.3; Aiy = 1.3
- rix = +9.8; riy = +8.8
- uix = 8.3; uiy = 9.3
- duix = 3.3; duiy = 4.3
-
- # vertex j: 8
- Ajx = 1.7; Ajy = 1.7
- rjx = -9.8; rjy = -8.8
- ujx = 8.7; ujy = 9.7
- dujx = 3.7; dujy = 4.7
-
elif vertex == 1:
- # Lagrange vertex 11, vertex i: 5, vertex j: 9
- dkp = 1.89546413727; lp = 3.0
- dkq = 0.08241148423; lq = 4.0
-
- # vertex i: 5
- Aix = 1.4; Aiy = 1.4
- rix = +9.0; riy = +8.0
- uix = 8.4; uiy = 9.4
- duix = 3.4; duiy = 4.4
-
- # vertex j: 9
- Ajx = 1.9; Ajy = 1.9
- rjx = -9.0; rjy = -8.0
- ujx = 8.9; ujy = 9.9
- dujx = 3.9; dujy = 4.9
+ # Lagrange vertex 11, vertex N: 5, vertex P: 9
+ jL = 1.0
+ rLx = 0.5*(8.7-8.3 + 8.9-8.4) - 0.5*(0.14794836271+0.08241148423)
+ rLy = 0.5*(9.7-9.3 + 9.9-9.4) - 0.5*(1.77538035254+1.89546413727)
+ jN = 1.4; duNx = 3.4; duNy = 4.4;
+ jP = 1.9; duPx = 3.9; duPy = 4.9;
- Sppi = Aix*Ajx / (Aix + Ajx)
- Sqqi = Aix*Ajx / (Aix + Ajx)
+ Sinv = (jN + jP) / (jL*jL)
+ duLx = Sinv * (-rLx + jL*(duPx - duNx))
+ duLy = Sinv * (-rLy + jL*(duPy - duNy))
- Arux = rix / Aix - rjx / Ajx + uix - ujx
- Aruy = riy / Aiy - rjy / Ajy + uiy - ujy
+ dduNx = jL / jN * duLx
+ dduNy = jL / jN * duLy
- Arup = Cpx*Arux + Cpy*Aruy
- Aruq = Cqx*Arux + Cqy*Aruy
- Arupslip = -Arup - dkp
- Aruqslip = -Aruq - dkq
+ dduPx = -jL / jP * duLx
+ dduPy = -jL / jP * duLy
+
+ print duNx+dduNx,duNy+dduNy
+ print duPx+dduPx,duPy+dduPy
+ print duLx,duLy
- dlp = Sppi * Arupslip
- dlq = Sqqi * Aruqslip
-
- dduix = +1.0/Aix * (Cpx*dlp + Cqx*dlq)
- dduiy = +1.0/Aiy * (Cpy*dlp + Cqy*dlq)
-
- ddujx = -1.0/Ajx * (Cpx*dlp + Cqx*dlq)
- ddujy = -1.0/Ajy * (Cpy*dlp + Cqy*dlq)
-
- print "Sppi",Sppi
- print "Sqqi",Sqqi
-
- print "Arup",Arup
- print "Aruq",Aruq
-
- print "Arupslip",Arupslip
- print "Aruqslip",Aruqslip
-
- print "dlp",dlp
- print "dlq",dlq
-
- print "dduix:",dduix
- print "dduiy:",dduiy
-
- print "ddujx:",ddujx
- print "ddujy:",ddujy
-
- print duix+dduix,duiy+dduiy
- print dujx+ddujx,dujy+ddujy
- print dlp,dlq
-
# ----------------------------------------------------------------------
elif test == "tet4" or test == "hex8":
if test == "tet4":
- Cpx = 0.0
- Cpy = 1.0
- Cpz = 0.0
- Cqx = 0.0
- Cqy = 0.0
- Cqz = 1.0
- Crx = 1.0
- Cry = 0.0
- Crz = 0.0
-
if vertex == 0:
- # Lagrange vertex 10, vertex i: 3, vertex j: 7
- dkp = 1.82575588523; lp = 3.7
- dkq = -0.55566483464; lq = 4.7
- dkr = 0.07938069066; lr = 5.7
-
- # vertex i: 3
- Aix = 1.2; Aiy = 1.2; Aiz = 1.2
- rix = +9.7; riy = +7.7; riz = +8.7
- uix = 7.2; uiy = 8.2; uiz = 9.2
- duix = 3.2; duiy = 4.2; duiz = 5.2
-
- # vertex j: 7
- Ajx = 1.6; Ajy = 1.6; Ajz = 1.6
- rjx = -9.7; rjy = -7.7; rjz = -8.7
- ujx = 7.6; ujy = 8.6; ujz = 9.6
- dujx = 3.6; dujy = 4.6; dujz = 5.6
-
+ # Lagrange vertex 10, vertex N: 3, vertex P: 7
+ jL = 1.0/3.0
+ rLx = (7.6-7.2 + 7.8-7.3 + 7.0-7.4)/9.0 + \
+ (-0.07938069066 + -0.14140241667 + -0.18205179147)/9.0
+ rLy = (8.6-8.2 + 8.8-8.3 + 9.0-9.4)/9.0 + \
+ (-1.69682900001 + -1.82575588523 + -1.51709826228)/9.0
+ rLz = (9.6-9.2 + 9.8-9.3 + 9.0-9.4)/9.0 + \
+ (0.55566483464 + 0.56560966667 + 0.54615537442)/9.0
+ jN = 1.2; duNx = 3.2; duNy = 4.2; duNz = 5.2;
+ jP = 1.6; duPx = 3.6; duPy = 4.6; duPz = 5.6;
+
elif vertex == 1:
- # Lagrange vertex 11, vertex i: 4, vertex j: 8
- dkp = 1.69682900001; lp = 3.9
- dkq = -0.56560966667; lq = 4.9
- dkr = 0.14140241667; lr = 5.9
+ # Lagrange vertex 11, vertex N: 4, vertex P: 8
+ jL = 1.0/3.0
+ rLx = (7.6-7.2 + 7.8-7.3 + 7.0-7.4)/9.0 + \
+ (-0.07938069066 + -0.14140241667 + -0.18205179147)/9.0
+ rLy = (8.6-8.2 + 8.8-8.3 + 9.0-9.4)/9.0 + \
+ (-1.69682900001 + -1.82575588523 + -1.51709826228)/9.0
+ rLz = (9.6-9.2 + 9.8-9.3 + 9.0-9.4)/9.0 + \
+ (0.55566483464 + 0.56560966667 + 0.54615537442)/9.0
+ jN = 1.3; duNx = 3.3; duNy = 4.3; duNz = 5.3;
+ jP = 1.8; duPx = 3.8; duPy = 4.8; duPz = 5.8;
- # vertex i: 4
- Aix = 1.3; Aiy = 1.3; Aiz = 1.3
- rix = +9.9; riy = +7.9; riz = +8.9
- uix = 7.3; uiy = 8.3; uiz = 9.3
- duix = 3.3; duiy = 4.3; duiz = 5.3
-
- # vertex j: 8
- Ajx = 1.8; Ajy = 1.8; Ajz = 1.8
- rjx = -9.9; rjy = -7.9; rjz = -8.9
- ujx = 7.8; ujy = 8.8; ujz = 9.8
- dujx = 3.8; dujy = 4.8; dujz = 5.8
-
elif vertex == 2:
- # Lagrange vertex 12, vertex i: 5, vertex j: 9
- dkp = 1.51709826228; lp = 3.1
- dkq = -0.54615537442; lq = 4.1
- dkr = 0.18205179147; lr = 5.1
+ # Lagrange vertex 12, vertex N: 5, vertex P: 9
+ jL = 1.0/3.0
+ rLx = (7.6-7.2 + 7.8-7.3 + 7.0-7.4)/9.0 + \
+ (-0.07938069066 + -0.14140241667 + -0.18205179147)/9.0
+ rLy = (8.6-8.2 + 8.8-8.3 + 9.0-9.4)/9.0 + \
+ (-1.69682900001 + -1.82575588523 + -1.51709826228)/9.0
+ rLz = (9.6-9.2 + 9.8-9.3 + 9.0-9.4)/9.0 + \
+ (0.55566483464 + 0.56560966667 + 0.54615537442)/9.0
+ jN = 1.4; duNx = 3.4; duNy = 4.4; duNz = 5.4;
+ jP = 1.0; duPx = 3.0; duPy = 4.0; duPz = 5.0;
- # vertex i: 5
- Aix = 1.4; Aiy = 1.4; Aiz = 1.4
- rix = +9.1; riy = +7.1; riz = +8.1
- uix = 7.4; uiy = 8.4; uiz = 9.4
- duix = 3.4; duiy = 4.4; duiz = 5.4
-
- # vertex j: 9
- Ajx = 1.0; Ajy = 1.0; Ajz = 1.0
- rjx = -9.1; rjy = -7.1; rjz = -8.1
- ujx = 7.0; ujy = 8.0; ujz = 9.0
- dujx = 3.0; dujy = 4.0; dujz = 5.0
-
elif test == "hex8":
- Cpx = 0.0
- Cpy = -1.0
- Cpz = 0.0
- Cqx = 0.0
- Cqy = 0.0
- Cqz = +1.0
- Crx = -1.0
- Cry = 0.0
- Crz = 0.0
-
if vertex == 0:
- # Lagrange vertex 18, vertex i: 6, vertex j: 14
- dkp = 1.82575588523; lp = 3.4
- dkq = -0.55566483464; lq = 4.4
- dkr = 0.07938069066; lr = 5.4
+ # Lagrange vertex 18, vertex N: 6, vertex P: 14
+ jL = 1.0
+ rLx = 0.62200847*(0.62200847*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)+0.04465820*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)) + \
+ 0.16666667*(0.16666667*(5.3-4.5+0.07938069066)+0.62200847*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)+0.04465820*(5.5-4.6+0.14140241667)) + \
+ 0.04465820*(0.04465820*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)+0.62200847*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)) + \
+ 0.16666667*(0.16666667*(5.3-4.5+0.07938069066)+0.04465820*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)+0.62200847*(5.5-4.6+0.14140241667))
+ rLy = 0.62200847*(0.62200847*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)+0.04465820*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)) + \
+ 0.16666667*(0.16666667*(7.3-6.5+1.82575588523)+0.62200847*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)+0.04465820*(7.5-6.6+1.69682900001)) + \
+ 0.04465820*(0.04465820*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)+0.62200847*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)) + \
+ 0.16666667*(0.16666667*(7.3-6.5+1.82575588523)+0.04465820*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)+0.62200847*(7.5-6.6+1.69682900001))
+ rLz = 0.62200847*(0.62200847*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)+0.04465820*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)) + \
+ 0.16666667*(0.16666667*(9.3-8.5+0.55566483464)+0.62200847*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)+0.04465820*(9.5-8.6+0.56560966667)) + \
+ 0.04465820*(0.04465820*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)+0.62200847*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)) + \
+ 0.16666667*(0.16666667*(9.3-8.5+0.55566483464)+0.04465820*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)+0.62200847*(9.5-8.6+0.56560966667))
+ jN = 1.5; duNx = 3.5; duNy = 4.5; duNz = 5.5;
+ jP = 1.3; duPx = 3.3; duPy = 4.3; duPz = 5.3;
- # vertex i: 6
- Aix = 1.5; Aiy = 1.5; Aiz = 1.5
- rix = -9.4; riy = -5.4; riz = +7.4
- uix = 4.5; uiy = 6.5; uiz = 8.5
- duix = 3.5; duiy = 4.5; duiz = 5.5
-
- # vertex j: 14
- Ajx = 1.3; Ajy = 1.3; Ajz = 1.3
- rjx = +9.4; rjy = +5.4; rjz = -7.4
- ujx = 5.3; ujy = 7.3; ujz = 9.3
- dujx = 3.3; dujy = 4.3; dujz = 5.3
-
elif vertex == 1:
- # Lagrange vertex 19, vertex i: 7, vertex j: 15
- dkp = 1.69682900001; lp = 3.6
- dkq = -0.56560966667; lq = 4.6
- dkr = 0.14140241667; lr = 5.6
+ # Lagrange vertex 19, vertex N: 7, vertex P: 15
+ jL = 1.0
+ rLx = 0.62200847*(0.62200847*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)+0.04465820*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)) + \
+ 0.16666667*(0.16666667*(5.5-4.6+0.14140241667)+0.62200847*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)+0.04465820*(5.9-4.8+0.19904410828)) + \
+ 0.04465820*(0.04465820*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)+0.62200847*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)) + \
+ 0.16666667*(0.16666667*(5.5-4.6+0.14140241667)+0.04465820*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)+0.62200847*(5.9-4.8+0.19904410828))
+ rLy = 0.62200847*(0.62200847*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)+0.04465820*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)) + \
+ 0.16666667*(0.16666667*(7.5-6.6+1.69682900001)+0.62200847*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)+0.04465820*(7.9-6.8+1.29378670385)) + \
+ 0.04465820*(0.04465820*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)+0.62200847*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)) + \
+ 0.16666667*(0.16666667*(7.5-6.6+1.69682900001)+0.04465820*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)+0.62200847*(7.9-6.8+1.29378670385))
+ rLz = 0.62200847*(0.62200847*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)+0.04465820*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)) + \
+ 0.16666667*(0.16666667*(9.5-8.6+0.56560966667)+0.62200847*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)+0.04465820*(9.9-8.8+0.49761027071)) + \
+ 0.04465820*(0.04465820*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)+0.62200847*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)) + \
+ 0.16666667*(0.16666667*(9.5-8.6+0.56560966667)+0.04465820*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)+0.62200847*(9.9-8.8+0.49761027071))
+ jN = 1.6; duNx = 3.6; duNy = 4.6; duNz = 5.6;
+ jP = 1.5; duPx = 3.5; duPy = 4.5; duPz = 5.5;
- # vertex i: 7
- Aix = 1.6; Aiy = 1.6; Aiz = 1.6
- rix = -9.6; riy = -5.6; riz = +7.6
- uix = 4.6; uiy = 6.6; uiz = 8.6
- duix = 3.6; duiy = 4.6; duiz = 5.6
-
- # vertex j: 15
- Ajx = 1.5; Ajy = 1.5; Ajz = 1.5
- rjx = +9.6; rjy = +5.6; rjz = -7.6
- ujx = 5.5; ujy = 7.5; ujz = 9.5
- dujx = 3.5; dujy = 4.5; dujz = 5.5
-
elif vertex == 2:
- # Lagrange vertex 20, vertex i: 8, vertex j: 16
- dkp = 1.51709826228; lp = 3.8
- dkq = -0.54615537442; lq = 4.8
- dkr = 0.18205179147; lr = 5.8
+ # Lagrange vertex 20, vertex N: 8, vertex P: 16
+ jL = 1.0
+ rLx = 0.62200847*(0.62200847*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)+0.04465820*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)) + \
+ 0.16666667*(0.16666667*(5.7-4.7+0.18205179147)+0.62200847*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)+0.04465820*(5.3-4.5+0.07938069066)) + \
+ 0.04465820*(0.04465820*(5.7-4.7+0.18205179147)+0.16666667*(5.9-4.8+0.19904410828)+0.62200847*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)) + \
+ 0.16666667*(0.16666667*(5.7-4.7+0.18205179147)+0.04465820*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)+0.62200847*(5.3-4.5+0.07938069066))
+ rLy = 0.62200847*(0.62200847*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)+0.04465820*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)) + \
+ 0.16666667*(0.16666667*(7.7-6.7+1.51709826228)+0.62200847*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)+0.04465820*(7.3-6.5+1.82575588523)) + \
+ 0.04465820*(0.04465820*(7.7-6.7+1.51709826228)+0.16666667*(7.9-6.8+1.29378670385)+0.62200847*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)) + \
+ 0.16666667*(0.16666667*(7.7-6.7+1.51709826228)+0.04465820*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)+0.62200847*(7.3-6.5+1.82575588523))
+ rLz = 0.62200847*(0.62200847*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)+0.04465820*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)) + \
+ 0.16666667*(0.16666667*(9.7-8.7+0.54615537442)+0.62200847*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)+0.04465820*(9.3-8.5+0.55566483464)) + \
+ 0.04465820*(0.04465820*(9.7-8.7+0.54615537442)+0.16666667*(9.9-8.8+0.49761027071)+0.62200847*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)) + \
+ 0.16666667*(0.16666667*(9.7-8.7+0.54615537442)+0.04465820*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)+0.62200847*(9.3-8.5+0.55566483464))
+ jN = 1.7; duNx = 3.7; duNy = 4.7; duNz = 5.7;
+ jP = 1.7; duPx = 3.7; duPy = 4.7; duPz = 5.7;
- # vertex i: 8
- Aix = 1.7; Aiy = 1.7; Aiz = 1.7
- rix = -9.8; riy = -5.8; riz = +7.8
- uix = 4.7; uiy = 6.7; uiz = 8.7
- duix = 3.7; duiy = 4.7; duiz = 5.7
-
- # vertex j: 16
- Ajx = 1.7; Ajy = 1.7; Ajz = 1.7
- rjx = +9.8; rjy = +5.8; rjz = -7.8
- ujx = 5.7; ujy = 7.7; ujz = 9.7
- dujx = 3.7; dujy = 4.7; dujz = 5.7
-
elif vertex == 3:
- # Lagrange vertex 21, vertex i: 9, vertex j: 17
- dkp = 1.29378670385; lp = 3.0
- dkq = -0.49761027071; lq = 4.0
- dkr = 0.19904410828; lr = 5.0
+ # Lagrange vertex 21, vertex N: 9, vertex P: 17
+ jL = 1.0
+ rLx = 0.62200847*(0.62200847*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)+0.04465820*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)) + \
+ 0.16666667*(0.16666667*(5.9-4.8+0.19904410828)+0.62200847*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)+0.04465820*(5.7-4.7+0.18205179147)) + \
+ 0.04465820*(0.04465820*(5.9-4.8+0.19904410828)+0.16666667*(5.5-4.6+0.14140241667)+0.62200847*(5.3-4.5+0.07938069066)+0.16666667*(5.7-4.7+0.18205179147)) + \
+ 0.16666667*(0.16666667*(5.9-4.8+0.19904410828)+0.04465820*(5.5-4.6+0.14140241667)+0.16666667*(5.3-4.5+0.07938069066)+0.62200847*(5.7-4.7+0.18205179147))
+ rLy = 0.62200847*(0.62200847*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)+0.04465820*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)) + \
+ 0.16666667*(0.16666667*(7.9-6.8+1.29378670385)+0.62200847*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)+0.04465820*(7.7-6.7+1.51709826228)) + \
+ 0.04465820*(0.04465820*(7.9-6.8+1.29378670385)+0.16666667*(7.5-6.6+1.69682900001)+0.62200847*(7.3-6.5+1.82575588523)+0.16666667*(7.7-6.7+1.51709826228)) + \
+ 0.16666667*(0.16666667*(7.9-6.8+1.29378670385)+0.04465820*(7.5-6.6+1.69682900001)+0.16666667*(7.3-6.5+1.82575588523)+0.62200847*(7.7-6.7+1.51709826228))
+ rLz = 0.62200847*(0.62200847*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)+0.04465820*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)) + \
+ 0.16666667*(0.16666667*(9.9-8.8+0.49761027071)+0.62200847*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)+0.04465820*(9.7-8.7+0.54615537442)) + \
+ 0.04465820*(0.04465820*(9.9-8.8+0.49761027071)+0.16666667*(9.5-8.6+0.56560966667)+0.62200847*(9.3-8.5+0.55566483464)+0.16666667*(9.7-8.7+0.54615537442)) + \
+ 0.16666667*(0.16666667*(9.9-8.8+0.49761027071)+0.04465820*(9.5-8.6+0.56560966667)+0.16666667*(9.3-8.5+0.55566483464)+0.62200847*(9.7-8.7+0.54615537442))
+ jN = 1.8; duNx = 3.8; duNy = 4.8; duNz = 5.8;
+ jP = 1.9; duPx = 3.9; duPy = 4.9; duPz = 5.9;
- # vertex i: 9
- Aix = 1.8; Aiy = 1.8; Aiz = 1.8
- rix = -9.0; riy = -5.0; riz = +7.0
- uix = 4.8; uiy = 6.8; uiz = 8.8
- duix = 3.8; duiy = 4.8; duiz = 5.8
-
- # vertex j: 17
- Ajx = 1.9; Ajy = 1.9; Ajz = 1.9
- rjx = +9.0; rjy = +5.0; rjz = -7.0
- ujx = 5.9; ujy = 7.9; ujz = 9.9
- dujx = 3.9; dujy = 4.9; dujz = 5.9
-
- Sppi = Aix*Ajx / (Aix + Ajx)
- Sqqi = Aiy*Ajy / (Aiy + Ajy)
- Srri = Aiz*Ajz / (Aiz + Ajz)
+ Sinv = (jN + jP) / (jL*jL)
+ duLx = Sinv * (-rLx + jL*(duPx - duNx))
+ duLy = Sinv * (-rLy + jL*(duPy - duNy))
+ duLz = Sinv * (-rLz + jL*(duPz - duNz))
- Arux = rix / Aix - rjx / Ajx + uix - ujx
- Aruy = riy / Aiy - rjy / Ajy + uiy - ujy
- Aruz = riz / Aiz - rjz / Ajz + uiz - ujz
+ dduNx = jL / jN * duLx
+ dduNy = jL / jN * duLy
+ dduNz = jL / jN * duLz
- Arup = Cpx*Arux + Cpy*Aruy + Cpz*Aruz
- Aruq = Cqx*Arux + Cqy*Aruy + Cqz*Aruz
- Arur = Crx*Arux + Cry*Aruy + Crz*Aruz
- Arupslip = -Arup - dkp
- Aruqslip = -Aruq - dkq
- Arurslip = -Arur - dkr
+ dduPx = -jL / jP * duLx
+ dduPy = -jL / jP * duLy
+ dduPz = -jL / jP * duLz
+
+ print duNx+dduNx,duNy+dduNy, duNz+dduNz
+ print duPx+dduPx,duPy+dduPy, duPz+dduPz
+ print duLx,duLy, duLz
- dlp = Sppi * Arupslip
- dlq = Sqqi * Aruqslip
- dlr = Srri * Arurslip
-
- dduix = +1.0/Aix * (Cpx*dlp + Cqx*dlq + Crx*dlr)
- dduiy = +1.0/Aiy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
- dduiz = +1.0/Aiy * (Cpz*dlp + Cqz*dlq + Crz*dlr)
-
- ddujx = -1.0/Ajx * (Cpx*dlp + Cqx*dlq + Crx*dlr)
- ddujy = -1.0/Ajy * (Cpy*dlp + Cqy*dlq + Cry*dlr)
- ddujz = -1.0/Ajz * (Cpz*dlp + Cqz*dlq + Crz*dlr)
-
-# print "Sppi",Sppi
-# print "Sqqi",Sqqi
-# print "Srri",Srri
-
-# print "Arup",Arup
-# print "Aruq",Aruq
-# print "Arur",Arur
-
-# print "Arupslip",Arupslip
-# print "Aruqslip",Aruqslip
-# print "Arurslip",Arurslip
-
-# print "dlp",dlp
-# print "dlq",dlq
-# print "dlr",dlr
-
-# print "dduix:",dduix
-# print "dduiy:",dduiy
-# print "dduiz:",dduiz
-
-# print "ddujx:",ddujx
-# print "ddujy:",ddujy
-# print "ddujz:",ddujz
-
- print duix+dduix,duiy+dduiy,duiz+dduiz
- print dujx+ddujx,dujy+ddujy,dujz+ddujz
- print dlp,dlq,dlr
-
More information about the CIG-COMMITS
mailing list