C++ hier::Patch类代码示例

// 输出网格到 VTK 文件
void MeshOpt::writeToVTK(hier::Patch<NDIM>& patch,
                         const double  time,
                         const double  dt,
                         const bool    initial_time)
{
    NULL_USE(dt);
    NULL_USE(time);
    NULL_USE(initial_time);

    const tbox::Pointer< hier::BlockPatchGeometry<NDIM> > pgeom =
            patch.getPatchGeometry();

    int block_index = pgeom->getBlockNumber();

    int patch_index = patch.getPatchNumber();

    std::stringstream bi, pi, df;
    bi << block_index;
    pi << patch_index;
    df << d_flag;

    std::string file_name = df.str() + "_block_ " + bi.str()+ "_patch_" +  pi.str()  + ".vtk";


    MsqError err;
    MeshImpl * mesh = createLocalMesh(patch);
    mesh->write_vtk(file_name.c_str(), err);

    return;
}

void PatchMultiblockTestStrategy::tagCellsInInputBoxes(
   hier::Patch& patch,
   int level_number,
   int tag_index)
{

   if (level_number < static_cast<int>(d_refine_level_boxes.size())) {

      std::shared_ptr<pdat::CellData<int> > tags(
         SAMRAI_SHARED_PTR_CAST<pdat::CellData<int>, hier::PatchData>(
            patch.getPatchData(tag_index)));
      TBOX_ASSERT(tags);
      tags->fillAll(0);

      const hier::Box pbox = patch.getBox();

      for (hier::BoxContainer::iterator k =
              d_refine_level_boxes[level_number].begin();
           k != d_refine_level_boxes[level_number].end(); ++k) {
         tags->fill(1, *k * pbox, 0);
      }

   }

}

void EdgeMultiblockTest::postprocessRefine(
   hier::Patch& fine,
   const hier::Patch& coarse,
   const std::shared_ptr<hier::VariableContext>& context,
   const hier::Box& fine_box,
   const hier::IntVector& ratio) const
{
   pdat::EdgeDoubleConstantRefine ref_op;

   hier::BoxContainer fine_box_list(fine_box);
   hier::BoxContainer empty_box_list;

   xfer::BoxGeometryVariableFillPattern fill_pattern;

   for (int i = 0; i < static_cast<int>(d_variables.size()); ++i) {

      int id = hier::VariableDatabase::getDatabase()->
         mapVariableAndContextToIndex(d_variables[i], context);

      std::shared_ptr<hier::PatchDataFactory> fine_pdf(
         fine.getPatchDescriptor()->getPatchDataFactory(id));

      std::shared_ptr<hier::BoxOverlap> fine_overlap(
         fill_pattern.computeFillBoxesOverlap(
            fine_box_list,
            empty_box_list,
            fine.getBox(),
            fine.getPatchData(id)->getGhostBox(),
            *fine_pdf));

      ref_op.refine(fine, coarse, id, id, *fine_overlap, ratio);
   }
}

// 扰动网格
void MeshOpt::disturbMesh(hier::Patch<NDIM>& patch,
                          const double  time,
                          const double  dt,
                          const bool    initial_time)
{
    NULL_USE(dt);
    NULL_USE(time);
    NULL_USE(initial_time);

    tbox::Pointer< pdat::NodeData<NDIM,double> > coords_current
            = patch.getPatchData(d_coords_current_id);

    tbox::Pointer< pdat::NodeData<NDIM,bool> > fixed_info
            = patch.getPatchData(d_fixed_info_id);

    int count = -1;
    double dist = 0.01;
    for(pdat::NodeIterator<NDIM> ic((*coords_current).getBox()); ic; ic++)
    {
        dist *= -1;
        if((*fixed_info)(ic(),0) == false)
        {
            (*coords_current)(ic(),0) += dist;
            (*coords_current)(ic(),1) -= dist;
        }
        ++count;
    }

    // 表示已扰动
    d_flag = 1;

}

/********************************************************************
 * 网格处理                                             *
 ********************************************************************/
 void MeshOpt::setNodeInfo(hier::Patch<NDIM>& patch)
{
    const hier::Index<NDIM> ifirst=patch.getBox().lower();
    const hier::Index<NDIM> ilast =patch.getBox().upper();

    int cols = ilast(0) - ifirst(0)+2;
    int rows = ilast(1) - ifirst(1)+2;

     int num_of_nodes = cols*rows;

    tbox::Array<bool> d_fixed_info(num_of_nodes,true);

    for(int row = 0; row < rows; row++)
    {
        for(int col = 0; col < cols; col++)
        {
            if(row == 0 || row == rows-1 || col == 0 || col == cols-1)
            {
                d_fixed_info[row*cols+col] = true;
            }
            else
                d_fixed_info[row*cols+col] = false;
        }
    }

    tbox::Pointer< pdat::NodeData<NDIM,bool> > fixed_info
            = patch.getPatchData(d_fixed_info_id);

    int count=-1;
    for(pdat::NodeIterator<NDIM> ic((*fixed_info).getBox()); ic; ic++)
    {
       (*fixed_info)(ic(),0) = d_fixed_info[++count];
    }
}

/*************************************************************************
* 步长构件: 计算并返回网格片上的稳定时间步长.  
*************************************************************************/
double MeshOpt::getPatchDt(hier::Patch<NDIM>& patch,
                           const double  time,
                           const bool    initial_time,
                           const int     flag_last_dt,
                           const double  last_dt,
                           const string& intc_name)
{

#ifdef DEBUG_CHECK_ASSERTIONS
    assert(intc_name=="TIME_STEP_SIZE");
#endif

    NULL_USE(flag_last_dt);
    NULL_USE(last_dt);

    tbox::Pointer< pdat::NodeData<NDIM,double> > coords_current =
            patch.getPatchData(d_coords_current_id);

#ifdef DEBUG_CHECK_ASSERTIONS
    assert(!coords_current.isNull());

#endif

#ifdef DEBUG_CHECK_ASSERTIONS
    assert(ghost_cells == d_zeroghosts);
#endif


    double stabdt = 0.1;


    return stabdt;

}

void EdgeMultiblockTest::initializeDataOnPatch(
   hier::Patch& patch,
   const std::shared_ptr<hier::PatchHierarchy> hierarchy,
   int level_number,
   const hier::BlockId& block_id,
   char src_or_dst)
{
   NULL_USE(hierarchy);
   NULL_USE(src_or_dst);

   if ((d_refine_option == "INTERIOR_FROM_SAME_LEVEL")
       || ((d_refine_option == "INTERIOR_FROM_COARSER_LEVEL")
           && (level_number < d_finest_level_number))) {

      for (int i = 0; i < static_cast<int>(d_variables.size()); ++i) {

         std::shared_ptr<pdat::EdgeData<double> > edge_data(
            SAMRAI_SHARED_PTR_CAST<pdat::EdgeData<double>, hier::PatchData>(
               patch.getPatchData(d_variables[i], getDataContext())));
         TBOX_ASSERT(edge_data);

         hier::Box dbox = edge_data->getGhostBox();

         edge_data->fillAll((double)block_id.getBlockValue());

      }
   }
}

/*
 *************************************************************************
 *
 * Tag cells for spherical octant problem
 *
 *************************************************************************
 */
void MblkGeometry::tagOctantCells(
   hier::Patch& patch,
   const int xyz_id,
   boost::shared_ptr<pdat::CellData<int> >& temp_tags,
   const double regrid_time,
   const int refine_tag_val)
{
   TBOX_ASSERT(d_geom_problem == "SPHERICAL_SHELL" &&
      d_sshell_type == "OCTANT");
   TBOX_ASSERT(temp_tags);

   boost::shared_ptr<pdat::NodeData<double> > xyz(
      BOOST_CAST<pdat::NodeData<double>, hier::PatchData>(
         patch.getPatchData(xyz_id)));
   TBOX_ASSERT(xyz);

   if (d_dim == tbox::Dimension(3)) {
      /*
       * Tag in X direction only
       */
      double xtag_loc_lo = d_sshell_rmin
         + (regrid_time * d_tag_velocity) - (0.5 * d_tag_width);
      double xtag_loc_hi = d_sshell_rmin
         + (regrid_time * d_tag_velocity) + (0.5 * d_tag_width);

      hier::Box pbox = patch.getBox();
      for (int k = pbox.lower(2); k <= pbox.upper(2) + 1; ++k) {
         for (int j = pbox.lower(1); j <= pbox.upper(1) + 1; ++j) {
            for (int i = pbox.lower(0); i <= pbox.upper(0) + 1; ++i) {
               hier::Index ic(i, j, k);
               pdat::NodeIndex node(ic, pdat::NodeIndex::LLL);
               hier::Index icm1(i - 1, j - 1, k - 1);
               pdat::CellIndex cell(icm1);

               double node_x_loc = (*xyz)(node, 0);

               if ((node_x_loc > xtag_loc_lo) &&
                   (node_x_loc < xtag_loc_hi)) {
                  (*temp_tags)(cell) = refine_tag_val;
               }
            }
         }
      }
   }

}

void PoissonGaussianSolution::setGridData(
   hier::Patch& patch,
   pdat::CellData<double>& exact_data,
   pdat::CellData<double>& source_data)
{
   boost::shared_ptr<geom::CartesianPatchGeometry> patch_geom(
      BOOST_CAST<geom::CartesianPatchGeometry, hier::PatchGeometry>(
         patch.getPatchGeometry()));
   TBOX_ASSERT(patch_geom);

   const double* h = patch_geom->getDx();
   const double* xl = patch_geom->getXLower();
   const int* il = &patch.getBox().lower()[0];

   {
      /* Set cell-centered data. */
      double sl[SAMRAI::MAX_DIM_VAL]; // Like XLower, except for cell.
      int j;
      for (j = 0; j < d_dim.getValue(); ++j) {
         sl[j] = xl[j] + 0.5 * h[j];
      }
      pdat::CellData<double>::iterator iter(pdat::CellGeometry::begin(patch.getBox()));
      pdat::CellData<double>::iterator iterend(pdat::CellGeometry::end(patch.getBox()));
      if (d_dim == tbox::Dimension(2)) {
         double x, y;
         for ( ; iter != iterend; ++iter) {
            const pdat::CellIndex& index = *iter;
            x = sl[0] + (index[0] - il[0]) * h[0];
            y = sl[1] + (index[1] - il[1]) * h[1];
            exact_data(index) = exactFcn(x, y);
            source_data(index) = sourceFcn(x, y);
         }
      } else if (d_dim == tbox::Dimension(3)) {
         double x, y, z;
         for ( ; iter != iterend; ++iter) {
            const pdat::CellIndex& index = *iter;
            x = sl[0] + (index[0] - il[0]) * h[0];
            y = sl[1] + (index[1] - il[1]) * h[1];
            z = sl[2] + (index[2] - il[2]) * h[2];
            exact_data(index) = exactFcn(x, y, z);
            source_data(index) = sourceFcn(x, y, z);
         }
      }
   }
}       // End patch loop.

void MblkGeometry::buildCartesianGridOnPatch(
   const hier::Patch& patch,
   const int xyz_id,
   const int level_number,
   const int block_number)
{

   boost::shared_ptr<pdat::NodeData<double> > xyz(
      BOOST_CAST<pdat::NodeData<double>, hier::PatchData>(
         patch.getPatchData(xyz_id)));

   TBOX_ASSERT(xyz);

   pdat::NodeIterator niend(pdat::NodeGeometry::end(patch.getBox()));
   for (pdat::NodeIterator ni(pdat::NodeGeometry::begin(patch.getBox()));
        ni != niend; ++ni) {
      pdat::NodeIndex node = *ni;
      if (d_block_rotation[block_number] == 0) {

         (*xyz)(node, 0) =
            d_cart_xlo[block_number][0] + node(0) * d_dx[level_number][block_number][0];
         (*xyz)(node, 1) =
            d_cart_xlo[block_number][1] + node(1) * d_dx[level_number][block_number][1];
         if (d_dim == tbox::Dimension(3)) {
            (*xyz)(node, 2) =
               d_cart_xlo[block_number][2] + node(2) * d_dx[level_number][block_number][2];
         }
      }
      if (d_block_rotation[block_number] == 1) { // I sideways, J down

         (*xyz)(node, 0) =
            d_cart_xlo[block_number][0] - node(0) * d_dx[level_number][block_number][0];
         (*xyz)(node, 1) =
            d_cart_xlo[block_number][1] + node(1) * d_dx[level_number][block_number][1];
         if (d_dim == tbox::Dimension(3)) {
            (*xyz)(node, 2) =
               d_cart_xlo[block_number][2] + node(2) * d_dx[level_number][block_number][2];
         }
      }

   }

}

boost::shared_ptr<hier::PatchData>
CellDataFactory<TYPE>::allocate(
   const hier::Patch& patch) const
{
   TBOX_ASSERT_OBJDIM_EQUALITY2(*this, patch);

   return boost::make_shared<CellData<TYPE> >(
             patch.getBox(),
             d_depth,
             d_ghosts);
}

/*************************************************************************
* 初值构件: 初始化网格结点坐标
*************************************************************************/
void MeshOpt::initializePatchData( hier::Patch<NDIM>& patch,
                                   const double  time,
                                   const bool    initial_time,
                                   const string& intc_name)
{
#ifdef DEBUG_CHECK_ASSERTIONS
    assert(intc_name=="INIT");
#endif

    (void) time;

    if (initial_time) {
        tbox::Pointer< pdat::NodeData<NDIM,double> > coords_current =
                patch.getPatchData(d_coords_current_id);

#ifdef DEBUG_CHECK_ASSERTIONS
        assert(!coords_current.isNull());
#endif

#ifdef DEBUG_CHECK_ASSERTIONS
        assert(ghost_cells == d_zeroghosts);
#endif

        // 生成初始网格.
        if(!d_grid_tool.isNull()) {
            d_grid_tool->generateDeformingMeshForDomain(patch,
                                                        d_coords_current_id);
        }

        tbox::Pointer< pdat::NodeData<NDIM,bool> > fixed_info
                = patch.getPatchData(d_fixed_info_id);

        // 赋初值
        fixed_info->fillAll(false);

        // 设定结点类型
        setNodeInfo(patch);
    }

}

void
CommTester::putCoordinatesToDatabase(
   std::shared_ptr<tbox::Database>& coords_db,
   const hier::Patch& patch)
{

   std::shared_ptr<geom::CartesianPatchGeometry> pgeom(
      SAMRAI_SHARED_PTR_CAST<geom::CartesianPatchGeometry, hier::PatchGeometry>(
         patch.getPatchGeometry()));
/*
   if (pgeom) {
      pgeom->putBlueprintCoords(coords_db, patch.getBox());
   }
*/
   const tbox::Dimension& dim(patch.getDim());

   pdat::NodeData<double> coords(patch.getBox(), dim.getValue(),
                                 hier::IntVector::getZero(dim));

   const hier::Index& box_lo = patch.getBox().lower();
   const double* x_lo = pgeom->getXLower();
   const double* dx = pgeom->getDx();
   
   pdat::NodeIterator nend = pdat::NodeGeometry::end(patch.getBox());
   for (pdat::NodeIterator itr(pdat::NodeGeometry::begin(patch.getBox())); itr != nend; ++itr) {
      const pdat::NodeIndex& ni = *itr;
      for (int d = 0; d < dim.getValue(); ++d) {
         coords(ni, d) = x_lo[d] + (ni(d)-box_lo(d))*dx[d];
      }
   }

   coords_db->putString("type", "explicit");

   std::shared_ptr<tbox::Database> values_db = coords_db->putDatabase("values");

   int data_size = coords.getArrayData().getBox().size();

   values_db->putDoubleArray("x", coords.getPointer(0), data_size);
   if (dim.getValue() > 1) {
      values_db->putDoubleArray("y", coords.getPointer(1), data_size);
   }
   if (dim.getValue() > 2) {
      values_db->putDoubleArray("z", coords.getPointer(2), data_size);
   }
}

void MeshOpt::transformMeshtoPatch(MeshImpl * mesh, hier::Patch<NDIM>& patch, MsqError& err)
{
    std::vector<Mesh::VertexHandle>  vertices;
    mesh->get_all_vertices(vertices,err);
    size_t num_of_vertices = vertices.size();
//    std::cout << num_of_vertices << std::endl;

    std::vector<MsqVertex> coords(num_of_vertices);
    mesh->vertices_get_coordinates(arrptr(vertices),arrptr(coords),num_of_vertices,err);

    tbox::Pointer< pdat::NodeData<NDIM,double> > coords_current
            = patch.getPatchData(d_coords_current_id);

    int count = 0;
    for(pdat::NodeIterator<NDIM> ic((*coords_current).getBox()); ic; ic++)
    {
        (*coords_current)(ic(),0) = coords[count][0];
        (*coords_current)(ic(),1) =  coords[count][1];
        ++count;
    }

    return;

}

/*
 *************************************************************************
 *
 * Verify results of communication operations.  This test must be
 * consistent with data initialization and boundary operations above.
 *
 *************************************************************************
 */
bool EdgeMultiblockTest::verifyResults(
   const hier::Patch& patch,
   const std::shared_ptr<hier::PatchHierarchy> hierarchy,
   const int level_number,
   const hier::BlockId& block_id)
{

   tbox::plog << "\nEntering EdgeMultiblockTest::verifyResults..." << endl;
   tbox::plog << "level_number = " << level_number << endl;
   tbox::plog << "Patch box = " << patch.getBox() << endl;

   hier::IntVector tgcw(d_dim, 0);
   for (int i = 0; i < static_cast<int>(d_variables.size()); ++i) {
      tgcw.max(patch.getPatchData(d_variables[i], getDataContext())->
         getGhostCellWidth());
   }
   hier::Box pbox = patch.getBox();

   std::shared_ptr<pdat::EdgeData<double> > solution(
      new pdat::EdgeData<double>(pbox, 1, tgcw));

   hier::Box tbox(pbox);
   tbox.grow(tgcw);

   std::shared_ptr<hier::BaseGridGeometry> grid_geom(
      hierarchy->getGridGeometry());

   hier::BoxContainer singularity(
      grid_geom->getSingularityBoxContainer(block_id));

   hier::IntVector ratio =
      hierarchy->getPatchLevel(level_number)->getRatioToLevelZero();

   singularity.refine(ratio);

   bool test_failed = false;

   for (int i = 0; i < static_cast<int>(d_variables.size()); ++i) {

      double correct = (double)block_id.getBlockValue();

      std::shared_ptr<pdat::EdgeData<double> > edge_data(
         SAMRAI_SHARED_PTR_CAST<pdat::EdgeData<double>, hier::PatchData>(
            patch.getPatchData(d_variables[i], getDataContext())));
      TBOX_ASSERT(edge_data);
      int depth = edge_data->getDepth();

      hier::Box interior_box(pbox);
      interior_box.grow(hier::IntVector(d_dim, -1));

      for (int axis = 0; axis < d_dim.getValue(); ++axis) {
         pdat::EdgeIterator ciend(pdat::EdgeGeometry::end(interior_box, axis));
         for (pdat::EdgeIterator ci(pdat::EdgeGeometry::begin(interior_box, axis));
              ci != ciend; ++ci) {
            for (int d = 0; d < depth; ++d) {
               double result = (*edge_data)(*ci, d);

               if (!tbox::MathUtilities<double>::equalEps(correct, result)) {
                  tbox::perr << "Test FAILED: ...."
                             << " : edge index = " << *ci << endl;
                  tbox::perr << "    Variable = " << d_variable_src_name[i]
                             << " : depth index = " << d << endl;
                  tbox::perr << "    result = " << result
                             << " : correct = " << correct << endl;
                  test_failed = true;
               }
            }
         }
      }

      hier::Box gbox = edge_data->getGhostBox();

      for (int axis = 0; axis < d_dim.getValue(); ++axis) {
         hier::Box patch_edge_box =
            pdat::EdgeGeometry::toEdgeBox(pbox, axis);

         hier::BoxContainer tested_neighbors;

         hier::BoxContainer sing_edge_boxlist;
         for (hier::BoxContainer::iterator si = singularity.begin();
              si != singularity.end(); ++si) {
            sing_edge_boxlist.pushFront(
               pdat::EdgeGeometry::toEdgeBox(*si, axis));
         }

         for (hier::BaseGridGeometry::ConstNeighborIterator ne(
                 grid_geom->begin(block_id));
              ne != grid_geom->end(block_id); ++ne) {

            const hier::BaseGridGeometry::Neighbor& nbr = *ne;
            if (nbr.isSingularity()) {
               continue;
//.........这里部分代码省略.........