[cig-commits] commit: Get edge viscosity by interpolating from cell viscosity, and coarsen cell viscosity using CONSERVATIVE_COARSEN

[Mercurial]

changeset:   180:949397d18c76
user:        Walter Landry <wlandry at caltech.edu>
date:        Wed Apr 20 10:33:08 2011 -0700
files:       FACStokes/fix_viscosity.C
description:
Get edge viscosity by interpolating from cell viscosity, and coarsen cell viscosity using CONSERVATIVE_COARSEN


diff -r 9a2347e0208a -r 949397d18c76 FACStokes/fix_viscosity.C
--- a/FACStokes/fix_viscosity.C	Wed Apr 20 10:30:44 2011 -0700
+++ b/FACStokes/fix_viscosity.C	Wed Apr 20 10:33:08 2011 -0700
@@ -25,7 +25,8 @@ void SAMRAI::FACStokes::fix_viscosity()
   vdb->mapIndexToVariable(cell_viscosity_id, variable);
   cell_viscosity_coarsen_operator =
     geometry->lookupCoarsenOperator(variable,
-                                    "CELL_VISCOSITY_COARSEN");
+                                    "CONSERVATIVE_COARSEN");
+                                    // "CELL_VISCOSITY_COARSEN");
 
   vdb->mapIndexToVariable(edge_viscosity_id, variable);
   edge_viscosity_coarsen_operator =
@@ -100,4 +101,67 @@ void SAMRAI::FACStokes::fix_viscosity()
 
   edge_viscosity_coarsen_algorithm.setNull();
   edge_viscosity_coarsen_schedules.setNull();
+
+  hier::Index ip(hier::Index::getZeroIndex(d_dim)), jp(ip), kp(ip);
+  ip[0]=1;
+  jp[1]=1;
+  if(d_dim.getValue()>2)
+    kp[2]=1;
+
+  for (int ln = 0; ln <= d_hierarchy->getFinestLevelNumber(); ++ln) {
+    tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
+    hier::PatchLevel::Iterator i_p(*level);
+    for ( ; i_p; i_p++) {
+      tbox::Pointer<hier::Patch> patch = *i_p;
+      tbox::Pointer<pdat::CellData<double> >
+        cell_viscosity_ptr = patch->getPatchData(cell_viscosity_id);
+      pdat::CellData<double> &cell_viscosity(*cell_viscosity_ptr);
+      if(d_dim.getValue()==2)
+        {
+          tbox::Pointer<pdat::NodeData<double> >
+            edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
+          pdat::NodeData<double> &edge_viscosity(*edge_viscosity_ptr);
+
+          for(pdat::NodeIterator ni(edge_viscosity.getBox()); ni; ni++)
+            {
+              pdat::NodeIndex e=ni();
+              pdat::CellIndex c(e);
+              cell_viscosity(c);
+              cell_viscosity(c-ip);
+              cell_viscosity(c-jp);
+              cell_viscosity(c-ip-jp);
+              edge_viscosity(e)=
+                1/(1/cell_viscosity(c)
+                   + 1/cell_viscosity(c-ip)
+                   + 1/cell_viscosity(c-jp)
+                   + 1/cell_viscosity(c-ip-jp));
+            }
+        }
+      else
+        {
+          tbox::Pointer<pdat::EdgeData<double> >
+            edge_viscosity_ptr = patch->getPatchData(edge_viscosity_id);
+          pdat::EdgeData<double> &edge_viscosity(*edge_viscosity_ptr);
+          for(int axis=0;axis<3;++axis)
+            {
+              for(pdat::EdgeIterator ni(edge_viscosity.getBox(),axis); ni; ni++)
+                {
+                  pdat::EdgeIndex e=ni();
+                  pdat::CellIndex c(e);
+                  edge_viscosity(e)=
+                    1/(1/cell_viscosity(c)
+                       + 1/cell_viscosity(c-ip)
+                       + 1/cell_viscosity(c-jp)
+                       + 1/cell_viscosity(c-ip-jp)
+                       + 1/cell_viscosity(c-kp)
+                       + 1/cell_viscosity(c-ip-kp)
+                       + 1/cell_viscosity(c-jp-kp)
+                       + 1/cell_viscosity(c-ip-jp-kp));
+                }
+            }
+        }
+    }
+  }
+
+
 }