[cig-commits] r16623 - short/3D/PyLith/trunk/libsrc/problems

[brad at geodynamics.org]

Author: brad
Date: 2010-05-04 12:58:46 -0700 (Tue, 04 May 2010)
New Revision: 16623

Modified:
   short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
   short/3D/PyLith/trunk/libsrc/problems/SolverLumped.cc
Log:
Use optimized restrict() with no array as arg.

Modified: short/3D/PyLith/trunk/libsrc/problems/Explicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/Explicit.cc	2010-05-04 18:51:01 UTC (rev 16622)
+++ short/3D/PyLith/trunk/libsrc/problems/Explicit.cc	2010-05-04 19:58:46 UTC (rev 16623)
@@ -56,15 +56,12 @@
   const int spaceDim = cs->spaceDim();
   
   // Get sections.
-  double_array dispIncrVertex(spaceDim);
   const ALE::Obj<RealSection>& dispIncrSection = dispIncr.section();
   assert(!dispIncrSection.isNull());
 	 
-  double_array dispTVertex(spaceDim);
   const ALE::Obj<RealSection>& dispTSection = _fields->get("disp(t)").section();
   assert(!dispTSection.isNull());
 
-  double_array dispTmdtVertex(spaceDim);
   const ALE::Obj<RealSection>& dispTmdtSection =
     _fields->get("disp(t-dt)").section();
   assert(!dispTmdtSection.isNull());
@@ -91,15 +88,21 @@
   for (SieveMesh::label_sequence::iterator v_iter=verticesBegin; 
        v_iter != verticesEnd;
        ++v_iter) {
-    dispIncrSection->restrictPoint(*v_iter, &dispIncrVertex[0],
-				   dispIncrVertex.size());
-    dispTSection->restrictPoint(*v_iter, &dispTVertex[0],
-				dispTVertex.size());
-    dispTmdtSection->restrictPoint(*v_iter, &dispTmdtVertex[0],
-				   dispTmdtVertex.size());
+    assert(spaceDim == dispIncrSection->getFiberDimension(*v_iter));
+    const double* dispIncrVertex = dispIncrSection->restrictPoint(*v_iter);
 
-    velVertex = (dispIncrVertex + dispTVertex - dispTmdtVertex) / twodt;
-    accVertex = (dispIncrVertex - dispTVertex + dispTmdtVertex) / dt2;
+    assert(spaceDim == dispTSection->getFiberDimension(*v_iter));
+    const double* dispTVertex = dispTSection->restrictPoint(*v_iter);
+
+    assert(spaceDim == dispTmdtSection->getFiberDimension(*v_iter));
+    const double* dispTmdtVertex = dispTmdtSection->restrictPoint(*v_iter);
+
+    for (int i=0; i < spaceDim; ++i) {
+      velVertex[i] = 
+	(dispIncrVertex[i] + dispTVertex[i] - dispTmdtVertex[i]) / twodt;
+      accVertex[i] = 
+	(dispIncrVertex[i] - dispTVertex[i] + dispTmdtVertex[i]) / dt2;
+    } // for
     
     assert(velSection->getFiberDimension(*v_iter) == spaceDim);
     velSection->updatePoint(*v_iter, &velVertex[0]);

Modified: short/3D/PyLith/trunk/libsrc/problems/SolverLumped.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/problems/SolverLumped.cc	2010-05-04 18:51:01 UTC (rev 16622)
+++ short/3D/PyLith/trunk/libsrc/problems/SolverLumped.cc	2010-05-04 19:58:46 UTC (rev 16623)
@@ -97,11 +97,9 @@
   const ALE::Obj<RealSection>& solutionSection = solution->section();
   assert(!solutionSection.isNull());
 	 
-  double_array jacobianVertex(spaceDim);
   const ALE::Obj<RealSection>& jacobianSection = jacobian.section();
   assert(!jacobianSection.isNull());
 
-  double_array residualVertex(spaceDim);
   const ALE::Obj<RealSection>& residualSection = residual.section();
   assert(!residualSection.isNull());
   
@@ -111,17 +109,16 @@
   for (SieveMesh::label_sequence::iterator v_iter=verticesBegin; 
        v_iter != verticesEnd;
        ++v_iter) {
-    jacobianSection->restrictPoint(*v_iter, &jacobianVertex[0],
-				   jacobianVertex.size());
-    residualSection->restrictPoint(*v_iter, &residualVertex[0],
-				   residualVertex.size());
+    assert(spaceDim == jacobianSection->getFiberDimension(*v_iter));
+    const double* jacobianVertex = jacobianSection->restrictPoint(*v_iter);
 
-#if !defined(NDEBUG)
-    for (int i=0; i < spaceDim; ++i)
-      assert(jacobianVertex[i] != 0.0);
-#endif
+    assert(spaceDim == residualSection->getFiberDimension(*v_iter));
+    const double* residualVertex = residualSection->restrictPoint(*v_iter);
 
-    solutionVertex = residualVertex / jacobianVertex;
+    for (int i=0; i < spaceDim; ++i) {
+      assert(jacobianVertex[i] != 0.0);
+      solutionVertex[i] = residualVertex[i] / jacobianVertex[i];
+    } // for
     
     assert(solutionSection->getFiberDimension(*v_iter) == spaceDim);
     solutionSection->updatePoint(*v_iter, &solutionVertex[0]);