[cig-commits] r4556 - short/3D/PyLith/trunk/playpen/integrate/src

[knepley at geodynamics.org]

Author: knepley
Date: 2006-09-15 11:17:25 -0700 (Fri, 15 Sep 2006)
New Revision: 4556

Modified:
   short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
   short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
   short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
Log:
Fixed integrator test


Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc	2006-09-15 17:41:14 UTC (rev 4555)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.cc	2006-09-15 18:17:25 UTC (rev 4556)
@@ -45,8 +45,21 @@
   }
 };
 
-void pylith::feassemble::Integrator::integrateFunction(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(value_type []))
+void pylith::feassemble::Integrator::fillSection(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const value_type []))
 {
+  const topology_type::patch_type patch = 0;
+  const section_type::chart_type& chart = field->getPatch(patch);
+
+  for(section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
+    const section_type::value_type *coords = coordinates->restrict(patch, *c_iter);
+    const section_type::value_type  value  = f(coords);
+
+    field->updatePoint(patch, *c_iter, &value);
+  }
+};
+
+void pylith::feassemble::Integrator::integrateFunction(const Obj<section_type>& field, const Obj<section_type>& coordinates, value_type (*f)(const value_type []))
+{
   const topology_type::patch_type               patch    = 0;
   const Obj<topology_type>&                     topology = field->getTopology();
   const Obj<topology_type::label_sequence>&     elements = topology->heightStratum(patch, 0);
@@ -56,6 +69,14 @@
   for(topology_type::label_sequence::iterator e_iter = elements->begin(); e_iter != end; ++e_iter) {
     computeElementGeometry(coordinates, *e_iter, _v0, _Jac, PETSC_NULL, detJ);
     // Element integral
+    for(int d = 0; d < _dim; d++) {
+      std::cout << "v0["<<d<<"]: "<<_v0[d]<<std::endl;
+    }
+    for(int d = 0; d < _dim; d++) {
+      for(int e = 0; e < _dim; e++) {
+        std::cout << "Jac["<<d<<", "<<e<<"]: " << _Jac[d*_dim+e] << std::endl;
+      }
+    }
     PetscMemzero(_elemVector, NUM_BASIS_FUNCTIONS*sizeof(value_type));
     for(int q = 0; q < NUM_QUADRATURE_POINTS; q++) {
       value_type xi  = _points[q*2+0] + 1.0;

Modified: short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh	2006-09-15 17:41:14 UTC (rev 4555)
+++ short/3D/PyLith/trunk/playpen/integrate/src/elemVector.hh	2006-09-15 18:17:25 UTC (rev 4556)
@@ -70,7 +70,8 @@
       };
     public:
       void computeElementGeometry(const Obj<section_type>&, const point_type&, value_type [], value_type [], value_type [], value_type&);
-      void integrateFunction(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(value_type []));
+      void fillSection(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
+      void integrateFunction(const Obj<section_type>&, const Obj<section_type>&, value_type (*)(const value_type []));
       void integrateLaplacianAction(const Obj<section_type>&, const Obj<section_type>&, const Obj<section_type>&);
     };
   }

Modified: short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc
===================================================================
--- short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc	2006-09-15 17:41:14 UTC (rev 4555)
+++ short/3D/PyLith/trunk/playpen/integrate/src/testintegrate.cc	2006-09-15 18:17:25 UTC (rev 4556)
@@ -17,7 +17,7 @@
 #include "Integration.hh"
 
 // ----------------------------------------------------------------------
-ALE::Mesh::section_type::value_type foo(ALE::Mesh::section_type::value_type coords[])
+ALE::Mesh::section_type::value_type foo(const ALE::Mesh::section_type::value_type coords[])
 {
   return 1.0;
 }
@@ -43,19 +43,22 @@
     iohandler.read(mesh, false);
 
     const ALE::Mesh::topology_type::patch_type patch = 0;
+    const Obj<ALE::Mesh::section_type>& coords = mesh->getSection("coordinates");
     const Obj<ALE::Mesh::section_type>& X = mesh->getSection("X");
     const Obj<ALE::Mesh::section_type>& F = mesh->getSection("F");
     pylith::feassemble::Integrator      integrator(mesh->getDimension(),
                                                    NUM_QUADRATURE_POINTS, points, weights,
                                                    NUM_BASIS_FUNCTIONS, Basis, BasisDerivatives);
 
+    coords->view("Coordinates");
     X->setFiberDimensionByDepth(patch, 0, 1);
     X->allocate();
     F->setFiberDimensionByDepth(patch, 0, 1);
     F->allocate();
-    integrator.integrateFunction(X, mesh->getSection("coordinates"), foo);
+    integrator.integrateFunction(X, coords, foo);
     X->view("Weak form of foo");
-    integrator.integrateLaplacianAction(X, F, mesh->getSection("coordinates"));
+    integrator.fillSection(X, coords, foo);
+    integrator.integrateLaplacianAction(X, F, coords);
     F->view("Weak form of \Delta foo");
   } catch(ALE::Exception e) {
     int rank;
@@ -70,4 +73,4 @@
 // version
 // $Id$
 
-// End of file 
+// End of file