[cig-commits] r3829 - in short/3D/PyLith/branches/pylith-0.8/pylith3d: libpylith3d module

[knepley at geodynamics.org]

Author: knepley
Date: 2006-06-20 18:42:06 -0700 (Tue, 20 Jun 2006)
New Revision: 3829

Modified:
   short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F
   short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/residu.F
   short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc
Log:
Fixed assembly of full_displacements from increments
Fixed residu() to not overwrite the solution


Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F	2006-06-21 00:34:54 UTC (rev 3828)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/iterate.F	2006-06-21 01:42:06 UTC (rev 3829)
@@ -338,7 +338,7 @@
 c
 c...compute the out-of-balance forces and convergence criteria
 c
-        call residu(rhs,sol,bextern,bintern,bresid,dispvec,
+        call residu(rhs,bextern,bintern,bresid,dispvec,
      &   gtol,gi,gprev,gcurr,
      &   id,idx,neq,nextflag,numnp,iter,itmaxp,idebug,idout,kto,kw,
      &   converge)

Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/residu.F
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/residu.F	2006-06-21 00:34:54 UTC (rev 3828)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/libpylith3d/residu.F	2006-06-21 01:42:06 UTC (rev 3829)
@@ -28,7 +28,7 @@
 c~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 c
 c
-      subroutine residu(rhs,sol,bextern,bintern,bresid,dispvec,gtol,gi,
+      subroutine residu(rhs,bextern,bintern,bresid,dispvec,gtol,gi,
      & gprev,
      & gcurr,id,idx,neq,nextflag,numnp,iter,itmaxp,idebug,idout,kto,kw,
      & converge)
@@ -44,7 +44,7 @@
 c
 c     PETSc objects
 c
-      Vec rhs, sol
+      Vec rhs
       Vec locRhs, locSol
       PetscErrorCode ierr
 c
@@ -116,9 +116,9 @@
         call dcopy(neq,bextern,ione,bresid,ione)
       end if
       call daxpy(neq,alpha,bintern,ione,bresid,ione)
-      call VecZeroEntries(sol, ierr)
-      call assembleVectorComplete(sol, locSol, ADD_VALUES, ierr)
-      call VecNorm(sol, NORM_2, tmp(1), ierr)
+      call VecZeroEntries(rhs, ierr)
+      call assembleVectorComplete(rhs, locSol, ADD_VALUES, ierr)
+      call VecNorm(rhs, NORM_2, tmp(1), ierr)
 cpetsc tmp(1)=dnrm2(neq,dispvec,ione)
       call VecZeroEntries(rhs, ierr)
       call assembleVectorComplete(rhs, locRhs, ADD_VALUES, ierr)

Modified: short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc
===================================================================
--- short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc	2006-06-21 00:34:54 UTC (rev 3828)
+++ short/3D/PyLith/branches/pylith-0.8/pylith3d/module/scanner.cc	2006-06-21 01:42:06 UTC (rev 3829)
@@ -502,7 +502,7 @@
         const ALE::Mesh::foliation_type::index_type& idx = boundaries->getIndex(bPatch, *e_iter);
 
         if (idx.index > 0) {
-          values[c] = boundaries->restrict(bPatch, *e_iter)[0];
+          values[c] = 0.0;
         } else if (dim > 0) {
           values[c] = array[v++];
         }
@@ -510,10 +510,28 @@
       if (v != dim) {
         std::cout << "ERROR: Invalid size " << v << " used for " << *e_iter << " with index " << displacement->getIndex(*p_iter, *e_iter) << std::endl;
       }
-      full_displacement->update(*p_iter, *e_iter, values);
+      full_displacement->updateAdd(*p_iter, *e_iter, values);
     }
   }
+  ALE::Obj<ALE::Mesh::foliation_type::order_type::baseSequence> bdPatches = boundaries->getPatches();
 
+  for(ALE::Mesh::foliation_type::order_type::baseSequence::iterator p_iter = bdPatches->begin(); p_iter != bdPatches->end(); ++p_iter) {
+    ALE::Obj<ALE::Mesh::foliation_type::order_type::coneSequence> elements = boundaries->getPatch(*p_iter);
+
+    for(ALE::Mesh::foliation_type::order_type::coneSequence::iterator e_iter = elements->begin(); e_iter != elements->end(); ++e_iter) {
+      const double *array = full_displacement->restrict((*p_iter).first, *e_iter);
+      double        values[3];
+
+      if (boundaries->getIndex(*p_iter, *e_iter).index > 0) {
+        for(int c = 0; c < 3; c++) {
+          values[c] = array[c];
+        }
+        values[(*p_iter).second-1] = boundaries->restrict(*p_iter, *e_iter)[0];
+        full_displacement->update((*p_iter).first, *e_iter, values);
+      }
+    }
+  }
+
   PetscViewerCreate(mesh->comm(), &viewer);
   PetscViewerSetType(viewer, PETSC_VIEWER_ASCII);
   PetscViewerSetFormat(viewer, PETSC_VIEWER_ASCII_VTK);