[cig-commits] r22564 - short/3D/PyLith/trunk/libsrc/pylith/feassemble

[knepley at geodynamics.org]

Author: knepley
Date: 2013-07-09 14:46:55 -0700 (Tue, 09 Jul 2013)
New Revision: 22564

Modified:
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
   short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
Log:
Converted explicit elasticity residual integration to use presistent arrays

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2013-07-09 18:51:21 UTC (rev 22563)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicit.cc	2013-07-09 21:46:55 UTC (rev 22564)
@@ -237,6 +237,11 @@
   _logger->eventBegin(computeEvent);
 #endif
 
+  PetscScalar *coordsCell = new PetscScalar[64];
+  PetscScalar *accCell    = new PetscScalar[64];
+  PetscScalar *velCell    = new PetscScalar[64];
+  PetscScalar *dispCell   = new PetscScalar[64];
+
   // Loop over cells
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
@@ -244,11 +249,10 @@
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventBegin(geometryEvent);
 #endif
-    PetscScalar *coordsCell = NULL;
-    PetscInt coordsSize = 0;
-    coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
+    PetscInt coordsSize = 64;
+    coordsVisitor.getClosure((PetscScalar **) &coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     _quadrature->computeGeometry(coordsCell, coordsSize, cell);
-    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+    //coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(geometryEvent);
@@ -267,17 +271,14 @@
     _resetCellVector();
 
     // Restrict input fields to cell
-    PetscScalar* accCell = NULL;
-    PetscInt accSize = 0;
-    accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
+    PetscInt accSize = 64;
+    accVisitor.getClosure((PetscScalar **) &accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
 
-    PetscScalar* velCell = NULL;
-    PetscInt velSize = 0;
-    velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
+    PetscInt velSize = 64;
+    velVisitor.getClosure((PetscScalar **) &velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
 
-    PetscScalar* dispCell = NULL;
-    PetscInt dispSize = 0;
-    dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
+    PetscInt dispSize = 64;
+    dispVisitor.getClosure((PetscScalar **) &dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
 #if defined(DETAILED_EVENT_LOGGING)
     _logger->eventEnd(restrictEvent);
@@ -306,7 +307,7 @@
         const int err = db->query(&gravVec[0], gravVec.size(), &quadPtsGlobal[0], spaceDim, cs);
         if (err) {
           throw std::runtime_error("Unable to get gravity vector for point.");
-	} // if
+        } // if
         _normalizer->nondimensionalize(&gravVec[0], gravVec.size(), gravityScale);
         const PylithScalar wt = quadWts[iQuad] * jacobianDet[iQuad] * density[iQuad];
         for (int iBasis = 0, iQ = iQuad * numBasis; iBasis < numBasis; ++iBasis) {
@@ -345,9 +346,9 @@
     for(PetscInt i = 0; i < dispSize; ++i) {
       dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
     } // for
-    accVisitor.restoreClosure(&accCell, &accSize, cell);
-    velVisitor.restoreClosure(&velCell, &velSize, cell);
-    dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+    //accVisitor.restoreClosure(&accCell, &accSize, cell);
+    //velVisitor.restoreClosure(&velCell, &velSize, cell);
+    //dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
 
 #if defined(DETAILED_EVENT_LOGGING)
     PetscLogFlops(numQuadPts*(4+numBasis*3));
@@ -379,6 +380,11 @@
 #endif
   } // for
 
+  delete [] coordsCell;
+  delete [] accCell;
+  delete [] velCell;
+  delete [] dispCell;
+
 #if !defined(DETAILED_EVENT_LOGGING)
   PetscLogFlops(numCells*numQuadPts*(4+numBasis*3));
   _logger->eventEnd(computeEvent);

Modified: short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2013-07-09 18:51:21 UTC (rev 22563)
+++ short/3D/PyLith/trunk/libsrc/pylith/feassemble/ElasticityExplicitTri3.cc	2013-07-09 21:46:55 UTC (rev 22564)
@@ -222,6 +222,11 @@
   _logger->eventBegin(computeEvent);
 #endif
 
+  PetscScalar *coordsCell = new PetscScalar[64];
+  PetscScalar *accCell    = new PetscScalar[64];
+  PetscScalar *velCell    = new PetscScalar[64];
+  PetscScalar *dispCell   = new PetscScalar[64];
+
   // Loop over cells
   for(PetscInt c = 0; c < numCells; ++c) {
     const PetscInt cell = cells[c];
@@ -230,16 +235,13 @@
 #endif
 
     // Restrict input fields to cell
-    PetscScalar* accCell = NULL;
-    PetscInt accSize = 0;
+    PetscInt accSize = 64;
     accVisitor.getClosure(&accCell, &accSize, cell);assert(accCell);assert(numBasis*spaceDim == accSize);
 
-    PetscScalar* velCell = NULL;
-    PetscInt velSize = 0;
+    PetscInt velSize = 64;
     velVisitor.getClosure(&velCell, &velSize, cell);assert(velCell);assert(numBasis*spaceDim == velSize);
 
-    PetscScalar* dispCell = NULL;
-    PetscInt dispSize = 0;
+    PetscInt dispSize = 64;
     dispVisitor.getClosure(&dispCell, &dispSize, cell);assert(dispCell);assert(numBasis*spaceDim == dispSize);
 
 #if defined(DETAILED_EVENT_LOGGING)
@@ -248,8 +250,7 @@
 #endif
 
     // Compute geometry information for current cell
-    PetscScalar *coordsCell = NULL;
-    PetscInt coordsSize = 0;
+    PetscInt coordsSize = 64;
     coordsVisitor.getClosure(&coordsCell, &coordsSize, cell);assert(coordsCell);assert(numBasis*spaceDim == coordsSize);
     const PylithScalar area = _area(coordsCell, coordsSize);assert(area > 0.0);
 
@@ -318,9 +319,9 @@
     for(PetscInt i = 0; i < cellVectorSize; ++i) {
       dispAdjCell[i] = dispCell[i] + viscosity * velCell[i];
     } // for
-    accVisitor.restoreClosure(&accCell, &accSize, cell);
-    velVisitor.restoreClosure(&velCell, &velSize, cell);
-    dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
+    //accVisitor.restoreClosure(&accCell, &accSize, cell);
+    //velVisitor.restoreClosure(&velCell, &velSize, cell);
+    //dispVisitor.restoreClosure(&dispCell, &dispSize, cell);
 
     // Compute B(transpose) * sigma, first computing strains
     const PylithScalar x0 = coordsCell[0];
@@ -373,7 +374,7 @@
     _logger->eventBegin(updateEvent);
 #endif
 
-    coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
+    //coordsVisitor.restoreClosure(&coordsCell, &coordsSize, cell);
 
     // Assemble cell contribution into field
     residualVisitor.setClosure(&_cellVector[0], _cellVector.size(), cell, ADD_VALUES);
@@ -383,6 +384,11 @@
 #endif
   } // for
 
+  delete [] coordsCell;
+  delete [] accCell;
+  delete [] velCell;
+  delete [] dispCell;
+
 #if !defined(DETAILED_EVENT_LOGGING)
   PetscLogFlops(numCells*(2 + numBasis*spaceDim*2 + 34+30));
   _logger->eventEnd(computeEvent);