[cig-commits] commit: Fix a bug in P_Boundary_Refine and extend it to 3D
Mercurial
hg at geodynamics.org
Sat Apr 23 00:05:15 PDT 2011
changeset: 193:d343cc39a066
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sat Apr 23 00:03:59 2011 -0700
files: src/P_Boundary_Refine.h src/P_Boundary_Refine/Update_P_2D.C src/P_Boundary_Refine/Update_P_3D.C src/P_Boundary_Refine/refine.C wscript
description:
Fix a bug in P_Boundary_Refine and extend it to 3D
diff -r 91b4296b8c4c -r d343cc39a066 src/P_Boundary_Refine.h
--- a/src/P_Boundary_Refine.h Fri Apr 22 19:48:29 2011 -0700
+++ b/src/P_Boundary_Refine.h Sat Apr 23 00:03:59 2011 -0700
@@ -119,9 +119,17 @@ private:
void Update_P_2D(const pdat::CellIndex &fine,
const hier::Index &ip, const hier::Index &jp,
- int &j,
- tbox::Pointer<SAMRAI::pdat::CellData<double> > &p,
- tbox::Pointer<SAMRAI::pdat::CellData<double> > &p_fine)
+ const int &j, const int &j_max,
+ SAMRAI::pdat::CellData<double> &p,
+ SAMRAI::pdat::CellData<double> &p_fine)
+ const;
+
+ void Update_P_3D(const pdat::CellIndex &fine,
+ const hier::Index &ip, const hier::Index &jp,
+ const hier::Index &kp,
+ const int &j, const int &k, const int &j_max, const int &k_max,
+ SAMRAI::pdat::CellData<double> &p,
+ SAMRAI::pdat::CellData<double> &p_fine)
const;
};
diff -r 91b4296b8c4c -r d343cc39a066 src/P_Boundary_Refine/Update_P_2D.C
--- a/src/P_Boundary_Refine/Update_P_2D.C Fri Apr 22 19:48:29 2011 -0700
+++ b/src/P_Boundary_Refine/Update_P_2D.C Sat Apr 23 00:03:59 2011 -0700
@@ -1,33 +1,37 @@
#include "P_Boundary_Refine.h"
-#include "quad_offset_interpolate.h"
/* Interpolate the pressure from the coarse (C) to the fine (f+/-)
coordinates using the intermediate fine points (F+/-, F_).
- i-1 i i+1
+ i-1 i i+1
- ------- -------
- | | |
+ ------- -------
+ | | |
j-1 | C | C | C
- | | |
- ------- -------
- | |f- F- |
+ | | |
+ ------- -------
+ | |f- F- |
j | C |F_ C | C
- | |f+ F+ |
- ------- -------
- | | |
+ | |f+ F+ |
+ ------- -------
+ | | |
j+1 | C | C | C
- | | |
- ------- -------
+ | | |
+ ------- -------
- F+/- = interpolate(C(i,j-1),C(i,j),C(i,j+1))
- F_ = interpolate(C(i-1,j),C(i,j),C(i+1,j))
+ C= a + b*x + c*x^2 + d*y + e*y^2 + f*x*y
- then
+ C(0,0)=a
+ C(+,0)=a+b+c
+ C(-,0)=a-b+c
+ C(0,+)=a+d+e
+ C(0,-)=a-d+e
+ C(-,-)=a-b+c-d+e+f
- f+/- = F+/- + F_ - C(i,j)
- + (C(i-1,j-1) + C(i,j) + C(i-1,j) + C(i,j-1))/16
-
+ f(-,-) = a - b/4 + c/16 - d/4 + e/16 + f/16
+ = C(-,-)/16 + (15/16)*C(0,0)
+ + (3/32)*(-C(+,0) - C(0,+) + C(-,0) + C(0,-))
+
This example show a boundary in the positive x direction. To reverse
the direction, pass in ip -> -ip. To do the y direction, switch ip
and jp, and replace j with i.
@@ -38,63 +42,33 @@ void SAMRAI::geom::P_Boundary_Refine::Up
void SAMRAI::geom::P_Boundary_Refine::Update_P_2D
(const pdat::CellIndex &fine,
const hier::Index &ip, const hier::Index &jp,
- int &j,
- tbox::Pointer<SAMRAI::pdat::CellData<double> > &p,
- tbox::Pointer<SAMRAI::pdat::CellData<double> > &p_fine) const
+ const int &j, const int &j_max,
+ SAMRAI::pdat::CellData<double> &p,
+ SAMRAI::pdat::CellData<double> &p_fine) const
{
- double p_plus, p_minus, p_offset;
pdat::CellIndex center(fine);
center.coarsen(hier::Index(2,2));
- quad_offset_interpolate((*p)(center+jp),(*p)(center),(*p)(center-jp),
- p_plus,p_minus);
- p_offset=
- quad_offset_interpolate((*p)(center-ip),(*p)(center),(*p)(center+ip));
-
- const double p_low=p_minus + p_offset - (*p)(center)
- + ((*p)(center-ip-jp) + (*p)(center)
- - (*p)(center-ip) - (*p)(center-jp))/16;
-
- const double p_high=p_plus + p_offset - (*p)(center)
- + ((*p)(center-ip+jp) + (*p)(center)
- - (*p)(center-ip) - (*p)(center+jp))/16;
-
-
+
/* If we are at an even index, update both of the elements in the cell */
if(j%2==0)
{
- (*p_fine)(fine)=p_low;
-
- (*p_fine)(fine+jp)=p_high;
-
- /* Since we update two points on j at once, we increment j again.
- This is ok, since the box in the 'i' direction is defined to be
- only one cell wide */
- ++j;
+ const double p_low=p(center-ip-jp)/16 + (15.0/16)*p(center)
+ + (3.0/32)*(-p(center+ip) - p(center+jp) + p(center-ip) + p(center-jp));
+ p_fine(fine)=p_low;
+ if(j<j_max)
+ {
+ const double p_high=p(center-ip+jp)/16 + (15.0/16)*p(center)
+ + (3.0/32)*(-p(center+ip) - p(center-jp)
+ + p(center-ip) + p(center+jp));
+ p_fine(fine+jp)=p_high;
+ }
}
else
{
- (*p_fine)(fine)=p_high;
+ const double p_high=p(center-ip+jp)/16 + (15.0/16)*p(center)
+ + (3.0/32)*(-p(center+ip) - p(center-jp)
+ + p(center-ip) + p(center+jp));
+ p_fine(fine)=p_high;
}
-
- // tbox::plog << "p bc "
- // << fine[0] << " "
- // << fine[1] << " "
- // << center[0] << " "
- // << center[1] << " "
- // << jp[0] << " "
- // << jp[1] << " "
- // << (*p_fine)(fine) << " "
- // << (*p_fine)(fine+jp) << " "
- // << p_minus << " "
- // << p_plus << " "
- // << p_offset << " "
- // << (*p)(center+jp) << " "
- // << (*p)(center) << " "
- // << (*p)(center-jp) << " "
- // << (*p)(center+ip) << " "
- // << (*p)(center-ip) << " "
- // << (*p)(center-ip+jp) << " "
- // << (*p)(center-ip-jp) << " "
- // << "\n";
}
diff -r 91b4296b8c4c -r d343cc39a066 src/P_Boundary_Refine/Update_P_3D.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/P_Boundary_Refine/Update_P_3D.C Sat Apr 23 00:03:59 2011 -0700
@@ -0,0 +1,113 @@
+#include "P_Boundary_Refine.h"
+
+/* Interpolate the pressure from the coarse (C) to the fine (f+/-)
+ coordinates using the intermediate fine points (F+/-, F_).
+
+ i-1 i i+1
+
+ ------- -------
+ | | |
+ j-1 | C | C | C
+ | | |
+ ------- -------
+ | |f- Y- |
+ j | C |X_ C | C
+ | |f+ Y+ |
+ ------- -------
+ | | |
+ j+1 | C | C | C
+ | | |
+ ------- -------
+
+
+ C= a + b*x + c*x^2 + d*y + e*y^2 + f*z + g*z^2
+ + h*x*y + i*x*z + j*y*z + k*x*y*z
+
+ C(0,0,0)=a
+ C(+,0,0)=a+b+c
+ C(-,0,0)=a-b+c
+ C(0,+,0)=a+d+e
+ C(0,-,0)=a-d+e
+ C(-,-,0)=a-b+c-d+e+h
+ C(-,-,-)=a-b+c-d+e-f+g+h+i+j-k
+
+ f(-,-) = a - b/4 + c/16 - d/4 + e/16 - f/4 + g/16 + h/16 + i/16 + j/16 - k/64
+ = C(-,-)/64 + (63/64)*C(0,0)
+ - (6/64)*(C(+,0,0) + C(0,+,0) + C(0,0,+))
+ + (3/64)*(C(-,0,0) + C(0,-,0) + C(0,0,-)
+ + C(-,-,0) + C(-,0,-) + C(0,-,-))
+
+ Note that if, for example, C is constant in the Z direction, then
+ this formula reduces to the one used in Update_P_2D.
+
+ This example show a boundary in the positive x direction. To
+ reverse the direction, pass in ip -> -ip. To do the y direction,
+ rotate [ip,jp,kp], and switch j with k. To do z, rotate again and
+ switch k with i.
+*/
+
+
+void SAMRAI::geom::P_Boundary_Refine::Update_P_3D
+(const pdat::CellIndex &fine,
+ const hier::Index &ip, const hier::Index &jp, const hier::Index &kp,
+ const int &j, const int &k, const int &j_max, const int &k_max,
+ SAMRAI::pdat::CellData<double> &p,
+ SAMRAI::pdat::CellData<double> &p_fine) const
+{
+ pdat::CellIndex center(fine);
+ center.coarsen(hier::Index(2,2));
+
+ const double p_mmm=p(center-ip-jp-kp)/64 + (63.0/64)*p(center)
+ - (3.0/32) * (p(center+ip) + p(center+jp) + p(center+kp))
+ + (3.0/64) * (p(center-ip) + p(center-jp) + p(center-kp)
+ + p(center-ip-jp) + p(center-jp-kp) + p(center-ip-kp));
+
+ const double p_mmp=p(center-ip-jp+kp)/64 + (63.0/64)*p(center)
+ - (3.0/32) * (p(center+ip) + p(center+jp) + p(center-kp))
+ + (3.0/64) * (p(center-ip) + p(center-jp) + p(center+kp)
+ + p(center-ip-jp) + p(center-jp+kp) + p(center-ip+kp));
+
+ const double p_mpm=p(center-ip+jp-kp)/64 + (63.0/64)*p(center)
+ - (3.0/32) * (p(center+ip) + p(center-jp) + p(center+kp))
+ + (3.0/64) * (p(center-ip) + p(center+jp) + p(center-kp)
+ + p(center-ip+jp) + p(center+jp-kp) + p(center-ip-kp));
+
+ const double p_mpp=p(center-ip+jp+kp)/64 + (63.0/64)*p(center)
+ - (3.0/32) * (p(center+ip) + p(center-jp) + p(center-kp))
+ + (3.0/64) * (p(center-ip) + p(center+jp) + p(center+kp)
+ + p(center-ip+jp) + p(center+jp+kp) + p(center-ip+kp));
+
+ /* If we are at an even index, update both of the elements in the cell */
+ if(j%2==0)
+ {
+ if(k%2==0)
+ {
+ p_fine(fine)=p_mmm;
+ if(j<j_max)
+ p_fine(fine+jp)=p_mpm;
+ if(k<k_max)
+ p_fine(fine+kp)=p_mmp;
+ if(j<j_max && k<k_max)
+ p_fine(fine+jp+kp)=p_mpp;
+ }
+ else
+ {
+ p_fine(fine)=p_mmp;
+ if(j<j_max)
+ p_fine(fine+jp)=p_mpp;
+ }
+ }
+ else
+ {
+ if(k%2==0)
+ {
+ p_fine(fine)=p_mpm;
+ if(k<k_max)
+ p_fine(fine+kp)=p_mpp;
+ }
+ else
+ {
+ p_fine(fine)=p_mpp;
+ }
+ }
+}
diff -r 91b4296b8c4c -r d343cc39a066 src/P_Boundary_Refine/refine.C
--- a/src/P_Boundary_Refine/refine.C Fri Apr 22 19:48:29 2011 -0700
+++ b/src/P_Boundary_Refine/refine.C Sat Apr 23 00:03:59 2011 -0700
@@ -1,23 +1,4 @@
-/*************************************************************************
- *
- * This file is part of the SAMRAI distribution. For full copyright
- * information, see COPYRIGHT and COPYING.LESSER.
- *
- * Copyright: (c) 1997-2010 Lawrence Livermore National Security, LLC
- * Description: Linear refine operator for cell-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
#include "P_Boundary_Refine.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/CellData.h"
-#include "SAMRAI/pdat/CellIndex.h"
-#include "SAMRAI/tbox/Utilities.h"
void SAMRAI::geom::P_Boundary_Refine::refine(hier::Patch& fine,
const hier::Patch& coarse,
@@ -89,22 +70,27 @@ void SAMRAI::geom::P_Boundary_Refine::re
if(dim==2)
{
- for(int j=p_min[1]; j<=p_max[1]; ++j)
- for(int i=p_min[0]; i<=p_max[0]; ++i)
+ /* This odd stride is because we handle all of the fine
+ * boundary cells in a coarse cell at once. However,
+ * sometimes there is only one fine cell in a coarse cell,
+ * so the starting point does not align with the coarse
+ * cell. The stride ensures that, if we start not aligned,
+ * the next step will be aligned. */
+
+ for(int j=p_min[1]; j<=p_max[1]; j=(j/2)*2+2)
+ for(int i=p_min[0]; i<=p_max[0]; i=(i/2)*2+2)
{
pdat::CellIndex fine(hier::Index(i,j));
switch(boundary_direction)
{
case 0:
- if(j<p_max[1])
- Update_P_2D(fine,boundary_positive ? ip : -ip,jp,j,
- p,p_fine);
+ Update_P_2D(fine,boundary_positive ? ip : -ip,jp,j,p_max[1],
+ *p,*p_fine);
break;
case 1:
- if(i<p_max[0])
- Update_P_2D(fine,boundary_positive ? jp : -jp,ip,i,
- p,p_fine);
+ Update_P_2D(fine,boundary_positive ? jp : -jp,ip,i,p_max[0],
+ *p,*p_fine);
break;
}
}
@@ -119,21 +105,18 @@ void SAMRAI::geom::P_Boundary_Refine::re
switch(boundary_direction)
{
- // case 0:
- // if(j<p_max[1] && k<p_max[2])
- // Update_P_3D(fine,boundary_positive ? ip : -ip,jp,kp,j,
- // p,p_fine);
- // break;
- // case 1:
- // if(i<p_max[0] && k<p_max[2])
- // Update_P_3D(fine,boundary_positive ? jp : -jp,kp,ip,k,
- // p,p_fine);
- // break;
- // case 1:
- // if(i<p_max[0] && j<p_max[1])
- // Update_P_3D(fine,boundary_positive ? kp : -kp,ip,jp,i,
- // p,p_fine);
- // break;
+ case 0:
+ Update_P_3D(fine,boundary_positive ? ip : -ip,jp,kp,
+ j,k,p_max[1],p_max[2],*p,*p_fine);
+ break;
+ case 1:
+ Update_P_3D(fine,boundary_positive ? jp : -jp,kp,ip,
+ k,i,p_max[2],p_max[0],*p,*p_fine);
+ break;
+ case 2:
+ Update_P_3D(fine,boundary_positive ? kp : -kp,ip,jp,
+ i,j,p_max[0],p_max[1],*p,*p_fine);
+ break;
}
}
}
diff -r 91b4296b8c4c -r d343cc39a066 wscript
--- a/wscript Fri Apr 22 19:48:29 2011 -0700
+++ b/wscript Sat Apr 23 00:03:59 2011 -0700
@@ -26,6 +26,7 @@ def build(bld):
'src/V_Coarsen/coarsen_3D.C',
'src/P_Boundary_Refine/refine.C',
'src/P_Boundary_Refine/Update_P_2D.C',
+ 'src/P_Boundary_Refine/Update_P_3D.C',
'src/V_Boundary_Refine/refine.C',
'src/V_Boundary_Refine/Update_V.C',
'src/V_Coarsen_Patch_Strategy/postprocessCoarsen_2D.C',
More information about the CIG-COMMITS
mailing list