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

brad at geodynamics.org brad at geodynamics.org
Mon Jan 1 13:23:51 PST 2007


Author: brad
Date: 2007-01-01 13:23:51 -0800 (Mon, 01 Jan 2007)
New Revision: 5630

Modified:
   short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
   short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
   short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
Log:
Initial implementation of using spatial database for density. Not tested.

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc	2007-01-01 21:05:23 UTC (rev 5629)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.cc	2007-01-01 21:23:51 UTC (rev 5630)
@@ -12,10 +12,11 @@
 
 #include <portinfo>
 
+#include "Integrator.hh" // implementation of class methods
+
 #include "Quadrature.hh"
+#include "spatialdata/spatialdb/SpatialDB.hh"
 
-#include "Integrator.hh" // implementation of class methods
-
 // ----------------------------------------------------------------------
 // Constructor
 pylith::feassemble::Integrator::Integrator(void) :

Modified: short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-01-01 21:05:23 UTC (rev 5629)
+++ short/3D/PyLith/trunk/libsrc/feassemble/Integrator.hh	2007-01-01 21:23:51 UTC (rev 5630)
@@ -31,6 +31,15 @@
   } // feassemble
 } // pylith
 
+namespace spatialdata {
+  namespace spatialdb {
+    class SpatialDB; // USES SpatialDB
+  } // spatialdb
+  namespace geocoords {
+    class CoordSys; // USES CoordSys
+  } // geocoords
+} // spatialdata
+
 class pylith::feassemble::Integrator
 { // Integrator
   friend class TestIntegrator; // unit testing
@@ -84,6 +93,22 @@
    */
   void quadrature(const Quadrature* q);
 
+  /** Set spatial database for material properties.
+   *
+   * @param db Pointer to spatial database
+   */
+  void database(const spatialdata::spatialdb::SpatialDB* db);
+
+  /** Initialize, get material property parameters from database.
+   *
+   * @param mesh PETSc mesh
+   * @param db Pointer to spatial database with material property parameters
+   * @param cs Pointer to coordinate system of vertices
+   */
+  virtual void initialize(ALE::Obj<ALE::Mesh>& mesh,
+			  spatialdata::spatialdb::SpatialDB* db,
+			  spatialdata::geocoords::CoordSys* cs) = 0;
+
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh	2007-01-01 21:05:23 UTC (rev 5629)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity3D.hh	2007-01-01 21:23:51 UTC (rev 5630)
@@ -63,6 +63,16 @@
 		 const ALE::Obj<ALE::Mesh::real_section_type>& fieldIn,
 		 const ALE::Obj<ALE::Mesh::real_section_type>& coordinates);
 
+  /** Initialize, get material property parameters from database.
+   *
+   * @param mesh PETSc mesh
+   * @param db Pointer to spatial database with material property parameters
+   * @param cs Pointer to coordinate system of vertices
+   */
+  void initialize(ALE::Obj<ALE::Mesh>& mesh,
+		  spatialdata::spatialdb::SpatialDB* db,
+		  spatialdata::geocoords::CoordSys* cs);
+
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc	2007-01-01 21:05:23 UTC (rev 5629)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.cc	2007-01-01 21:23:51 UTC (rev 5630)
@@ -17,6 +17,7 @@
 #include "Quadrature.hh" // USES Quadrature
 
 #include "petscmat.h" // USES PetscMat
+#include "spatialdata/spatialdb/SpatialDB.hh"
 
 #include <assert.h> // USES assert()
 #include <stdexcept> // USES std::runtime_error
@@ -84,8 +85,8 @@
     const int spaceDim = _quadrature->spaceDim();
 
     // Restrict density from material database to quadrature points for this cell
-    const real_section_type::value_type* _density->restrict(patch, *cellIter);
-    // ALSO RUN tractest IN PARALLEL (PROBLEM WITH SURFACE MESH DISTRIBUTION) pylith-0.8/pylith3d/examples/lintet/tractest
+    const real_section_type::value_type* density = 
+      _density->restrict(patch, *cellIter);
 
     // Compute action for cell
     for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -260,5 +261,70 @@
   } // for
 } // integrateLumped
 
+// ----------------------------------------------------------------------
+// Initialize, get material property parameters from database.
+void
+pylith::feassemble::IntegratorInertia::initialize(
+				     ALE::Obj<ALE::Mesh>& mesh,
+				     spatialdata::spatialdb::SpatialDB* db,
+				     spatialdata::geocoords::CoordSys* cs)
+{ // initialize
+  typedef ALE::Mesh::real_section_type real_section_type;
+  typedef ALE::Mesh::topology_type topology_type;
 
+  assert(0 != db);
+  assert(0 != cs);
+
+  // Create density section
+  const ALE::Mesh::int_section_type::patch_type patch = 0;
+  _density = mesh->getRealSection("density");
+  const int fiberDim = 1; // number of values in field per cell
+  _density->setName("fieldOut");
+  _density->setFiberDimensionByDepth(patch, 0, fiberDim);
+  _density->allocate();
+
+  // Open database
+  db->open();
+  const int numVals = 1;
+  const char* names[numVals];
+  names[0] = "density";
+  db->queryVals(names, numVals);
+  
+  const ALE::Obj<real_section_type>& coordinates = 
+    mesh->getRealSection("coordinates");
+  const ALE::Obj<topology_type>& topology = coordinates->getTopology();
+  const ALE::Obj<topology_type::label_sequence>& cells = 
+    topology->heightStratum(patch, 0);
+  const topology_type::label_sequence::iterator cellsEnd = cells->end();
+
+  // Loop over cells
+  const int numQuadPts = _quadrature->numQuadPts();
+  double* cellDensity = (numQuadPts > 0) ? new double[numQuadPts] : 0;
+  for (topology_type::label_sequence::iterator cellIter=cells->begin();
+       cellIter != cellsEnd;
+       ++cellIter) {
+    // Compute geometry information for current cell
+    _quadrature->computeGeometry(coordinates, *cellIter);
+
+    const double* quadPts = _quadrature->quadPts();
+    const int spaceDim = _quadrature->spaceDim();
+
+    // Loop over quadrature points in cell and query database
+    for (int iQuadPt=0, index=0; 
+	 iQuadPt < numQuadPts; 
+	 ++iQuadPt, index+=spaceDim)
+      const int err = db->query(&cellDensity[iQuadPt], numVals, 
+				quadPts[index], 
+				quadPts[index+1], 
+				quadPts[index+2], cs);
+    // Assemble cell contribution into field
+    _density->updateAdd(patch, *cellIter, cellDensity);
+  } // for
+  delete[] cellDensity; cellDensity = 0;
+
+  // Close database
+  db->close();
+} // initialize
+
+
 // End of file 

Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh	2007-01-01 21:05:23 UTC (rev 5629)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorInertia.hh	2007-01-01 21:23:51 UTC (rev 5630)
@@ -86,9 +86,16 @@
   void integrateLumped(const ALE::Obj<real_section_type>& fieldOut,
 		       const ALE::Obj<real_section_type>& coordinates);
 
+  /** Initialize, get material property parameters from database.
+   *
+   * @param mesh PETSc mesh
+   * @param db Pointer to spatial database with material property parameters
+   * @param cs Pointer to coordinate system of vertices
+   */
+  void initialize(ALE::Obj<ALE::Mesh>& mesh,
+		  spatialdata::spatialdb::SpatialDB* db,
+		  spatialdata::geocoords::CoordSys* cs);
 
-  void setDensity(const ALE::Obj<real_section_type>& density) {_density = density};
-
 // PROTECTED METHODS ////////////////////////////////////////////////////
 protected :
 
@@ -98,14 +105,19 @@
    */
   IntegratorInertia(const IntegratorInertia& i);
 
-  ALE::Obj<real_section_type> _density;
-
 // PRIVATE METHODS //////////////////////////////////////////////////////
 private :
 
   /// Not implemented
   const IntegratorInertia& operator=(const IntegratorInertia&);
 
+// PRIVATE MEMBERS //////////////////////////////////////////////////////
+private :
+
+  /// Field with density at quadrature points
+  /// Fiber dimension is number of quadrature points per cell
+  ALE::Obj<real_section_type> _density;
+
 }; // IntegratorInertia
 
 #include "IntegratorInertia.icc" // inline methods



More information about the cig-commits mailing list