[cig-commits] r7775 - short/3D/PyLith/trunk/libsrc/bc

[willic3 at geodynamics.org]

Author: willic3
Date: 2007-08-03 13:45:12 -0700 (Fri, 03 Aug 2007)
New Revision: 7775

Modified:
   short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
   short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
Log:
Fixed a few problems.  The main problem fixed was confusing the
dimension of the reference cell quadrature (cellDim) with the dimension
of the global geometry (spaceDim).



Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.cc	2007-08-03 19:01:02 UTC (rev 7774)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.cc	2007-08-03 20:45:12 UTC (rev 7775)
@@ -16,9 +16,12 @@
 
 #include "pylith/feassemble/Quadrature.hh" // USES Quadrature
 #include "pylith/feassemble/CellGeometry.hh" // USES CellGeometry
+#include "pylith/topology/FieldsManager.hh" // USES FieldsManager
+#include "pylith/utils/array.hh" // USES double_array
+#include <petscmat.h> // USES PETSc Mat
 
+#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
 #include "spatialdata/spatialdb/SpatialDB.hh" // USES SpatialDB
-#include "spatialdata/geocoords/CoordSys.hh" // USES CoordSys
 
 #include <Selection.hh> // USES submesh algorithms
 
@@ -157,7 +160,8 @@
   // reference geometry.
   double_array tractionDataLocal(spaceDim);
   const double_array& quadPtsRef = _quadrature->quadPtsRef();
-  double_array quadPtRef(spaceDim);
+  // double_array quadPtRef(spaceDim);
+  const double_array& quadPts = _quadrature->quadPts();
 
   // Container for cell tractions rotated to global coordinates.
   double_array cellTractionsGlobal(fiberDim);
@@ -170,18 +174,17 @@
       ++c_iter) {
     _quadrature->computeGeometry(_boundaryMesh, coordinates, *c_iter);
     mesh->restrict(coordinates, *c_iter, &cellVertices[0], cellVertices.size());
-    const double_array& quadPts = _quadrature->quadPts();
 
-    for(int iQuad = 0, index=0; iQuad < numQuadPts; ++iQuad, index+=spaceDim) {
+    for(int iQuad = 0, iRef=0, iSpace=0; iQuad < numQuadPts; ++iQuad, iRef+=cellDim, iSpace+=spaceDim) {
       // Get traction vector in local coordinate system at quadrature point
       const int err = _db->query(&tractionDataLocal[0], spaceDim,
-				 &quadPts[index], spaceDim, cs);
+				 &quadPts[iSpace], spaceDim, cs);
       if (err) {
 	std::ostringstream msg;
 	msg << "Could not find traction values at \n"
 	    << "(";
 	for (int i=0; i < spaceDim; ++i)
-	  msg << " " << quadPts[index+spaceDim];
+	  msg << " " << quadPts[i+iSpace];
 	msg << ") for traction boundary condition " << _label << "\n"
 	    << "using spatial database " << _db->label() << ".";
 	throw std::runtime_error(msg.str());
@@ -189,7 +192,7 @@
 
       // Compute Jacobian and determinant at quadrature point, then get
       // orientation.
-      memcpy(&quadPtRef[0], &quadPtsRef[index], spaceDim*sizeof(double));
+      double_array quadPtRef(&quadPtsRef[iRef], cellDim);
       cellGeometry.jacobian(&jacobian, &jacobianDet, cellVertices, quadPtRef);
       cellGeometry.orientation(&orientation, jacobian, jacobianDet, upDir);
 
@@ -198,7 +201,7 @@
       cellTractionsGlobal = 0.0;
       for(int iDim = 0; iDim < spaceDim; ++iDim) {
 	for(int jDim = 0; jDim < spaceDim; ++jDim)
-	  cellTractionsGlobal[iDim+index] +=
+	  cellTractionsGlobal[iDim+iSpace] +=
 	    orientation[iDim*spaceDim+jDim] * tractionDataLocal[jDim];
       } // for
     } // for

Modified: short/3D/PyLith/trunk/libsrc/bc/Neumann.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/bc/Neumann.hh	2007-08-03 19:01:02 UTC (rev 7774)
+++ short/3D/PyLith/trunk/libsrc/bc/Neumann.hh	2007-08-03 20:45:12 UTC (rev 7775)
@@ -22,8 +22,8 @@
 #include "BoundaryCondition.hh" // ISA BoundaryCondition
 #include "pylith/feassemble/Integrator.hh" // ISA Integrator
 
-#include "pylith/utils/array.hh" // USES std::vector, double_array, int_array
-#include "pylith/utils/sievetypes.hh" // USES real_section_type
+// #include "pylith/utils/array.hh" // USES std::vector, double_array, int_array
+// #include "pylith/utils/sievetypes.hh" // USES real_section_type
 
 /// Namespace for pylith package
 namespace pylith {