[cig-commits] commit: Generalize V_Boundary_Refine::refine to 3D somewhat. Also skip j==j_min in Update_V_2D if we are not doing j_min+1.
Mercurial
hg at geodynamics.org
Sun Apr 24 06:02:24 PDT 2011
changeset: 196:3a20848262b4
user: Walter Landry <wlandry at caltech.edu>
date: Sat Apr 23 14:01:05 2011 -0700
files: src/V_Boundary_Refine.h src/V_Boundary_Refine/Update_V.C src/V_Boundary_Refine/Update_V_2D.C src/V_Boundary_Refine/refine.C wscript
description:
Generalize V_Boundary_Refine::refine to 3D somewhat. Also skip j==j_min in Update_V_2D if we are not doing j_min+1.
diff -r 1597fc3e19a6 -r 3a20848262b4 src/V_Boundary_Refine.h
--- a/src/V_Boundary_Refine.h Sat Apr 23 13:49:14 2011 -0700
+++ b/src/V_Boundary_Refine.h Sat Apr 23 14:01:05 2011 -0700
@@ -133,7 +133,7 @@ private:
private:
std::string d_name_id;
- void Update_V
+ void Update_V_2D
(const int &axis,
const int &boundary_direction,
const bool &boundary_positive,
@@ -142,6 +142,7 @@ private:
int &i,
int &j,
const int &i_max,
+ const int &j_min,
const int &j_max,
tbox::Pointer<pdat::SideData<double> > &v,
tbox::Pointer<pdat::SideData<double> > &v_fine) const;
diff -r 1597fc3e19a6 -r 3a20848262b4 src/V_Boundary_Refine/Update_V.C
--- a/src/V_Boundary_Refine/Update_V.C Sat Apr 23 13:49:14 2011 -0700
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,157 +0,0 @@
-/*************************************************************************
- *
- * 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 side-centered double data on
- * a Cartesian mesh.
- *
- ************************************************************************/
-
-#include "V_Boundary_Refine.h"
-#include "quad_offset_interpolate.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/pdat/SideVariable.h"
-#include "SAMRAI/tbox/Utilities.h"
-#include "SAMRAI/pdat/CellData.h"
-
-#include "Boundary.h"
-
-/* This is written from the perspective of axis==x. For axis==y, we
- switch i and j and everything works out. */
-void SAMRAI::geom::V_Boundary_Refine::Update_V
-(const int &axis,
- const int &boundary_direction,
- const bool &boundary_positive,
- const pdat::SideIndex &fine,
- const hier::Index &ip, const hier::Index &jp,
- int &i,
- int &j,
- const int &i_max,
- const int &j_max,
- tbox::Pointer<SAMRAI::pdat::SideData<double> > &v,
- tbox::Pointer<SAMRAI::pdat::SideData<double> > &v_fine) const
-{
- pdat::SideIndex center(fine);
- center.coarsen(hier::Index(2,2));
-
- // tbox::plog << "Boundary Update V "
- // << axis << " "
- // << fine[0] << " "
- // << fine[1] << " "
- // << center[0] << " "
- // << center[1] << " ";
- /* Set the derivative for the normal direction */
- if(boundary_direction==axis)
- {
- /* Return early if we are at j==j_max, because that is a corner
- that we do not care about. */
- if(j==j_max)
- return;
- /* Compute the derivative at the nearest three coarse points and
- then interpolate */
-
- const double dv_plus=(*v)(center+jp+ip)-(*v)(center+jp-ip);
- const double dv_minus=(*v)(center-jp+ip)-(*v)(center-jp-ip);
- const double dv_center=(*v)(center+ip)-(*v)(center-ip);
-
- double dv_fine_minus, dv_fine_plus;
-
- quad_offset_interpolate(dv_plus,dv_minus,dv_center,
- dv_fine_plus,dv_fine_minus);
-
- hier::Index offset(ip);
-
- if(boundary_positive)
- {
- offset[axis]=2;
- }
- else
- {
- offset[axis]=-2;
- dv_fine_minus=-dv_fine_minus;
- dv_fine_plus=-dv_fine_plus;
- }
-
- if(j%2==0)
- {
- (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_minus/2;
- (*v_fine)(fine+jp)=(*v_fine)(fine-offset+jp) + dv_fine_plus/2;
- /* 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;
- }
- else
- {
- (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_plus/2;
- }
- // tbox::plog << offset[0] << " "
- // << offset[1] << " "
- // << (*v_fine)(fine) << " "
- // << (*v_fine)(fine+jp) << " "
- // << (*v_fine)(fine-offset) << " "
- // << (*v_fine)(fine-offset+jp) << " "
- // << (*v)(center) << " "
- // << (*v)(center+jp) << " "
- // << (*v)(center-jp) << " "
- // // << (&(*v)(center-jp)) << " "
- // // << v_plus << " "
- // // << v_minus << " "
- // << dv_fine_plus << " "
- // << dv_fine_minus << " "
- // << "\n";
-
- }
- /* Set the value for the tangential direction */
- else
- {
- double v_center, v_plus;
-
- v_center=
- quad_offset_interpolate((*v)(center-jp),(*v)(center),(*v)(center+jp));
-
- if(i%2==0)
- {
- (*v_fine)(fine)=v_center;
-
- // tbox::plog << (*v_fine)(fine) << " ";
- // // << (*v_fine)(fine+jp) << " "
- // // << (*v_fine)(fine-offset) << " "
- // // << (*v_fine)(fine-offset+jp) << " "
- // // << (*v)(center) << " "
- // // << (*v)(center+jp) << " "
- // // << (*v)(center-jp) << " "
- // // << dv_fine_plus << " "
- // // << dv_fine_minus << " "
-
- if(i<i_max)
- {
- v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
- (*v)(center+ip+jp));
- (*v_fine)(fine+ip)=(v_center+v_plus)/2;
-
- // tbox::plog << (*v_fine)(fine+ip) << " ";
-
- /* Since we update two points on 'i' at once, we increment 'i' again.
- This is ok, since the box in the 'j' direction is defined to be
- only one cell wide */
- ++i;
- }
- }
- else
- {
- v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
- (*v)(center+ip+jp));
- (*v_fine)(fine)=(v_center+v_plus)/2;
- // tbox::plog << (*v_fine)(fine) << " ";
- }
- // tbox::plog << "\n";
- }
-}
diff -r 1597fc3e19a6 -r 3a20848262b4 src/V_Boundary_Refine/Update_V_2D.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/V_Boundary_Refine/Update_V_2D.C Sat Apr 23 14:01:05 2011 -0700
@@ -0,0 +1,106 @@
+#include "V_Boundary_Refine.h"
+#include "quad_offset_interpolate.h"
+#include "Boundary.h"
+
+/* This is written from the perspective of axis==x. For axis==y, we
+ switch i and j and everything works out. */
+void SAMRAI::geom::V_Boundary_Refine::Update_V_2D
+(const int &axis,
+ const int &boundary_direction,
+ const bool &boundary_positive,
+ const pdat::SideIndex &fine,
+ const hier::Index &ip, const hier::Index &jp,
+ int &i,
+ int &j,
+ const int &i_max,
+ const int &j_min,
+ const int &j_max,
+ tbox::Pointer<SAMRAI::pdat::SideData<double> > &v,
+ tbox::Pointer<SAMRAI::pdat::SideData<double> > &v_fine) const
+{
+ pdat::SideIndex center(fine);
+ center.coarsen(hier::Index(2,2));
+
+ /* Set the derivative for the normal direction */
+ if(boundary_direction==axis)
+ {
+ /* Return early if we are at j==j_max, because that is a corner
+ that we do not care about. We could also skip if j==j_min as
+ long as we do not have to do j_min+1. We do not really have
+ to skip these since we are guaranteed to have valid data for
+ those "past the end" points since they are needed for
+ pressure refinement. */
+ // if(j==j_max || (j==j_min && j%2!=0))
+ if(j==j_max)
+ return;
+ /* Compute the derivative at the nearest three coarse points and
+ then interpolate */
+
+ const double dv_plus=(*v)(center+jp+ip)-(*v)(center+jp-ip);
+ const double dv_minus=(*v)(center-jp+ip)-(*v)(center-jp-ip);
+ const double dv_center=(*v)(center+ip)-(*v)(center-ip);
+
+ double dv_fine_minus, dv_fine_plus;
+
+ quad_offset_interpolate(dv_plus,dv_minus,dv_center,
+ dv_fine_plus,dv_fine_minus);
+
+ hier::Index offset(ip);
+
+ if(boundary_positive)
+ {
+ offset[axis]=2;
+ }
+ else
+ {
+ offset[axis]=-2;
+ dv_fine_minus=-dv_fine_minus;
+ dv_fine_plus=-dv_fine_plus;
+ }
+
+ if(j%2==0)
+ {
+ (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_minus/2;
+ (*v_fine)(fine+jp)=(*v_fine)(fine-offset+jp) + dv_fine_plus/2;
+ /* 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;
+ }
+ else
+ {
+ (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_plus/2;
+ }
+ }
+ /* Set the value for the tangential direction */
+ else
+ {
+ double v_center, v_plus;
+
+ v_center=
+ quad_offset_interpolate((*v)(center-jp),(*v)(center),(*v)(center+jp));
+
+ if(i%2==0)
+ {
+ (*v_fine)(fine)=v_center;
+
+ if(i<i_max)
+ {
+ v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
+ (*v)(center+ip+jp));
+ (*v_fine)(fine+ip)=(v_center+v_plus)/2;
+
+ /* Since we update two points on 'i' at once, we increment 'i' again.
+ This is ok, since the box in the 'j' direction is defined to be
+ only one cell wide */
+ ++i;
+ }
+ }
+ else
+ {
+ v_plus=quad_offset_interpolate((*v)(center+ip-jp),(*v)(center+ip),
+ (*v)(center+ip+jp));
+ (*v_fine)(fine)=(v_center+v_plus)/2;
+ }
+ }
+}
diff -r 1597fc3e19a6 -r 3a20848262b4 src/V_Boundary_Refine/refine.C
--- a/src/V_Boundary_Refine/refine.C Sat Apr 23 13:49:14 2011 -0700
+++ b/src/V_Boundary_Refine/refine.C Sat Apr 23 14:01:05 2011 -0700
@@ -10,14 +10,6 @@
************************************************************************/
#include "V_Boundary_Refine.h"
-
-#include <float.h>
-#include <math.h>
-#include "SAMRAI/geom/CartesianPatchGeometry.h"
-#include "SAMRAI/hier/Index.h"
-#include "SAMRAI/pdat/SideData.h"
-#include "SAMRAI/pdat/SideVariable.h"
-#include "SAMRAI/tbox/Utilities.h"
void SAMRAI::geom::V_Boundary_Refine::refine(
hier::Patch& fine,
@@ -50,8 +42,9 @@ void SAMRAI::geom::V_Boundary_Refine::re
const hier::IntVector& ratio,
const int &axis) const
{
- const tbox::Dimension& dim(getDim());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, overlap_box, ratio);
+ const tbox::Dimension& dimension(getDim());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dimension, fine, coarse, overlap_box, ratio);
+ const int dim(dimension.getValue());
tbox::Pointer<pdat::SideData<double> >
v = coarse.getPatchData(src_component);
@@ -64,109 +57,92 @@ void SAMRAI::geom::V_Boundary_Refine::re
TBOX_ASSERT(v->getDepth() == 1);
#endif
- hier::Box coarse_box=coarse.getBox();
hier::Box fine_box=fine.getBox();
/* We have to infer where the boundary is from the boxes */
int boundary_direction;
bool boundary_positive(false);
- if(std::abs(overlap_box.lower(0)-overlap_box.upper(0))==(axis==0 ? 1 : 0))
+
+ for(int d=0;d<dim;++d)
{
- boundary_direction=0;
- if(fine_box.upper(0)<=overlap_box.lower(0))
- boundary_positive=true;
- else if(fine_box.lower(0)>=overlap_box.upper(0))
- boundary_positive=false;
+ if(std::abs(overlap_box.lower(d)-overlap_box.upper(d))==(axis==d ? 1 : 0))
+ {
+ boundary_direction=d;
+ if(fine_box.upper(d)<=overlap_box.lower(d))
+ boundary_positive=true;
+ else if(fine_box.lower(d)>=overlap_box.upper(d))
+ boundary_positive=false;
+ else
+ abort();
+ break;
+ }
+ }
+
+ hier::Index p_min(overlap_box.lower()), p_max(overlap_box.upper());
+
+ if(boundary_direction==axis)
+ {
+ if(boundary_positive)
+ {
+ p_min[axis]=p_max[axis];
+ }
else
- abort();
+ {
+ p_max[axis]=p_min[axis];
+ }
}
- else if(std::abs(overlap_box.lower(1)-overlap_box.upper(1))==
- (axis==1 ? 1 : 0))
+
+ hier::Index ip(hier::Index::getZeroIndex(dimension)), jp(ip), kp(ip);
+ ip[0]=1;
+ jp[1]=1;
+ if(dim>2)
+ kp[2]=1;
+ // hier::Index pp[]={ip,jp,kp};
+
+ if(dim==2)
{
- boundary_direction=1;
- if(fine_box.upper(1)<=overlap_box.lower(1))
- boundary_positive=true;
- else if(fine_box.lower(1)>=overlap_box.upper(1))
- boundary_positive=false;
- else
- abort();
+ for(int j=p_min[1]; j<=p_max[1]; ++j)
+ for(int i=p_min[0]; i<=p_max[0]; ++i)
+ {
+ pdat::SideIndex fine(hier::Index(i,j),axis,pdat::SideIndex::Lower);
+ switch(axis)
+ {
+ case 0:
+ Update_V_2D(axis,boundary_direction,boundary_positive,fine,
+ ip,jp,i,j,p_max[0],p_min[1],p_max[1],v,v_fine);
+ break;
+ case 1:
+ Update_V_2D(axis,boundary_direction,boundary_positive,fine,
+ jp,ip,j,i,p_max[1],p_min[0],p_max[0],v,v_fine);
+ break;
+ default:
+ abort();
+ break;
+ }
+ }
}
else
{
- abort();
+ for(int k=p_min[2]; k<=p_max[2]; k=(k/2)*2+2)
+ 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::SideIndex fine(hier::Index(i,j,k),axis,
+ pdat::SideIndex::Lower);
+ switch(axis)
+ {
+ case 0:
+ // Update_V_3D(axis,boundary_direction,boundary_positive,fine,
+ // ip,jp,i,j,p_max[0],p_max[1],v,v_fine);
+ break;
+ case 1:
+ // Update_V_3D(axis,boundary_direction,boundary_positive,fine,
+ // jp,ip,j,i,p_max[1],p_max[0],v,v_fine);
+ break;
+ default:
+ abort();
+ break;
+ }
+ }
}
-
- int i_min(overlap_box.lower(0)), i_max(overlap_box.upper(0)),
- j_min(overlap_box.lower(1)), j_max(overlap_box.upper(1));
- if(axis==0)
- {
- if(boundary_direction==0)
- {
- if(boundary_positive)
- {
- i_min=i_max;
- }
- else
- {
- i_max=i_min;
- }
- }
- }
- else if(axis==1)
- {
- if(boundary_direction==1)
- {
- if(boundary_positive)
- {
- j_min=j_max;
- }
- else
- {
- j_max=j_min;
- }
- }
- }
-
- // tbox::plog << "VBR "
- // << fine.getPatchLevelNumber() << " "
- // << axis << " "
- // << boundary_direction << " "
- // << std::boolalpha
- // << boundary_positive << " "
- // << coarse_box.lower(0) << " "
- // << coarse_box.upper(0) << " "
- // << coarse_box.lower(1) << " "
- // << coarse_box.upper(1) << " "
- // << fine_box.lower(0) << " "
- // << fine_box.upper(0) << " "
- // << fine_box.lower(1) << " "
- // << fine_box.upper(1) << " "
-
- // << overlap_box.lower(0) << " "
- // << overlap_box.upper(0) << " "
- // << overlap_box.lower(1) << " "
- // << overlap_box.upper(1) << " "
- // << i_min << " "
- // << i_max << " "
- // << j_min << " "
- // << j_max << " "
- // << "\n";
-
- for(int j=j_min; j<=j_max; ++j)
- for(int i=i_min; i<=i_max; ++i)
- {
- pdat::SideIndex fine(hier::Index(i,j),axis,pdat::SideIndex::Lower);
- hier::Index ip(1,0), jp(0,1);
-
- if(axis==0)
- {
- Update_V(axis,boundary_direction,boundary_positive,fine,
- ip,jp,i,j,i_max,j_max,v,v_fine);
- }
- else if(axis==1)
- {
- Update_V(axis,boundary_direction,boundary_positive,fine,
- jp,ip,j,i,j_max,i_max,v,v_fine);
- }
- }
}
diff -r 1597fc3e19a6 -r 3a20848262b4 wscript
--- a/wscript Sat Apr 23 13:49:14 2011 -0700
+++ b/wscript Sat Apr 23 14:01:05 2011 -0700
@@ -28,7 +28,7 @@ def build(bld):
'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_Boundary_Refine/Update_V_2D.C',
'src/V_Coarsen_Patch_Strategy/postprocessCoarsen_2D.C',
'src/V_Coarsen_Patch_Strategy/postprocessCoarsen_3D.C',
'src/Cell_Viscosity_Coarsen.C',
More information about the CIG-COMMITS
mailing list