[cig-commits] commit: Adaptive seems to be working
Mercurial
hg at geodynamics.org
Fri Feb 25 14:16:08 PST 2011
changeset: 73:ac5b255663cf
user: Walter Landry <wlandry at caltech.edu>
date: Wed Feb 09 06:44:45 2011 -0800
files: Makefile P_Boundary_Refine.C P_Boundary_Refine.h P_Refine.C P_Refine_Patch_Strategy.C P_Refine_Patch_Strategy.h StokesFACOps.h StokesFACOps/StokesFACOps.C StokesFACOps/Update_V.C StokesFACOps/computeCompositeResidualOnLevel.C StokesFACOps/initializeOperatorState.C StokesFACOps/prolongErrorAndCorrect.C StokesFACOps/restrictSolution.C StokesFACOps/set_boundaries.C StokesFACOps/smoothErrorByRedBlack.C V_Boundary_Refine.h V_Boundary_Refine/Update_V.C V_Boundary_Refine/refine.C V_Coarsen.C V_Coarsen_Patch_Strategy.C V_Coarsen_Patch_Strategy.h V_Refine.C V_Refine_Patch_Strategy.C example_inputs/const_refine.2d.input main.C quad_offset_interpolate.h set_V_boundary.C set_V_boundary.h
description:
Adaptive seems to be working
diff -r 78f870959148 -r ac5b255663cf Makefile
--- a/Makefile Wed Feb 02 16:33:30 2011 -0800
+++ b/Makefile Wed Feb 09 06:44:45 2011 -0800
@@ -31,10 +31,13 @@ CXX_OBJS = main.o FACStokes/FACStok
FACStokes/setupPlotter.o \
FACStokes/solveStokes.o \
P_Refine.o V_Refine.o V_Coarsen.o \
+ P_Boundary_Refine.o \
V_Boundary_Refine/refine.o \
V_Boundary_Refine/Update_V.o \
+ P_Refine_Patch_Strategy.o \
V_Refine_Patch_Strategy.o \
V_Coarsen_Patch_Strategy.o \
+ set_V_boundary.o \
StokesFACOps/StokesFACOps.o \
StokesFACOps/checkInputPatchDataIndices.o \
StokesFACOps/computeCompositeResidualOnLevel.o \
diff -r 78f870959148 -r ac5b255663cf P_Boundary_Refine.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Boundary_Refine.C Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,210 @@
+/*************************************************************************
+ *
+ * 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 "quad_offset_interpolate.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,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& overlap,
+ const hier::IntVector& ratio)
+ const
+{
+ const pdat::CellOverlap* t_overlap =
+ dynamic_cast<const pdat::CellOverlap *>(&overlap);
+ const hier::BoxList& boxes = t_overlap->getDestinationBoxList();
+ for (hier::BoxList::Iterator b(boxes); b; b++)
+ {
+ hier::Box &overlap_box=b();
+
+ const tbox::Dimension& dim(getDim());
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, overlap_box, ratio);
+
+ tbox::Pointer<pdat::CellData<double> >
+ p = coarse.getPatchData(src_component);
+ tbox::Pointer<pdat::CellData<double> >
+ p_fine = fine.getPatchData(dst_component);
+#ifdef DEBUG_CHECK_ASSERTIONS
+ TBOX_ASSERT(!p.isNull());
+ TBOX_ASSERT(!p_fine.isNull());
+ TBOX_ASSERT(p->getDepth() == p_fine->getDepth());
+ TBOX_ASSERT(p->getDepth() == 1);
+#endif
+
+ hier::Box coarse_box=coarse.getBox();
+ hier::Box fine_box=fine.getBox();
+ hier::Box gbox=p_fine->getGhostBox();
+ hier::Box coarse_gbox=p->getGhostBox();
+
+ /* We have to infer where the boundary is from the boxes */
+ int boundary_direction;
+ bool boundary_positive;
+ if(overlap_box.lower(0)-overlap_box.upper(0)==0)
+ {
+ 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;
+ else
+ abort();
+ }
+ else if(overlap_box.lower(1)-overlap_box.upper(1)==0)
+ {
+ 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();
+ }
+ else
+ {
+ abort();
+ }
+
+ /* We have to infer where the boundary is from the boxes */
+ tbox::plog << "PBR "
+ << fine.getPatchLevelNumber() << " "
+ << 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) << " "
+ << "\n";
+
+ int i_min(std::max(overlap_box.lower(0),gbox.lower(0))),
+ i_max(std::min(overlap_box.upper(0),gbox.upper(0))),
+ j_min(std::max(overlap_box.lower(1),gbox.lower(1))),
+ j_max(std::min(overlap_box.upper(1),gbox.upper(1)));
+
+ for(int j=j_min; j<=j_max; ++j)
+ for(int i=i_min; i<=i_max; ++i)
+ {
+ pdat::CellIndex fine(hier::Index(i,j));
+ hier::Index ip(1,0), jp(0,1);
+
+ switch(boundary_direction)
+ {
+ case 0:
+ if(j<j_max)
+ Update_P(fine,boundary_positive ? ip : -ip,jp,j,p,p_fine);
+ break;
+ case 1:
+ if(i<i_max)
+ Update_P(fine,boundary_positive ? jp : -jp,ip,i,p,p_fine);
+ }
+ }
+ }
+}
+
+/* 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- F- |
+ j | C |F_ C | C
+ | |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))
+
+ then
+
+ f+/- = F+/- + F_ - C(i,j)
+ + (C(i-1,j-1) + C(i,j) + C(i-1,j) + C(i,j-1))/16
+
+ 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.
+
+*/
+
+
+void SAMRAI::geom::P_Boundary_Refine::Update_P
+(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
+{
+ 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));
+
+ (*p_fine)(fine)=p_minus + p_offset - (*p)(center)
+ + ((*p)(center-ip-jp) + (*p)(center)
+ - (*p)(center-ip) - (*p)(center-jp))/16;
+
+ (*p_fine)(fine+jp)=p_plus + p_offset - (*p)(center)
+ + ((*p)(center-ip+jp) + (*p)(center)
+ - (*p)(center-ip) - (*p)(center+jp))/16;
+
+ /* 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;
+
+ 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) << " "
+ << "\n";
+}
diff -r 78f870959148 -r ac5b255663cf P_Boundary_Refine.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Boundary_Refine.h Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,129 @@
+/*************************************************************************
+ *
+ * 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.
+ *
+ ************************************************************************/
+
+#ifndef included_geom_P_Boundary_Refine
+#define included_geom_P_Boundary_Refine
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefineOperator.h"
+#include "SAMRAI/hier/Box.h"
+#include "SAMRAI/hier/IntVector.h"
+#include "SAMRAI/hier/Patch.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/CellVariable.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+
+#include <string>
+
+namespace SAMRAI {
+namespace geom {
+
+/**
+ * Class P_Boundary_Refine implements the special interpolation needed
+ * for the boundary elements on the coarse-fine interface.
+ *
+ * The findRefineOperator() operator function returns true if the input
+ * variable is cell-centered double, and the std::string is "P_BOUNDARY_REFINE".
+ *
+ * @see xfer::RefineOperator
+ */
+
+class P_Boundary_Refine:
+ public xfer::RefineOperator
+{
+public:
+ /**
+ * Uninteresting default constructor.
+ */
+ explicit P_Boundary_Refine(const tbox::Dimension& dim):
+ xfer::RefineOperator(dim, "P_BOUNDARY_REFINE")
+ {
+ d_name_id = "P_BOUNDARY_REFINE";
+ }
+
+
+ /**
+ * Uninteresting virtual destructor.
+ */
+ virtual ~P_Boundary_Refine(){}
+
+ /**
+ * Return true if the variable and name std::string match cell-centered
+ * double linear interpolation; otherwise, return false.
+ */
+ bool findRefineOperator(const tbox::Pointer<hier::Variable>& var,
+ const std::string& op_name) const
+ {
+ TBOX_DIM_ASSERT_CHECK_ARGS2(*this, *var);
+
+ const tbox::Pointer<pdat::CellVariable<double> > cast_var(var);
+ if (!cast_var.isNull() && (op_name == d_name_id)) {
+ return true;
+ } else {
+ return false;
+ }
+ }
+ /**
+ * Return name std::string identifier of this refinement operator.
+ */
+ const std::string& getOperatorName() const
+ {
+ return d_name_id;
+ }
+
+ /**
+ * The priority of cell-centered double linear interpolation is 0.
+ * It will be performed before any user-defined interpolation operations.
+ */
+ int getOperatorPriority() const
+ {
+ return 0;
+ }
+
+ /**
+ * The stencil width of the linear interpolation operator is the vector
+ * of ones. That is, its stencil extends one cell outside the fine box.
+ */
+ hier::IntVector getStencilWidth() const
+ {
+ return hier::IntVector::getOne(getDim());
+ }
+
+ /**
+ * Refine the source component on the coarse patch to the destination
+ * component on the fine patch using the cell-centered double linear
+ * interpolation operator. Interpolation is performed on the intersection
+ * of the destination patch and the boxes contained in fine_overlap
+ * It is assumed that the coarse patch contains sufficient data for the
+ * stencil width of the refinement operator.
+ */
+ virtual void refine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const int dst_component,
+ const int src_component,
+ const hier::BoxOverlap& fine_overlap,
+ const hier::IntVector& ratio) const;
+
+
+private:
+ std::string d_name_id;
+
+ void Update_P(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;
+};
+
+}
+}
+#endif
diff -r 78f870959148 -r ac5b255663cf P_Refine.C
--- a/P_Refine.C Wed Feb 02 16:33:30 2011 -0800
+++ b/P_Refine.C Wed Feb 09 06:44:45 2011 -0800
@@ -47,29 +47,30 @@ void SAMRAI::geom::P_Refine::refine(
}
void SAMRAI::geom::P_Refine::refine(
- hier::Patch& fine,
- const hier::Patch& coarse,
+ hier::Patch& fine_patch,
+ const hier::Patch& coarse_patch,
const int dst_component,
const int src_component,
const hier::Box& fine_box,
const hier::IntVector& ratio) const
{
const tbox::Dimension& dim(getDim());
- TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine, coarse, fine_box, ratio);
+ TBOX_DIM_ASSERT_CHECK_DIM_ARGS4(dim, fine_patch, coarse_patch,
+ fine_box, ratio);
tbox::Pointer<pdat::CellData<double> >
- p = coarse.getPatchData(src_component);
+ p = coarse_patch.getPatchData(src_component);
tbox::Pointer<pdat::CellData<double> >
- p_fine = fine.getPatchData(dst_component);
+ p_fine = fine_patch.getPatchData(dst_component);
#ifdef DEBUG_CHECK_ASSERTIONS
TBOX_ASSERT(!p.isNull());
TBOX_ASSERT(!p_fine.isNull());
TBOX_ASSERT(p->getDepth() == p_fine->getDepth());
#endif
- hier::Box coarse_box=coarse.getBox();
+ hier::Box coarse_box=coarse_patch.getBox();
tbox::Pointer<geom::CartesianPatchGeometry>
- geom = coarse.getPatchGeometry();
+ geom = coarse_patch.getPatchGeometry();
for(int j=fine_box.lower(1); j<=fine_box.upper(1); ++j)
for(int i=fine_box.lower(0); i<=fine_box.upper(0); ++i)
@@ -129,14 +130,15 @@ void SAMRAI::geom::P_Refine::refine(
+ ((i%2==0) ? (-dp_dx) : dp_dx)
+ ((j%2==0) ? (-dp_dy) : dp_dy);
- // tbox::plog << "P_Refine "
- // << i << " "
- // << j << " "
- // << (*p_fine)(fine) << " "
- // << (*p)(center) << " "
- // << dp_dx << " "
- // << dp_dy << " "
- // << "\n";
+ tbox::plog << "P_Refine "
+ << fine_patch.getPatchLevelNumber() << " "
+ << i << " "
+ << j << " "
+ << (*p_fine)(fine) << " "
+ << (*p)(center) << " "
+ << dp_dx << " "
+ << dp_dy << " "
+ << "\n";
}
}
diff -r 78f870959148 -r ac5b255663cf P_Refine_Patch_Strategy.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Refine_Patch_Strategy.C Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,93 @@
+#include "P_Refine_Patch_Strategy.h"
+
+void
+SAMRAI::solv::P_Refine_Patch_Strategy::preprocessRefine
+(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio)
+{
+ tbox::Pointer<pdat::CellData<double> >
+ p = coarse.getPatchData(p_id);
+
+ hier::Box pbox=coarse.getBox();
+ hier::Box gbox=p->getGhostBox();
+
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = coarse.getPatchGeometry();
+ tbox::plog << "P Refine Patch preprocess "
+ << gbox.lower(0) << " "
+ << gbox.upper(0) << " "
+ << gbox.lower(1) << " "
+ << gbox.upper(1) << " "
+ << pbox.lower(0) << " "
+ << pbox.upper(0) << " "
+ << pbox.lower(1) << " "
+ << pbox.upper(1) << " "
+ << "\n";
+ for(int j=gbox.lower(1); j<=gbox.upper(1); ++j)
+ for(int i=gbox.lower(0); i<=gbox.upper(0); ++i)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+ hier::Index ip(1,0), jp(0,1);
+
+ /* This is a bit complicated because we have to deal with
+ corners. I have a feeling that corners are impossible on
+ these boundary boxes, but this is just to be safe. */
+
+ /* x=0 */
+ if(i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
+ {
+ /* y=0 */
+ if(j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ {
+ (*p)(center)=4*(*p)(center+ip+jp) + (*p)(center+ip*2+jp*2)
+ - 2*((*p)(center+ip*2+jp) + (*p)(center+ip+jp*2));
+ }
+ /* y=y_max */
+ else if(j>pbox.upper(1) && geom->getTouchesRegularBoundary(1,1))
+ {
+ (*p)(center)=4*(*p)(center+ip-jp) + (*p)(center+ip*2-jp*2)
+ - 2*((*p)(center+ip*2-jp) + (*p)(center+ip-jp*2));
+ }
+ else
+ {
+ (*p)(center)=2*(*p)(center+ip) - (*p)(center+ip*2);
+ }
+ }
+ /* x=x_max */
+ else if(i>pbox.upper(0) && geom->getTouchesRegularBoundary(0,1))
+ {
+ /* y=0 */
+ if(j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ {
+ (*p)(center)=4*(*p)(center-ip+jp) + (*p)(center-ip*2+jp*2)
+ - 2*((*p)(center-ip*2+jp) + (*p)(center-ip+jp*2));
+ }
+ /* y=y_max */
+ else if(j>pbox.upper(1) && geom->getTouchesRegularBoundary(1,1))
+ {
+ (*p)(center)=4*(*p)(center-ip-jp) + (*p)(center-ip*2-jp*2)
+ - 2*((*p)(center-ip*2-jp) + (*p)(center-ip-jp*2));
+ }
+ else
+ {
+ (*p)(center)=2*(*p)(center-ip) - (*p)(center-ip*2);
+ }
+ }
+ /* These boundaries are simpler because we have already
+ handled the corners */
+ /* y=0 */
+ else if(j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ {
+ (*p)(center)=2*(*p)(center+jp) - (*p)(center+jp*2);
+ }
+ /* y=y_max */
+ else if(j>pbox.upper(1) && geom->getTouchesRegularBoundary(1,1))
+ {
+ (*p)(center)=2*(*p)(center-jp) - (*p)(center-jp*2);
+ }
+ }
+}
diff -r 78f870959148 -r ac5b255663cf P_Refine_Patch_Strategy.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/P_Refine_Patch_Strategy.h Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,404 @@
+/*************************************************************************
+ *
+ * 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: Robin boundary condition support on cartesian grids.
+ *
+ ************************************************************************/
+#ifndef included_solv_P_Refine_Patch_Strategy
+#define included_solv_P_Refine_Patch_Strategy
+
+#include "SAMRAI/SAMRAI_config.h"
+
+#include "SAMRAI/xfer/RefinePatchStrategy.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/pdat/CellData.h"
+#include "SAMRAI/pdat/CellIndex.h"
+#include "Boundary.h"
+
+namespace SAMRAI {
+namespace solv {
+
+/*!
+ * @brief Helper utility for setting boundary conditions on P.
+ *
+ * This class inherits and implements virtual functions from
+ * xfer::RefinePatchStrategy so it may be used to help create
+ * communication schedules if desired.
+ */
+class P_Refine_Patch_Strategy:
+ public xfer::RefinePatchStrategy
+{
+
+public:
+ /*!
+ * @brief Constructor using.
+ *
+ * @param object_name Name of the object, for general referencing.
+ * @param coef_strategy Coefficients strategy being helped.
+ */
+ P_Refine_Patch_Strategy(
+ const tbox::Dimension& dim,
+ std::string object_name = std::string()):
+ xfer::RefinePatchStrategy(dim), d_dim(dim),d_object_name(object_name) {}
+
+ /*!
+ * @brief Destructor.
+ */
+ virtual ~P_Refine_Patch_Strategy(void) {}
+
+ //@{ @name xfer::RefinePatchStrategy virtuals
+
+ virtual void
+ setPhysicalBoundaryConditions(
+ hier::Patch& patch,
+ const double fill_time,
+ const hier::IntVector& ghost_width_to_fill) {}
+ hier::IntVector
+ getRefineOpStencilWidth() const
+ { return hier::IntVector::getOne(d_dim); }
+ // virtual void
+ // preprocessRefineBoxes(
+ // hier::Patch& fine,
+ // const hier::Patch& coarse,
+ // const hier::BoxList& fine_boxes,
+ // const hier::IntVector& ratio) {}
+ virtual void
+ preprocessRefine(hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio);
+
+ virtual void
+ postprocessRefineBoxes(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::BoxList& fine_boxes,
+ const hier::IntVector& ratio) {}
+ virtual void
+ postprocessRefine(
+ hier::Patch& fine,
+ const hier::Patch& coarse,
+ const hier::Box& fine_box,
+ const hier::IntVector& ratio) {}
+
+ //@}
+
+ //@{
+
+ /*!
+ * @name Functions to set boundary condition values
+ */
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the first ghost cells.
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::RefinePatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to cell-centered variables.
+ * -# Only one ghost cell width is filled. Setting a Robin
+ * boundary condition for cell-centered quantities requires
+ * only one ghost cell to be set.
+ * (More ghost cells can be filled by continuing the linear
+ * distribution of data beyond the first cell, but that is
+ * not implemented at this time.)
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to
+ * get the coefficients, then it sets the values in the first
+ * ghost cell on the boundary.
+ *
+ * To determine the value for the ghost cell,
+ * a @em linear approximation in the direction normal to the
+ * boundary is assumed. We write the following discrete
+ * approximations:
+ * @f[ u_b = \frac{ u_i + u_o }{2} @f]
+ * @f[ [u_n]_b = \frac{ u_o - u_i }{h} @f]
+ * where the subscript b stands for the the boundary,
+ * i stands for the first cell inside the boundary and
+ * o stands for the first cell outside the boundary
+ * and h is the grid spacing normal to the boundary.
+ * Applying this to the Robin formula gives
+ * @f[ u_o = \frac{ h\gamma + u_i( \beta - \frac{h}{2} \alpha ) }
+ * { \beta + \frac{h}{2} \alpha } @f] or equivalently
+ * @f[ u_o = \frac{ hg + u_i (1-a(1+\frac{h}{2})) }{ 1-a(1-\frac{h}{2}) } @f]
+ *
+ * After setting the edge (face in 3D) boundary conditions,
+ * linear approximations are used to set the boundary conditions
+ * of higher boundary types (nodes in 2D, edges and nodes in 3D).
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required in solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @internal To be more general to other data types,
+ * there could be versions for other data types also,
+ * such as ...InNodes, ...InFaces, etc. However, I'm not
+ * sure exactly how to implement the Robin boundary condition
+ * on the faces and nodes when m != 1. Should the boundary
+ * value be set or should the first ghost value be set?
+ *
+ * @internal I have not addressed possibility of differences
+ * in chosing the discrete formulation with which to compute
+ * the boundary value. The above formulation is obviously
+ * one specific approximation, but there could be others.
+ * If anoter approximation is required, there should be
+ * another class like this or this class can offer a choice
+ * to be set by the user. I favor another class.
+ *
+ * @internal Since the data alignment can be found through
+ * the target_data_id, these types of functions may be changed
+ * to just plain setBoundaryValues or setBoundaryValuesInBoundaryBoxes
+ * since it does assume boundary boxes. This may have to be
+ * expanded to later include coarse-fine boundary boxes
+ * for more generality.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set ghost cells for an entire level.
+ *
+ * Loop through all patches on the given level and call
+ * setBoundaryValuesInCells(hier::Patch &patch,
+ * const double fill_time ,
+ * const hier::IntVector &ghost_width_to_fill ,
+ * int target_data_id ,
+ * bool homogeneous_bc=false )
+ * for each.
+ *
+ * @param level PatchLevel on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param ghost_width_to_fill Max ghost width requiring fill
+ * @param target_data_id hier::Patch data index of data to be set.
+ * This data must be a cell-centered double.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesInCells(
+ // hier::PatchLevel& level,
+ // const double fill_time,
+ // const hier::IntVector& ghost_width_to_fill,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ /*!
+ * @brief Set the physical boundary condition by setting the
+ * value of the boundary nodes.
+ *
+ * This function is not yet implemented!
+ *
+ * There are some decisions that must be made before
+ * the implementation can be written.
+ * -# Do we set the values on the boundary or one cell
+ * away from the boundary?
+ * -# What is the discrete formulation we should use
+ * to compute the value to be set?
+ *
+ * This function has an interface similar to the virtual function
+ * xfer::RefinePatchStrategy::setPhysicalBoundaryConditions(),
+ * and it may be used to help implement that function,
+ * but it does not serve the same purpose. The primary
+ * differences are:
+ * -# It specializes to node-centered variables.
+ * -# User must specify the index of the data whose ghost
+ * cells need to be filled. This index is used to determine
+ * the variable for which to set the boundary coefficients
+ * and to get the data to be set.
+ *
+ * This function calls RobinBcStrategy::setBcCoefs() to get the
+ * coefficients, then it sets the values at the boundary nodes.
+ *
+ * In some cases, the calling function wants to set the
+ * boundary condition homogeneously, with g=0.
+ * This is useful in problems where the the solution of the
+ * homogeneous problem is required to solving the inhomogeneous
+ * problem. This function respects such requests specified
+ * through the argument @c homogeneous_bc.
+ *
+ * @param patch hier::Patch on which to set boundary condition
+ * @param fill_time Solution time corresponding to filling
+ * @param target_data_id hier::Patch data index of data to be set.
+ * @param homogeneous_bc Set a homogeneous boundary condition.
+ * This means g=0 for the boundary.
+ */
+ // void
+ // setBoundaryValuesAtNodes(
+ // hier::Patch& patch,
+ // const double fill_time,
+ // int target_data_id,
+ // bool homogeneous_bc = false) const;
+
+ //@}
+
+ //@{
+ /*!
+ * @name Ways to provide the Robin bc coefficients
+ */
+
+ /*!
+ * @brief Provide an implementation of the RobinBcCoefStrategy
+ * for determining the boundary coefficients.
+ *
+ * Provide the implementation that can be used to set the
+ * Robin bc coefficients.
+ *
+ * @param coef_strategy tbox::Pointer to a concrete inmplementation of
+ * the coefficient strategy.
+ */
+ // void
+ // setCoefImplementation(
+ // const RobinBcCoefStrategy* coef_strategy);
+
+ /*!
+ * @brief Set the data id that should be filled when setting
+ * physical boundary conditions.
+ *
+ * When setPhysicalBoundaryConditions is called, the data
+ * specified will be set. This information is required because
+ * the it is not passed in through the argument list of
+ * setPhysicalBounaryConditions.
+ */
+ void setTargetDataId(int id)
+ {
+ p_id=id;
+ }
+
+ /*!
+ * @brief Set whether boundary filling should assume homogeneous
+ * conditions.
+ *
+ * In certain circumstances, only the value of a is needed, while
+ * the value of g is temporarily not required and taken to be zero.
+ * (An example is in setting the boundary condition for error
+ * value in an iterative method.) In such cases, use this function
+ * to set a flag that will cause a null pointer to be given to
+ * setBcCoefs() to indicate that fact.
+ */
+ // void
+ // setHomogeneousBc(
+ // bool homogeneous_bc);
+
+ //@}
+
+private:
+ /*!
+ * @brief Trim a boundary box so that it does not stick out
+ * past a given box.
+ *
+ * Certain boundary-related operations occur on patch data that
+ * do not or cannot extend past the edgr or corner of a patch.
+ * This function is used to trim down the parts of the boundary box
+ * that extend past those points so that a suitable index range
+ * is achieved.
+ *
+ * The boundary box trimmed must be of type 1 or 2.
+ *
+ * @param boundary_box Boundary box to be trimmed.
+ * @param limit_box hier::Box to not stick past
+ *
+ * @return New trimmed boundary box.
+ */
+ // hier::BoundaryBox
+ // trimBoundaryBox(
+ // const hier::BoundaryBox& boundary_box,
+ // const hier::Box& limit_box) const;
+
+ /*!
+ * @brief Return box describing the index space of boundary nodes
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the nodes corresponding
+ * to the input boundary box. These nodes lie on the boundary
+ * itself.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the node indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeNodeBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ /*!
+ * @brief Return box describing the index space of faces
+ * defined by a boundary box.
+ *
+ * Define a box describing the indices of the codimension 1
+ * surface corresponding to the input boundary box.
+ *
+ * The input boundary_box must be of type 1
+ * (see hier::BoundaryBox::getBoundaryType()).
+ *
+ * This is a utility function for working with the
+ * indices coresponding to a boundary box but coincide
+ * with the patch boundary.
+ *
+ * @param boundary_box input boundary box
+ * @return a box to define the face indices corresponding to
+ * boundary_box
+ */
+ // hier::Box
+ // makeFaceBoundaryBox(
+ // const hier::BoundaryBox& boundary_box) const;
+
+ const tbox::Dimension d_dim;
+
+ std::string d_object_name;
+
+ /*!
+ * @brief Coefficient strategy giving a way to get to
+ * Robin bc coefficients.
+ */
+ // const RobinBcCoefStrategy* d_coef_strategy;
+
+ /*!
+ * @brief hier::Index of target patch data when filling ghosts.
+ */
+ int p_id;
+
+ /*!
+ * @brief Whether to assumg g=0 when filling ghosts.
+ */
+ // bool d_homogeneous_bc;
+
+ /*!
+ * @brief Timers for performance measurement.
+ */
+ // tbox::Pointer<tbox::Timer> t_set_boundary_values_in_cells;
+ // tbox::Pointer<tbox::Timer> t_use_set_bc_coefs;
+};
+
+}
+}
+
+#endif // included_solv_P_Refine_Patch_Strategy
diff -r 78f870959148 -r ac5b255663cf StokesFACOps.h
--- a/StokesFACOps.h Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps.h Wed Feb 09 06:44:45 2011 -0800
@@ -42,6 +42,7 @@
#include "SAMRAI/tbox/Database.h"
#include "SAMRAI/tbox/Pointer.h"
#include "SAMRAI/tbox/Timer.h"
+#include "P_Refine_Patch_Strategy.h"
#include "V_Refine_Patch_Strategy.h"
#include "V_Coarsen_Patch_Strategy.h"
@@ -1019,6 +1020,7 @@ private:
* to set and whether we are setting data with homogeneous
* boundary condition.
*/
+ P_Refine_Patch_Strategy p_refine_patch_strategy;
V_Refine_Patch_Strategy v_refine_patch_strategy;
V_Coarsen_Patch_Strategy v_coarsen_patch_strategy;
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/StokesFACOps.C
--- a/StokesFACOps/StokesFACOps.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/StokesFACOps.C Wed Feb 09 06:44:45 2011 -0800
@@ -146,6 +146,8 @@ namespace SAMRAI {
v_nocoarse_refine_operator(),
v_nocoarse_refine_algorithm(),
v_nocoarse_refine_schedules(),
+ p_refine_patch_strategy(dim,
+ d_object_name + "::refine patch strategy"),
v_refine_patch_strategy(dim,
d_object_name + "::refine patch strategy"),
v_coarsen_patch_strategy(dim,
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/Update_V.C
--- a/StokesFACOps/Update_V.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/Update_V.C Wed Feb 09 06:44:45 2011 -0800
@@ -86,14 +86,14 @@ void SAMRAI::solv::StokesFACOps::Update_
/* If y==0 */
hier::Index offset(0,0);
bool set_boundary(false);
- if(center[off_axis]==pbox.lower(off_axis))
+ if(center[axis]==pbox.lower(axis)+1)
{
- offset[off_axis]=-1;
+ offset[axis]=-1;
set_boundary=true;
}
- else if(center[off_axis]==pbox.upper(off_axis))
+ else if(center[axis]==pbox.upper(axis)-1)
{
- offset[off_axis]=1;
+ offset[axis]=1;
set_boundary=true;
}
@@ -101,7 +101,7 @@ void SAMRAI::solv::StokesFACOps::Update_
if(set_boundary)
{
dv=(*v)(pdat::SideIndex
- (center,
+ (center-offset,
axis,
pdat::SideIndex::Lower))
- (*v)(pdat::SideIndex
@@ -143,6 +143,10 @@ void SAMRAI::solv::StokesFACOps::Update_
tbox::plog << "v " << axis << " "
+ << (*v)(pdat::SideIndex(center,
+ axis,
+ pdat::SideIndex::Lower))
+ << " "
<< (*v_rhs)(pdat::SideIndex(center,
axis,
pdat::SideIndex::Lower))
@@ -185,7 +189,7 @@ void SAMRAI::solv::StokesFACOps::Update_
(*v)(pdat::SideIndex
(center+offset,axis,
pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center,axis,
+ (*v)(pdat::SideIndex(center-offset,axis,
pdat::SideIndex::Lower))
-dv;
}
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/computeCompositeResidualOnLevel.C
--- a/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/computeCompositeResidualOnLevel.C Wed Feb 09 06:44:45 2011 -0800
@@ -74,6 +74,7 @@ void SAMRAI::solv::StokesFACOps::compute
*/
const int p_id = solution.getComponentDescriptorIndex(0);
const int v_id = solution.getComponentDescriptorIndex(1);
+ p_refine_patch_strategy.setTargetDataId(p_id);
v_refine_patch_strategy.setTargetDataId(v_id);
// v_refine_patch_strategy.setHomogeneousBc(error_equation_indicator);
@@ -264,6 +265,7 @@ void SAMRAI::solv::StokesFACOps::compute
tbox::plog << "resid "
+ << ln << " "
<< i << " "
<< j << " "
// << (*p_resid)(center) << " "
@@ -274,15 +276,15 @@ void SAMRAI::solv::StokesFACOps::compute
<< (*v_resid)(pdat::SideIndex(center,pdat::SideIndex::X,
pdat::SideIndex::Lower))
<< " "
- // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // pdat::SideIndex::Upper))
- // << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Upper))
+ << " "
// << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
// pdat::SideIndex::Lower))
// << " "
- // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- // pdat::SideIndex::Upper)))
- // << " "
+ << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Upper)))
+ << " "
<< "\n";
}
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/initializeOperatorState.C
--- a/StokesFACOps/initializeOperatorState.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/initializeOperatorState.C Wed Feb 09 06:44:45 2011 -0800
@@ -201,6 +201,8 @@ void SAMRAI::solv::StokesFACOps::initial
ln,
max_gcw);
}
+
+ v_coarsen_patch_strategy.coarse_fine=d_cf_boundary;
// #ifdef HAVE_HYPRE
// if (d_coarse_solver_choice == "hypre") {
// d_hypre_solver.initializeSolverState(d_hierarchy, d_ln_min);
@@ -252,7 +254,7 @@ void SAMRAI::solv::StokesFACOps::initial
vdb->mapIndexToVariable(d_cell_scratch_id, variable);
p_ghostfill_refine_operator =
geometry->lookupRefineOperator(variable,
- "LINEAR_REFINE");
+ "P_BOUNDARY_REFINE");
vdb->mapIndexToVariable(d_side_scratch_id, variable);
v_ghostfill_refine_operator =
@@ -434,8 +436,8 @@ void SAMRAI::solv::StokesFACOps::initial
p_ghostfill_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
dest_ln - 1,
- d_hierarchy);
- // &v_refine_patch_strategy);
+ d_hierarchy,
+ &p_refine_patch_strategy);
if (!p_ghostfill_refine_schedules[dest_ln]) {
TBOX_ERROR(d_object_name
<< ": Cannot create a refine schedule for ghost filling!\n");
@@ -444,7 +446,6 @@ void SAMRAI::solv::StokesFACOps::initial
v_ghostfill_refine_algorithm->
createSchedule(d_hierarchy->getPatchLevel(dest_ln),
dest_ln - 1,
- // d_hierarchy);
d_hierarchy,
&v_refine_patch_strategy);
if (!v_ghostfill_refine_schedules[dest_ln]) {
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/prolongErrorAndCorrect.C
--- a/StokesFACOps/prolongErrorAndCorrect.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/prolongErrorAndCorrect.C Wed Feb 09 06:44:45 2011 -0800
@@ -87,6 +87,7 @@ void SAMRAI::solv::StokesFACOps::prolong
* interior in the scratch space, then use that refined data
* to correct the fine level error.
*/
+ p_refine_patch_strategy.setTargetDataId(d_cell_scratch_id);
v_refine_patch_strategy.setTargetDataId(d_side_scratch_id);
// v_refine_patch_strategy.setHomogeneousBc(true);
xeqScheduleProlongation(d_cell_scratch_id,
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/restrictSolution.C
--- a/StokesFACOps/restrictSolution.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/restrictSolution.C Wed Feb 09 06:44:45 2011 -0800
@@ -69,6 +69,7 @@ void SAMRAI::solv::StokesFACOps::restric
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(dest_ln);
set_boundaries(v_dst,level);
// v_refine_patch_strategy.setHomogeneousBc(false);
+ p_refine_patch_strategy.setTargetDataId(d.getComponentDescriptorIndex(0));
v_refine_patch_strategy.setTargetDataId(d.getComponentDescriptorIndex(1));
if (dest_ln == d_ln_min) {
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/set_boundaries.C
--- a/StokesFACOps/set_boundaries.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/set_boundaries.C Wed Feb 09 06:44:45 2011 -0800
@@ -41,12 +41,9 @@
#include "SAMRAI/xfer/PatchLevelFullFillPattern.h"
#include "Boundary.h"
-/*
-********************************************************************
-* Workhorse function to smooth error using red-black *
-* Gauss-Seidel iterations. *
-********************************************************************
-*/
+#include "set_V_boundary.h"
+
+/* Set the physical boundaries for the velocity. */
void SAMRAI::solv::StokesFACOps::set_boundaries
(const int &v_id, tbox::Pointer<hier::PatchLevel> &level)
@@ -54,15 +51,16 @@ void SAMRAI::solv::StokesFACOps::set_bou
for (hier::PatchLevel::Iterator pi(*level); pi; pi++)
{
tbox::Pointer<hier::Patch> patch = *pi;
+ set_V_boundary(*patch,v_id);
- tbox::Pointer<pdat::SideData<double> >
- v = patch->getPatchData(v_id);
+ // tbox::Pointer<pdat::SideData<double> >
+ // v = patch->getPatchData(v_id);
- hier::Box pbox=patch->getBox();
- tbox::Pointer<geom::CartesianPatchGeometry>
- geom = patch->getPatchGeometry();
+ // hier::Box pbox=patch->getBox();
+ // tbox::Pointer<geom::CartesianPatchGeometry>
+ // geom = patch->getPatchGeometry();
- hier::Box gbox=v->getGhostBox();
+ // hier::Box gbox=v->getGhostBox();
// tbox::plog << "boundary "
// << gbox.lower(0) << " "
@@ -74,87 +72,89 @@ void SAMRAI::solv::StokesFACOps::set_bou
// << pbox.lower(1) << " "
// << pbox.upper(1) << " "
// << "\n";
- for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
- for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
- {
- pdat::CellIndex center(tbox::Dimension(2));
- center[0]=i;
- center[1]=j;
- hier::Index ip(1,0), jp(0,1);
+ // for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ // for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ // {
+ // pdat::CellIndex center(tbox::Dimension(2));
+ // center[0]=i;
+ // center[1]=j;
+ // hier::Index ip(1,0), jp(0,1);
- /* vx */
- if(j<=gbox.upper(1))
- {
- /* Set a sentinel value */
- if(i<pbox.lower(0) || i>pbox.upper(0)+1)
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- boundary_value;
- }
- /* Set the value so the derivative=0 */
- else if(j<pbox.lower(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
- }
- else if(j>pbox.upper(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
- }
- tbox::plog << "set boundary vx "
- << level->getLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
- }
- /* vy */
- if(i<=gbox.upper(0))
- {
- if(j<pbox.lower(1) || j>pbox.upper(1)+1)
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- boundary_value;
- }
- else if(i<pbox.lower(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- else if(i>pbox.upper(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- tbox::plog << "set boundary vy "
- << level->getLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
- }
- }
+ // /* vx */
+ // if(j<=gbox.upper(1))
+ // {
+ // /* Set a sentinel value */
+ // if((i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
+ // || (i>pbox.upper(0)+1
+ // && geom->getTouchesRegularBoundary(0,1)))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // /* Set the value so the derivative=0 */
+ // else if(j<pbox.lower(0)
+ // && geom->getTouchesRegularBoundary(1,0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // }
+ // else if(j>pbox.upper(0)
+ // && geom->getTouchesRegularBoundary(1,1))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // }
+ // // tbox::plog << "set bc vx "
+ // // << level->getLevelNumber() << " "
+ // // << i << " "
+ // // << j << " "
+ // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // // pdat::SideIndex::Lower))
+ // // << " "
+ // // << "\n";
+ // }
+ // /* vy */
+ // if(i<=gbox.upper(0))
+ // {
+ // if((j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ // || (j>pbox.upper(1)+1
+ // && geom->getTouchesRegularBoundary(1,1)))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // else if(i<pbox.lower(0)
+ // && geom->getTouchesRegularBoundary(0,0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // else if(i>pbox.upper(0)
+ // && geom->getTouchesRegularBoundary(0,1))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // // tbox::plog << "set bc vy "
+ // // << level->getLevelNumber() << " "
+ // // << i << " "
+ // // << j << " "
+ // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // // pdat::SideIndex::Lower))
+ // // << " "
+ // // << "\n";
+ // }
+ // }
}
}
diff -r 78f870959148 -r ac5b255663cf StokesFACOps/smoothErrorByRedBlack.C
--- a/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 02 16:33:30 2011 -0800
+++ b/StokesFACOps/smoothErrorByRedBlack.C Wed Feb 09 06:44:45 2011 -0800
@@ -72,28 +72,29 @@ void SAMRAI::solv::StokesFACOps::smoothE
#endif
tbox::Pointer<hier::PatchLevel> level = d_hierarchy->getPatchLevel(ln);
- {
- hier::PatchLevel::Iterator pi(*level);
- for (pi.initialize(*level); pi; pi++) {
- tbox::Pointer<hier::Patch> patch = *pi;
- tbox::Pointer<pdat::CellData<double> > p_rhs_data =
- patch->getPatchData(p_rhs_id);
- for(pdat::CellIterator ci(p_rhs_data->getBox()); ci; ci++)
- {
- pdat::CellIndex cc=ci();
- tbox::plog << "p_rhs "
- << cc[0] << " "
- << cc[1] << " "
- << (*p_rhs_data)(cc) << " "
- // << (&(*p_rhs_data)(cc)) << " "
- << "\n";
- }
- }
- }
+ // {
+ // hier::PatchLevel::Iterator pi(*level);
+ // for (pi.initialize(*level); pi; pi++) {
+ // tbox::Pointer<hier::Patch> patch = *pi;
+ // tbox::Pointer<pdat::CellData<double> > p_rhs_data =
+ // patch->getPatchData(p_rhs_id);
+ // for(pdat::CellIterator ci(p_rhs_data->getBox()); ci; ci++)
+ // {
+ // pdat::CellIndex cc=ci();
+ // tbox::plog << "p_rhs "
+ // << cc[0] << " "
+ // << cc[1] << " "
+ // << (*p_rhs_data)(cc) << " "
+ // // << (&(*p_rhs_data)(cc)) << " "
+ // << "\n";
+ // }
+ // }
+ // }
/* Only need to sync the rhs once. This sync is needed because
calculating a new pressure update requires computing in the ghost
region so that the update for the velocity inside the box will be
correct. */
+ p_refine_patch_strategy.setTargetDataId(p_id);
v_refine_patch_strategy.setTargetDataId(v_id);
xeqScheduleGhostFillNoCoarse(p_rhs_id,v_rhs_id,ln);
set_boundaries(v_id,level);
@@ -149,6 +150,25 @@ void SAMRAI::solv::StokesFACOps::smoothE
double dx = *(geom->getDx());
double dy = *(geom->getDx());
+ /* Set an array of bools that tells me whether a point
+ should set the pressure or just let it be. This is
+ needed at coarse/fine boundaries where the pressure
+ is fixed. */
+ hier::Box gbox=p->getGhostBox();
+ std::vector<bool> set_p(gbox.size(),true);
+
+ const tbox::Array<hier::BoundaryBox >&boundaries
+ =d_cf_boundary[ln]->getEdgeBoundaries(patch->getGlobalId());
+ for(int mm=0; mm<boundaries.size(); ++mm)
+ for(int j=boundaries[mm].getBox().lower(1);
+ j<=boundaries[mm].getBox().upper(1); ++j)
+ for(int i=boundaries[mm].getBox().lower(0);
+ i<=boundaries[mm].getBox().upper(0); ++i)
+ {
+ set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)+1)*(j-gbox.lower(1))]=false;
+ }
+
for(int j=pbox.lower(1); j<=pbox.upper(1)+1; ++j)
{
/* Do the red-black skip */
@@ -169,6 +189,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
tbox::plog << "smooth "
<< ln << " "
+ << sweep << " "
+ << rb << " "
<< i << " "
<< j << " "
<< pbox.lower(0) << " "
@@ -177,18 +199,8 @@ void SAMRAI::solv::StokesFACOps::smoothE
<< pbox.upper(1) << " ";
/* Update p */
- if(!((*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- ==boundary_value
- || (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Upper))
- ==boundary_value
- || (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- ==boundary_value
- || (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Upper))
- ==boundary_value))
+ if(set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)+1)*(j-gbox.lower(1))])
{
double dvx_dx=
((*v)(pdat::SideIndex(center,pdat::SideIndex::X,
@@ -215,14 +227,28 @@ void SAMRAI::solv::StokesFACOps::smoothE
<< (*p_rhs)(center) << " "
<< dvx_dx << " "
<< dvy_dy << " "
- << delta_R_continuity << " ";
+ << delta_R_continuity << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Upper)) << " "
+ << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower)) << " "
+ << dx << " ";
}
/* Update v */
- Update_V(0,j,pbox,center,left,right,down,up,p,v,v_rhs,
- maxres,dx,dy,viscosity,theta_momentum);
- Update_V(1,i,pbox,center,down,up,left,right,p,v,v_rhs,
- maxres,dy,dx,viscosity,theta_momentum);
-
+ if(set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)+1)*(j-gbox.lower(1))]
+ || j<pbox.upper(1)+1)
+ {
+ Update_V(0,j,pbox,center,left,right,down,up,p,v,v_rhs,
+ maxres,dx,dy,viscosity,theta_momentum);
+ }
+ if(set_p[(i-gbox.lower(0))
+ + (gbox.upper(0)+1)*(j-gbox.lower(1))]
+ || i<pbox.upper(0)+1)
+ {
+ Update_V(1,i,pbox,center,down,up,left,right,p,v,v_rhs,
+ maxres,dy,dx,viscosity,theta_momentum);
+ }
// << (*v)(pdat::SideIndex
// (center,
// pdat::SideIndex::X,
diff -r 78f870959148 -r ac5b255663cf V_Boundary_Refine.h
--- a/V_Boundary_Refine.h Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Boundary_Refine.h Wed Feb 09 06:44:45 2011 -0800
@@ -139,7 +139,9 @@ private:
const bool &boundary_positive,
const pdat::SideIndex &fine,
const hier::Index &ip, const hier::Index &jp,
+ int &i,
int &j,
+ const int &j_max,
tbox::Pointer<pdat::SideData<double> > &v,
tbox::Pointer<pdat::SideData<double> > &v_fine) const;
diff -r 78f870959148 -r ac5b255663cf V_Boundary_Refine/Update_V.C
--- a/V_Boundary_Refine/Update_V.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Boundary_Refine/Update_V.C Wed Feb 09 06:44:45 2011 -0800
@@ -10,6 +10,7 @@
************************************************************************/
#include "V_Boundary_Refine.h"
+#include "quad_offset_interpolate.h"
#include <float.h>
#include <math.h>
@@ -30,84 +31,85 @@ void SAMRAI::geom::V_Boundary_Refine::Up
const bool &boundary_positive,
const pdat::SideIndex &fine,
const hier::Index &ip, const hier::Index &jp,
+ int &i,
int &j,
+ 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/=2;
+ center.coarsen(hier::Index(2,2));
const int off_axis= (axis==0) ? 1 : 0;
+
+ /* Set the derivative for the normal direction */
if(boundary_direction==axis)
{
- /* Interpolate in the y direction */
- const double dv_plus=(*v)(center+jp)-(*v)(center);
- const double dv_minus=(*v)(center)-(*v)(center-jp);
+ /* 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 */
- /* Quadratic interpolation */
- double v_plus=(*v)(center)
- + (5.0/32)*dv_plus - (3.0/32)*dv_minus;
- double v_minus=(*v)(center)
- - (5.0/32)*dv_minus + (3.0/32)*dv_plus;
+ 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);
- const double epsilon(1.0e-8);
- if(std::abs(v_plus+v_minus)<=epsilon*std::abs((*v)(center)))
- {
- (*v_fine)(fine)=(*v_fine)(fine+jp)=(*v)(center);
- }
- else
- {
- (*v_fine)(fine)=v_minus*(2*(*v)(center))/(v_plus + v_minus);
- (*v_fine)(fine+jp)=v_plus*(2*(*v)(center))/(v_plus + v_minus);
- }
+ double dv_fine_minus, dv_fine_plus;
- tbox::plog << "Update V "
- << axis << " "
- << fine[0] << " "
- << fine[1] << " "
- << center[0] << " "
- << center[1] << " "
- << (*v_fine)(fine) << " "
- << (*v_fine)(fine+jp) << " "
- << (*v)(center) << " "
- << (*v)(center+jp) << " "
- << (*v)(center-jp) << " "
- << (&(*v)(center-jp)) << " "
- << v_plus << " "
- << v_minus << " "
- << dv_plus << " "
- << dv_minus << " "
- << "\n";
+ quad_offset_interpolate(dv_plus,dv_minus,dv_center,
+ dv_fine_plus,dv_fine_minus);
- /* Set the point outside of the boundary to be max_double. This
- give us a marker for whether the current value is a boundary
- condition. */
hier::Index offset(ip);
if(!boundary_positive)
{
- offset[axis]=-1;
+ offset[axis]=-2;
}
- (*v_fine)(fine+offset)=(*v_fine)(fine+offset+jp)=
- boundary_value;
+
+ (*v_fine)(fine)=(*v_fine)(fine-offset) + dv_fine_minus;
+ (*v_fine)(fine+jp)=(*v_fine)(fine-offset+jp) + dv_fine_plus;
+
+ // tbox::plog << "Update V "
+ // << axis << " "
+ // << fine[0] << " "
+ // << fine[1] << " "
+ // << center[0] << " "
+ // << center[1] << " "
+ // << (*v_fine)(fine) << " "
+ // << (*v_fine)(fine+jp) << " "
+ // << (*v)(center) << " "
+ // << (*v)(center+jp) << " "
+ // << (*v)(center-jp) << " "
+ // << (&(*v)(center-jp)) << " "
+ // << v_plus << " "
+ // << v_minus << " "
+ // << dv_plus << " "
+ // << dv_minus << " "
+ // << "\n";
+
+ /* 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;
}
+ /* Set the value for the tangential direction */
else
{
- double dv;
- /* Compute derivatives and use that to set the new vx */
- if(fine[axis]%2==0)
+ double v_center, v_plus;
+
+ v_center=
+ quad_offset_interpolate((*v)(center-jp),(*v)(center),(*v)(center+jp));
+ (*v_fine)(fine)=v_center;
+
+ if(j<j_max)
{
- dv=((*v)(center+jp) - (*v)(center))/2;
+ 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;
}
- else
- {
- dv=((*v)(center+jp) - (*v)(center)
- + (*v)(center+jp+ip) - (*v)(center+ip))/4;
- }
- hier::Index offset(jp);
- if(!boundary_positive)
- {
- offset[off_axis]=-1;
- }
- (*v_fine)(fine)=(*v_fine)(fine-offset)+dv;
+ /* 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;
}
}
diff -r 78f870959148 -r ac5b255663cf V_Boundary_Refine/refine.C
--- a/V_Boundary_Refine/refine.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Boundary_Refine/refine.C Wed Feb 09 06:44:45 2011 -0800
@@ -104,22 +104,12 @@ void SAMRAI::geom::V_Boundary_Refine::re
{
if(boundary_positive)
{
- i_max=i_min;
+ i_min=i_max;
}
else
{
- i_min=i_max;
+ i_max=i_min;
}
- j_min=std::max(j_min,fine_box.lower(1));
- j_max=std::min(j_max,fine_box.upper(1));
- }
- /* We need to shrink the box because we do not want the edges.
- Those are points that are either covered by other patches or
- are ghost points that we do not care about. */
- else
- {
- --i_max;
- ++i_min;
}
}
else if(axis==1)
@@ -128,19 +118,12 @@ void SAMRAI::geom::V_Boundary_Refine::re
{
if(boundary_positive)
{
- j_max=j_min;
+ j_min=j_max;
}
else
{
- j_min=j_max;
+ j_max=j_min;
}
- i_min=std::max(i_min,fine_box.lower(0));
- i_max=std::min(i_max,fine_box.upper(0));
- }
- else
- {
- --j_max;
- ++j_min;
}
}
@@ -178,12 +161,12 @@ void SAMRAI::geom::V_Boundary_Refine::re
if(axis==0)
{
Update_V(axis,boundary_direction,boundary_positive,fine,
- ip,jp,j,v,v_fine);
+ ip,jp,i,j,j_max,v,v_fine);
}
else if(axis==1)
{
Update_V(axis,boundary_direction,boundary_positive,fine,
- jp,ip,i,v,v_fine);
+ jp,ip,j,i,i_max,v,v_fine);
}
}
}
diff -r 78f870959148 -r ac5b255663cf V_Coarsen.C
--- a/V_Coarsen.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Coarsen.C Wed Feb 09 06:44:45 2011 -0800
@@ -22,6 +22,19 @@
#include "SAMRAI/pdat/SideVariable.h"
#include "SAMRAI/tbox/Utilities.h"
#include "Boundary.h"
+
+using namespace SAMRAI;
+
+inline void coarsen_v(const pdat::SideIndex &coarse,
+ const hier::Index &ip, const hier::Index &jp,
+ tbox::Pointer<pdat::SideData<double> > &v,
+ tbox::Pointer<pdat::SideData<double> > &v_fine )
+{
+ pdat::SideIndex center(coarse*2);
+ (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/4
+ + ((*v_fine)(center-ip) + (*v_fine)(center-ip+jp)
+ + (*v_fine)(center+ip) + (*v_fine)(center+jp+ip))/8;
+}
void SAMRAI::geom::V_Coarsen::coarsen(hier::Patch& coarse,
@@ -88,7 +101,9 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
case, (i*2,j*2), (i*2,j*2+1), (i*2-1,j*2),
(i*2-1,j*2+1), (i*2+1,j*2), (i*2+1,j*2+1).
+ The boundaries get fixed up later.
*/
+ hier::Index ip(1,0), jp(0,1);
for(int j=coarse_box.lower(1); j<=coarse_box.upper(1)+1; ++j)
for(int i=coarse_box.lower(0); i<=coarse_box.upper(0)+1; ++i)
{
@@ -97,30 +112,32 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
<< i << " "
<< j << " ";
- hier::Index ip(1,0), jp(0,1);
if(directions(0) && j!=coarse_box.upper(1)+1)
{
pdat::SideIndex coarse(hier::Index(i,j),0,
pdat::SideIndex::Lower);
pdat::SideIndex center(coarse*2);
if((i==coarse_box.lower(0)
- && (*v_fine)(center-ip)==boundary_value)
+ && cgeom->getTouchesRegularBoundary(0,0))
|| (i==coarse_box.upper(0)+1
- && (*v_fine)(center+ip)==boundary_value))
+ && cgeom->getTouchesRegularBoundary(0,1)))
{
(*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
}
else
{
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/4
- + ((*v_fine)(center-ip) + (*v_fine)(center-ip+jp)
- + (*v_fine)(center+ip) + (*v_fine)(center+ip+jp))/8;
+ coarsen_v(coarse,ip,jp,v,v_fine);
}
tbox::plog << "vx "
<< (*v)(coarse) << " "
- << coarse_box.lower(0) << " "
+ << (*v_fine)(center) << " "
+ << (*v_fine)(center+jp) << " "
<< (*v_fine)(center-ip) << " "
- << &((*v_fine)(center-ip)) << " ";
+ << (*v_fine)(center-ip+jp) << " "
+ << (*v_fine)(center+ip) << " "
+ << (*v_fine)(center+ip+jp) << " "
+ << &((*v_fine)(center+ip)) << " ";
+
}
if(directions(1) && i!=coarse_box.upper(0)+1)
{
@@ -128,17 +145,15 @@ void SAMRAI::geom::V_Coarsen::coarsen(hi
pdat::SideIndex::Lower);
pdat::SideIndex center(coarse*2);
if((j==coarse_box.lower(1)
- && (*v_fine)(center-jp)==boundary_value)
+ && cgeom->getTouchesRegularBoundary(1,0))
|| (j==coarse_box.upper(1)+1
- && (*v_fine)(center+jp)==boundary_value))
+ && cgeom->getTouchesRegularBoundary(0,0)))
{
(*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
}
else
{
- (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/4
- + ((*v_fine)(center-jp) + (*v_fine)(center-jp+ip)
- + (*v_fine)(center+jp) + (*v_fine)(center+ip+jp))/8;
+ coarsen_v(coarse,jp,ip,v,v_fine);
}
tbox::plog << "vy "
<< (*v)(coarse) << " ";
diff -r 78f870959148 -r ac5b255663cf V_Coarsen_Patch_Strategy.C
--- a/V_Coarsen_Patch_Strategy.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Coarsen_Patch_Strategy.C Wed Feb 09 06:44:45 2011 -0800
@@ -1,4 +1,5 @@
#include "V_Coarsen_Patch_Strategy.h"
+#include "set_V_boundary.h"
void
SAMRAI::solv::V_Coarsen_Patch_Strategy::preprocessCoarsen
@@ -7,115 +8,231 @@ SAMRAI::solv::V_Coarsen_Patch_Strategy::
const hier::Box& coarse_box,
const hier::IntVector& ratio)
{
+ set_V_boundary(fine,v_id);
+
+ // tbox::Pointer<pdat::SideData<double> >
+ // v = fine.getPatchData(v_id);
+
+ // hier::Box pbox=fine.getBox();
+
+ // hier::Box gbox=v->getGhostBox();
+
+ // // tbox::plog << "boundary "
+ // // << gbox.lower(0) << " "
+ // // << gbox.upper(0) << " "
+ // // << gbox.lower(1) << " "
+ // // << gbox.upper(1) << " "
+ // // << pbox.lower(0) << " "
+ // // << pbox.upper(0) << " "
+ // // << pbox.lower(1) << " "
+ // // << pbox.upper(1) << " "
+ // // << "\n";
+ // for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ // for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ // {
+ // pdat::CellIndex center(tbox::Dimension(2));
+ // center[0]=i;
+ // center[1]=j;
+ // hier::Index ip(1,0), jp(0,1);
+
+ // /* vx */
+ // if(j<=gbox.upper(1))
+ // {
+ // /* Set a sentinel value */
+ // if(i<pbox.lower(0) || i>pbox.upper(0)+1)
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // /* Set the value so the derivative=0 */
+ // else if(j<pbox.lower(0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // tbox::plog << "V Coarsen Patch vx lower "
+ // << fine.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << gbox.lower(0) << " "
+ // << pbox.lower(0) << " "
+ // << "\n";
+ // }
+ // else if(j>pbox.upper(0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // }
+ // tbox::plog << "V Coarsen Patch vx "
+ // << fine.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower)))
+ // << " "
+ // << "\n";
+ // }
+ // /* vy */
+ // if(i<=gbox.upper(0))
+ // {
+ // if(j<pbox.lower(1) || j>pbox.upper(1)+1)
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // else if(i<pbox.lower(0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // else if(i>pbox.upper(0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // tbox::plog << "V Coarsen Patch vy "
+ // << fine.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower)))
+ // << " "
+ // << "\n";
+ // }
+ // }
+}
+
+
+void
+SAMRAI::solv::V_Coarsen_Patch_Strategy::postprocessCoarsen
+(hier::Patch& coarse,
+ const hier::Patch& fine,
+ const hier::Box& coarse_box,
+ const hier::IntVector& ratio)
+{
+ /* Fix up the boundary elements by iterating through the boundary
+ boxes */
+
+ /* We only care about edges, not corners, so we only iterate over
+ edge boundary boxes. */
+ const tbox::Array<hier::BoundaryBox>
+ &boundaries=coarse_fine[fine.getPatchLevelNumber()]->getEdgeBoundaries(coarse.getGlobalId());
+
tbox::Pointer<pdat::SideData<double> >
- v = fine.getPatchData(v_id);
+ v_fine = fine.getPatchData(v_id);
+ tbox::Pointer<pdat::SideData<double> >
+ v = coarse.getPatchData(v_id);
- hier::Box pbox=fine.getBox();
+ TBOX_ASSERT(!v.isNull());
+ TBOX_ASSERT(!v_fine.isNull());
+ TBOX_ASSERT(v_fine->getDepth() == v->getDepth());
+ TBOX_ASSERT(v->getDepth() == 1);
- hier::Box gbox=v->getGhostBox();
+ tbox::plog << "V Coarsen Patch Strategy "
+ << boundaries.size() << " "
+ << "\n";
- // tbox::plog << "boundary "
- // << gbox.lower(0) << " "
- // << gbox.upper(0) << " "
- // << gbox.lower(1) << " "
- // << gbox.upper(1) << " "
- // << pbox.lower(0) << " "
- // << pbox.upper(0) << " "
- // << pbox.lower(1) << " "
- // << pbox.upper(1) << " "
- // << "\n";
- for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
- for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
- {
- pdat::CellIndex center(tbox::Dimension(2));
- center[0]=i;
- center[1]=j;
- hier::Index ip(1,0), jp(0,1);
+ hier::Index ip(1,0), jp(0,1);
+ for(int mm=0; mm<boundaries.size(); ++mm)
+ {
+ hier::Box bbox=boundaries[mm].getBox();
+ int location_index=boundaries[mm].getLocationIndex();
- /* vx */
- if(j<=gbox.upper(1))
+ hier::Index lower=hier::Index::coarsen(bbox.lower(),hier::Index(2,2)),
+ upper=hier::Index::coarsen(bbox.upper(),hier::Index(2,2));
+
+ tbox::plog << "BBox "
+ << mm << " "
+ << bbox.lower(0) << " "
+ << bbox.upper(0) << " "
+ << bbox.lower(1) << " "
+ << bbox.upper(1) << " "
+ << location_index << " "
+ << lower(0) << " "
+ << upper(0) << " "
+ << lower(1) << " "
+ << upper(1) << " "
+ << "\n";
+
+ for(int j=lower(1); j<=upper(1); ++j)
+ for(int i=lower(0); i<=upper(0); ++i)
{
- /* Set a sentinel value */
- if(i<pbox.lower(0) || i>pbox.upper(0)+1)
+ tbox::plog << "VCPS "
+ << i << " "
+ << j << " ";
+
+ /* Fix vx */
+ if(location_index==0)
{
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- boundary_value;
+ // if(!cgeom->getTouchesRegularBoundary(0,0))
+ {
+ pdat::SideIndex coarse(hier::Index(i,j),0,
+ pdat::SideIndex::Upper);
+ pdat::SideIndex center(coarse*2);
+ (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
+ tbox::plog << (*v)(coarse) << " ";
+ }
}
- /* Set the value so the derivative=0 */
- else if(j<pbox.lower(0))
+ else if(location_index==1)
{
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
- tbox::plog << "V Coarsen Patch vx lower "
- << fine.getPatchLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- << " "
- << (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- << " "
- << gbox.lower(0) << " "
- << pbox.lower(0) << " "
- << "\n";
+ // if(!cgeom->getTouchesRegularBoundary(0,1))
+ {
+ pdat::SideIndex coarse(hier::Index(i,j),0,
+ pdat::SideIndex::Lower);
+ pdat::SideIndex center(coarse*2);
+ (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+jp))/2;
+ tbox::plog << (*v)(coarse) << " ";
+ }
}
- else if(j>pbox.upper(0))
+ /* Fix vy */
+ else if(location_index==2)
{
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
+ // if(!cgeom->getTouchesRegularBoundary(1,0))
+ {
+ pdat::SideIndex coarse(hier::Index(i,j),1,
+ pdat::SideIndex::Upper);
+ pdat::SideIndex center(coarse*2);
+ (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
+ tbox::plog << (*v)(coarse) << " ";
+ }
}
- tbox::plog << "V Coarsen Patch vx "
- << fine.getPatchLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
+ else if(location_index==3)
+ {
+ // if(!cgeom->getTouchesRegularBoundary(1,1))
+ {
+ pdat::SideIndex coarse(hier::Index(i,j),1,
+ pdat::SideIndex::Lower);
+ pdat::SideIndex center(coarse*2);
+ (*v)(coarse)=((*v_fine)(center) + (*v_fine)(center+ip))/2;
+ tbox::plog << (*v)(coarse) << " ";
+ }
+ }
+ else
+ {
+ abort();
+ }
+ tbox::plog << "\n";
}
- /* vy */
- if(i<=gbox.upper(0))
- {
- if(j<pbox.lower(1) || j>pbox.upper(1)+1)
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- boundary_value;
- }
- else if(i<pbox.lower(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- else if(i>pbox.upper(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- tbox::plog << "V Coarsen Patch vy "
- << fine.getPatchLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
- }
- }
+ }
}
diff -r 78f870959148 -r ac5b255663cf V_Coarsen_Patch_Strategy.h
--- a/V_Coarsen_Patch_Strategy.h Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Coarsen_Patch_Strategy.h Wed Feb 09 06:44:45 2011 -0800
@@ -14,6 +14,7 @@
#include "SAMRAI/xfer/CoarsenPatchStrategy.h"
#include "SAMRAI/geom/CartesianPatchGeometry.h"
+#include "SAMRAI/hier/CoarseFineBoundary.h"
#include "SAMRAI/pdat/SideData.h"
#include "SAMRAI/pdat/CellIndex.h"
#include "Boundary.h"
@@ -72,7 +73,7 @@ public:
hier::Patch& coarse,
const hier::Patch& fine,
const hier::Box& coarse_box,
- const hier::IntVector& ratio) {}
+ const hier::IntVector& ratio);
//@}
@@ -297,6 +298,8 @@ public:
//@}
+ tbox::Array<tbox::Pointer<hier::CoarseFineBoundary> > coarse_fine;
+
private:
/*!
* @brief Trim a boundary box so that it does not stick out
diff -r 78f870959148 -r ac5b255663cf V_Refine.C
--- a/V_Refine.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Refine.C Wed Feb 09 06:44:45 2011 -0800
@@ -78,14 +78,17 @@ void SAMRAI::geom::V_Refine::refine(hier
hier::Index ip(1,0), jp(0,1);
pdat::SideIndex center(fine);
- center[0]/=2;
- center[1]/=2;
+ center.coarsen(hier::Index(2,2));
/* This assumes that the levels are always properly nested,
so that we always have an extra grid space for
interpolation. So we only have to have a special case for
physical boundaries, where we do not have an extra grid
space. */
+
+ tbox::plog << "V refine "
+ << i << " "
+ << j << " ";
if(axis==0)
{
@@ -187,6 +190,10 @@ void SAMRAI::geom::V_Refine::refine(hier
+ ((i%2==0) ? (-dvy_dx) : dvy_dx))/2;
}
}
+ tbox::plog << axis << " "
+ << (*v_fine)(fine) << " "
+ << "\n";
+
}
}
#endif
diff -r 78f870959148 -r ac5b255663cf V_Refine_Patch_Strategy.C
--- a/V_Refine_Patch_Strategy.C Wed Feb 02 16:33:30 2011 -0800
+++ b/V_Refine_Patch_Strategy.C Wed Feb 09 06:44:45 2011 -0800
@@ -1,4 +1,5 @@
#include "V_Refine_Patch_Strategy.h"
+#include "set_V_boundary.h"
void
SAMRAI::solv::V_Refine_Patch_Strategy::preprocessRefine
@@ -7,13 +8,15 @@ SAMRAI::solv::V_Refine_Patch_Strategy::p
const hier::Box& fine_box,
const hier::IntVector& ratio)
{
- tbox::Pointer<pdat::SideData<double> >
- v = coarse.getPatchData(v_id);
+ // tbox::Pointer<pdat::SideData<double> >
+ // v = coarse.getPatchData(v_id);
- hier::Box pbox=coarse.getBox();
- hier::Box gbox=v->getGhostBox();
+ // hier::Box pbox=coarse.getBox();
+ // hier::Box gbox=v->getGhostBox();
- // tbox::plog << "boundary "
+ // tbox::Pointer<geom::CartesianPatchGeometry>
+ // geom = coarse->getPatchGeometry();
+ // tbox::plog << "V Refine Patch preprocess "
// << gbox.lower(0) << " "
// << gbox.upper(0) << " "
// << gbox.lower(1) << " "
@@ -23,85 +26,91 @@ SAMRAI::solv::V_Refine_Patch_Strategy::p
// << pbox.lower(1) << " "
// << pbox.upper(1) << " "
// << "\n";
- for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
- for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
- {
- pdat::CellIndex center(tbox::Dimension(2));
- center[0]=i;
- center[1]=j;
- hier::Index ip(1,0), jp(0,1);
- /* vx */
- if(j<=gbox.upper(1))
- {
- /* Set a sentinel value */
- if(i<pbox.lower(0) || i>pbox.upper(0)+1)
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- boundary_value;
- }
- /* Set the value so the derivative=0 */
- else if(j<pbox.lower(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
- }
- else if(j>pbox.upper(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
- pdat::SideIndex::Lower));
- }
- tbox::plog << "V Refine Patch vx "
- << coarse.getPatchLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
- }
- /* vy */
- if(i<=gbox.upper(0))
- {
- if(j<pbox.lower(1) || j>pbox.upper(1)+1)
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- boundary_value;
- }
- else if(i<pbox.lower(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- else if(i>pbox.upper(0))
- {
- (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))=
- (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
- pdat::SideIndex::Lower));
- }
- tbox::plog << "V Refine Patch vy "
- << coarse.getPatchLevelNumber() << " "
- << i << " "
- << j << " "
- << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower))
- << " "
- << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
- pdat::SideIndex::Lower)))
- << " "
- << "\n";
- }
- }
+ set_V_boundary(coarse,v_id);
+
+
+ // for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ // for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ // {
+ // pdat::CellIndex center(tbox::Dimension(2));
+ // center[0]=i;
+ // center[1]=j;
+ // hier::Index ip(1,0), jp(0,1);
+
+ // /* vx */
+ // if(j<=gbox.upper(1))
+ // {
+ // /* Set a sentinel value */
+ // if((i<pbox.lower(0) && geom->getTouchesPhysicalBoundary(0,0))
+ // || (i>pbox.upper(0)+1 && geom->getTouchesPhysicalBoundary(0,1)))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // /* Set the value so the derivative=0 */
+ // else if(j<pbox.lower(1) && geom->getTouchesPhysicalBoundary(1,0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // }
+ // else if(j>pbox.upper(1) && geom->getTouchesPhysicalBoundary(1,1))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower));
+ // }
+ // // tbox::plog << "V Refine Patch vx "
+ // // << coarse.getPatchLevelNumber() << " "
+ // // << i << " "
+ // // << j << " "
+ // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // // pdat::SideIndex::Lower))
+ // // << " "
+ // // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // // pdat::SideIndex::Lower)))
+ // // << " "
+ // // << "\n";
+ // }
+ // /* vy */
+ // if(i<=gbox.upper(0))
+ // {
+ // if((j<pbox.lower(1) && geom->getTouchesPhysicalBoundary(1,0))
+ // || (j>pbox.upper(1)+1 && geom->getTouchesPhysicalBoundary(1,1)))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // boundary_value;
+ // }
+ // else if(i<pbox.lower(0) && geom->getTouchesPhysicalBoundary(0,0))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // else if(i>pbox.upper(0) && geom->getTouchesPhysicalBoundary(0,1))
+ // {
+ // (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))=
+ // (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower));
+ // }
+ // // tbox::plog << "V Refine Patch vy "
+ // // << coarse.getPatchLevelNumber() << " "
+ // // << i << " "
+ // // << j << " "
+ // // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // // pdat::SideIndex::Lower))
+ // // << " "
+ // // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // // pdat::SideIndex::Lower)))
+ // // << " "
+ // // << "\n";
+ // }
+ // }
}
diff -r 78f870959148 -r ac5b255663cf example_inputs/const_refine.2d.input
--- a/example_inputs/const_refine.2d.input Wed Feb 02 16:33:30 2011 -0800
+++ b/example_inputs/const_refine.2d.input Wed Feb 09 06:44:45 2011 -0800
@@ -32,11 +32,11 @@ FACStokes {
// This is the input for the cell-centered Stokes FAC solver
// class in the SAMRAI library.
enable_logging = TRUE // Bool flag to switch logging on/off
- max_cycles = 10 // Max number of FAC cycles to use
+ max_cycles = 20 // Max number of FAC cycles to use
residual_tol = 1e-8 // Residual tolerance to solve for
- num_pre_sweeps = 0 // Number of presmoothing sweeps to use
- num_post_sweeps = 0 // Number of postsmoothing sweeps to use
- coarse_solver_max_iterations = 200
+ num_pre_sweeps = 5 // Number of presmoothing sweeps to use
+ num_post_sweeps = 5 // Number of postsmoothing sweeps to use
+ coarse_solver_max_iterations = 100
coarse_solver_tolerance = 1e-10
p_prolongation_method = "P_REFINE" // Type of refinement
// used in prolongation.
diff -r 78f870959148 -r ac5b255663cf main.C
--- a/main.C Wed Feb 02 16:33:30 2011 -0800
+++ b/main.C Wed Feb 09 06:44:45 2011 -0800
@@ -31,6 +31,7 @@ using namespace std;
#include "SAMRAI/xfer/RefineOperator.h"
#include "P_Refine.h"
#include "V_Refine.h"
+#include "P_Boundary_Refine.h"
#include "V_Boundary_Refine.h"
#include "V_Coarsen.h"
@@ -159,6 +160,9 @@ int main(
grid_geometry->addSpatialRefineOperator
(tbox::Pointer<SAMRAI::xfer::RefineOperator>
(new SAMRAI::geom::V_Refine(dim)));
+ grid_geometry->addSpatialRefineOperator
+ (tbox::Pointer<SAMRAI::xfer::RefineOperator>
+ (new SAMRAI::geom::P_Boundary_Refine(dim)));
grid_geometry->addSpatialRefineOperator
(tbox::Pointer<SAMRAI::xfer::RefineOperator>
(new SAMRAI::geom::V_Boundary_Refine(dim)));
diff -r 78f870959148 -r ac5b255663cf quad_offset_interpolate.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/quad_offset_interpolate.h Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,31 @@
+#ifndef QUAD_OFFSET_INTERPOLATE_H
+#define QUAD_OFFSET_INTERPOLATE_H
+
+/* Interpolate up to quadratic order from three coarse points (C) to
+ two fine points (f) with the setup
+
+ C--|--|--f--C--f--|--|--C
+
+*/
+
+inline void quad_offset_interpolate(const double &plus, const double ¢er,
+ const double &minus,
+ double &fine_plus, double &fine_minus)
+{
+ const double d_plus=plus-center;
+ const double d_minus=center-minus;
+
+ fine_plus=center + (5*d_plus - 3*d_minus)/32;
+ fine_minus=center + (3*d_plus - 5*d_minus)/32;
+}
+
+inline double quad_offset_interpolate(const double &plus, const double ¢er,
+ const double &minus)
+{
+ const double d_plus=plus-center;
+ const double d_minus=center-minus;
+
+ return center + (5*d_plus - 3*d_minus)/32;
+}
+
+#endif
diff -r 78f870959148 -r ac5b255663cf set_V_boundary.C
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/set_V_boundary.C Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,104 @@
+#include "set_V_boundary.h"
+#include "Boundary.h"
+#include "SAMRAI/tbox/Pointer.h"
+#include "SAMRAI/pdat/SideData.h"
+#include "SAMRAI/pdat/CellIndex.h"
+#include "SAMRAI/geom/CartesianPatchGeometry.h"
+
+using namespace SAMRAI;
+
+void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &v_id)
+{
+ tbox::Pointer<pdat::SideData<double> >
+ v = patch.getPatchData(v_id);
+
+ hier::Box pbox=patch.getBox();
+ hier::Box gbox=v->getGhostBox();
+
+ tbox::Pointer<geom::CartesianPatchGeometry>
+ geom = patch.getPatchGeometry();
+
+ for(int j=gbox.lower(1); j<=gbox.upper(1)+1; ++j)
+ for(int i=gbox.lower(0); i<=gbox.upper(0)+1; ++i)
+ {
+ pdat::CellIndex center(tbox::Dimension(2));
+ center[0]=i;
+ center[1]=j;
+ hier::Index ip(1,0), jp(0,1);
+
+ /* vx */
+ if(j<=gbox.upper(1))
+ {
+ /* Set a sentinel value */
+ if((i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
+ || (i>pbox.upper(0)+1 && geom->getTouchesRegularBoundary(0,1)))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ /* Set the value so the derivative=0 */
+ else if(j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ }
+ else if(j>pbox.upper(1) && geom->getTouchesRegularBoundary(1,1))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-jp,pdat::SideIndex::X,
+ pdat::SideIndex::Lower));
+ }
+ // tbox::plog << "V Refine Patch vx "
+ // << coarse.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::X,
+ // pdat::SideIndex::Lower)))
+ // << " "
+ // << "\n";
+ }
+ /* vy */
+ if(i<=gbox.upper(0))
+ {
+ if((j<pbox.lower(1) && geom->getTouchesRegularBoundary(1,0))
+ || (j>pbox.upper(1)+1 && geom->getTouchesRegularBoundary(1,1)))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ boundary_value;
+ }
+ else if(i<pbox.lower(0) && geom->getTouchesRegularBoundary(0,0))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center+ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ else if(i>pbox.upper(0) && geom->getTouchesRegularBoundary(0,1))
+ {
+ (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower))=
+ (*v)(pdat::SideIndex(center-ip,pdat::SideIndex::Y,
+ pdat::SideIndex::Lower));
+ }
+ // tbox::plog << "V Refine Patch vy "
+ // << coarse.getPatchLevelNumber() << " "
+ // << i << " "
+ // << j << " "
+ // << (*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower))
+ // << " "
+ // << (&(*v)(pdat::SideIndex(center,pdat::SideIndex::Y,
+ // pdat::SideIndex::Lower)))
+ // << " "
+ // << "\n";
+ }
+ }
+}
diff -r 78f870959148 -r ac5b255663cf set_V_boundary.h
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/set_V_boundary.h Wed Feb 09 06:44:45 2011 -0800
@@ -0,0 +1,8 @@
+#ifndef SET_V_BOUNDARY_H
+#define SET_V_BOUNDARY_H
+
+#include "SAMRAI/hier/Patch.h"
+
+void set_V_boundary(const SAMRAI::hier::Patch& patch, const int &id);
+
+#endif
More information about the CIG-COMMITS
mailing list