vtkCSGGrid.h 25.6 KB
Newer Older
hrchilds's avatar
hrchilds committed
1 2
/*****************************************************************************
*
3
* Copyright (c) 2000 - 2012, Lawrence Livermore National Security, LLC
hrchilds's avatar
hrchilds committed
4
* Produced at the Lawrence Livermore National Laboratory
5
* LLNL-CODE-442911
hrchilds's avatar
hrchilds committed
6 7
* All rights reserved.
*
8
* This file is  part of VisIt. For  details, see https://visit.llnl.gov/.  The
hrchilds's avatar
hrchilds committed
9 10 11 12 13 14 15 16 17 18
* full copyright notice is contained in the file COPYRIGHT located at the root
* of the VisIt distribution or at http://www.llnl.gov/visit/copyright.html.
*
* Redistribution  and  use  in  source  and  binary  forms,  with  or  without
* modification, are permitted provided that the following conditions are met:
*
*  - Redistributions of  source code must  retain the above  copyright notice,
*    this list of conditions and the disclaimer below.
*  - Redistributions in binary form must reproduce the above copyright notice,
*    this  list of  conditions  and  the  disclaimer (as noted below)  in  the
19 20 21
*    documentation and/or other materials provided with the distribution.
*  - Neither the name of  the LLNS/LLNL nor the names of  its contributors may
*    be used to endorse or promote products derived from this software without
hrchilds's avatar
hrchilds committed
22 23 24 25 26
*    specific prior written permission.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT  HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR  IMPLIED WARRANTIES, INCLUDING,  BUT NOT  LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND  FITNESS FOR A PARTICULAR  PURPOSE
27 28 29
* ARE  DISCLAIMED. IN  NO EVENT  SHALL LAWRENCE  LIVERMORE NATIONAL  SECURITY,
* LLC, THE  U.S.  DEPARTMENT OF  ENERGY  OR  CONTRIBUTORS BE  LIABLE  FOR  ANY
* DIRECT,  INDIRECT,   INCIDENTAL,   SPECIAL,   EXEMPLARY,  OR   CONSEQUENTIAL
hrchilds's avatar
hrchilds committed
30 31 32 33 34 35 36 37 38
* DAMAGES (INCLUDING, BUT NOT  LIMITED TO, PROCUREMENT OF  SUBSTITUTE GOODS OR
* SERVICES; LOSS OF  USE, DATA, OR PROFITS; OR  BUSINESS INTERRUPTION) HOWEVER
* CAUSED  AND  ON  ANY  THEORY  OF  LIABILITY,  WHETHER  IN  CONTRACT,  STRICT
* LIABILITY, OR TORT  (INCLUDING NEGLIGENCE OR OTHERWISE)  ARISING IN ANY  WAY
* OUT OF THE  USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
* DAMAGE.
*
*****************************************************************************/

39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84
// .NAME vtkCSGGrid - Constructive Solid Geometry representation of a grid 
// .SECTION Description
// vtkCSGGrid is a derived class of vtkDataSet used to represent variables
// on a Constructive Solid Geometry "grid." The word "grid" is used loosely
// here because, of course, a CSG representation is continuous as in not
// discrete. So, the notion of a "grid" is a bit misleading. Nontheless,
// one of the important methods of this class, Discretize, produces a 
// vtkUnstructuredGrid discrete representation from the CSG
// representation.
//
// A further potential confusion here is that the base class we are
// deriving from, vtkDataSet, is pre-disposed to think in discrete
// terms. For example, it has a method to GetNumberOfPoints() which
// returns the number of discrete points in the representation.
//
// Nonetheless, since all of VTK's useful grid types are derived from
// vtkDataSet, it is still the best place to add the new CSG grid type.
// Also, since the CSG grid type will need to be discretized before it is
// used in most pipelines, that makes it even a better candidate to be
// tucked under vtkDataSet. Lastly, since VisIt's interface to get meshes 
// from database plugins expects a vtkDataSet, that pretty much clinches
// the need to define the vtkCSGGrid class in this way.
//
// There are multiple ways around vtkDataSet's discrete notions of points.
// One is to accept the analogy that in a CSG grid, the analytic functions
// play a role similar to "nodes" of an unstructured grid. For example,
// in an unstructured grid, a list of points is specified along with a 
// a list of cells. The list of cells defines how points are combined to
// form each of the cells in the grid. In a CSG grid, a list of analytic
// functions defining boundaries is specified. Then, cells are constructed
// from boolean, point-set expressions involving the boundaries. So, the
// boundaries of a CSG grid play a role similar to the nodes of an
// unstructured grid. Likewise, just as we have cell-oriented or point-
// oriented variables in an unstructured grid, we could have cell-oriented
// or boundary-oriented variables in a CSG grid. If we take this approach,
// any method in vtkDataSet's interface that refers to points can be,
// by analogy, viewed as refering to the analytic bounaries of a CSG grid.
//
// Another way around vtkDataSet's discrete notion of points is to argue
// that a CSG grid involves an infinity of points. It is the union of all
// the points satisfying all the analytic equations of the boundaries.
// What should a method like GetNumberOfPoints() return in this case?
// We can't return infinity even though that might be the correct. This
// approach has more difficulty. So, we take the former approach and treat
// the word "point" in vtkDataSet's interface as "boundary" in vtkCSGGrid's
// interface.
85 86 87 88
//
// Modifications:
//    Mark C. Miller, Tue Feb 17 17:54:34 PST 2009
//    Added operator==
89 90 91 92 93 94
//
//    Jeremy Meredith, Fri Feb 26 13:29:56 EST 2010
//    Added a new "multi-pass" approach which splits a starting data set
//    once per boundary to increase accuracy at edges/corners and improve
//    support for thin shells.
//
95 96 97 98 99 100 101 102
//    Eric Brugger, Wed Jul 25 09:58:59 PDT 2012
//    Increase the number of boundaries that can be handled by the mulit-pass
//    CSG discretization from 128 to 512.
//    Modified the multi-pass CSG discretization to perform the partitions
//    against all the boundaries and then create a vtkDataSet at the end
//    rather than creating a new vtkDataSet after partitioning with each
//    boundary.
//
103 104 105 106 107 108 109 110

// .SECTION See Also
// vtkImplicitFunction, vtkQuadric, vtkUnstructuredGrid, vtkDataSet
#ifndef __vtkCSGGrid_h
#define __vtkCSGGrid_h
#include <visit_vtk_exports.h>

#include <map>
hrchilds's avatar
hrchilds committed
111
#include <vector>
112 113 114 115 116 117 118 119

#include "vtkDataArray.h"
#include "vtkIdTypeArray.h"
#include "vtkDataSet.h"
#include "vtkImplicitFunctionCollection.h"
#include "vtkPlanes.h"
#include "vtkStructuredData.h"

120 121
#include <FixedLengthBitField.h>

hrchilds's avatar
hrchilds committed
122 123
#include <float.h>

124
class vtkPolyData;
hrchilds's avatar
hrchilds committed
125
class vtkUnstructuredGrid;
126 127 128 129 130 131 132 133

#define VTK_CSG_GRID 20

class VISIT_VTK_API vtkCSGGrid : public vtkDataSet
{
public:
  static vtkCSGGrid *New();

134
  vtkTypeMacro(vtkCSGGrid,vtkDataSet);
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
  void PrintSelf(ostream& os, vtkIndent indent);

  // Description:
  // Create a similar type object.
  vtkDataObject *MakeObject() {return vtkCSGGrid::New();};

  // Description:
  // Return what type of dataset this is.
  int GetDataObjectType() {return VTK_CSG_GRID;};
  virtual const char *GetClassName() {return "vtkCSGGrid";};

  // Description:
  // Copy the geometric and topological structure of an input rectilinear grid
  // object.
  void CopyStructure(vtkDataSet *ds);

  // Description:
  // Restore object to initial state. Release memory back to system.
  void Initialize();

  // Description:
  // Standard vtkDataSet API methods. See vtkDataSet for more information.
  vtkIdType GetNumberOfCells();
  vtkIdType GetNumberOfPoints();
  vtkIdType GetNumberOfBoundaries() const;
hrchilds's avatar
hrchilds committed
160 161 162
  double *GetPoint(vtkIdType ptId);
  double *GetBoundary(vtkIdType bndId) const;
  void GetPoint(vtkIdType id, double x[3]);
163 164 165

  vtkCell *GetCell(vtkIdType cellId);
  void GetCell(vtkIdType cellId, vtkGenericCell *cell);
hrchilds's avatar
hrchilds committed
166 167 168 169 170 171 172 173 174 175 176
  void GetCellBounds(vtkIdType cellId, double bounds[6]);
  int FindPoint(double x, double y, double z) { return this->vtkDataSet::FindPoint(x, y, z);};
  vtkIdType FindPoint(double x[3]);
  vtkIdType FindCell(double x[3], vtkCell *cell, vtkIdType cellId, double tol2,
                     int& subId, double pcoords[3], double *weights);
  vtkIdType FindCell(double x[3], vtkCell *cell, vtkGenericCell *gencell,
                     vtkIdType cellId, double tol2, int& subId, 
                     double pcoords[3], double *weights);
  vtkCell *FindAndGetCell(double x[3], vtkCell *cell, vtkIdType cellId, 
                          double tol2, int& subId, double pcoords[3],
                          double *weights);
177 178 179 180
  int GetCellType(vtkIdType cellId);
  void GetCellPoints(vtkIdType cellId, vtkIdList *ptIds);
  void GetPointCells(vtkIdType ptId, vtkIdList *cellIds);
  void ComputeBounds();
hrchilds's avatar
hrchilds committed
181 182 183
  void SetBounds(double minX, double maxX,
                 double minY, double maxY,
                 double minZ, double maxZ)
hrchilds's avatar
hrchilds committed
184 185 186
      {Bounds[0] = minX; Bounds[1] = maxX;
       Bounds[2] = minY; Bounds[3] = maxY;
       Bounds[4] = minZ; Bounds[5] = maxZ;};
187 188 189 190
  int GetMaxCellSize();
  void GetCellNeighbors(vtkIdType cellId, vtkIdList *ptIds,
                        vtkIdList *cellIds);

hrchilds's avatar
hrchilds committed
191
  void BuildVTKImplicitFunction(int zoneId, vtkImplicitFunction **func) const;
192 193 194 195

  //
  // A discretize method that returns the surfaces only
  //
hrchilds's avatar
hrchilds committed
196 197 198 199 200 201 202 203 204
  vtkPolyData  *DiscretizeSurfaces(int specificZone = -1, double tol = 0.01,
                                   double minX = -10.0, double maxX = 10.0,
                                   double minY = -10.0, double maxY = 10.0,
                                   double minZ = -10.0, double maxZ = 10.0);

  //
  // A discretize method that returns the volumetric mesh, uniformally
  // sampled to a specific number of samples in x, y and z
  //
205 206
  vtkUnstructuredGrid *GetMultiPassDiscretization(int specificZone = -1);

207 208
  bool                 DiscretizeSpaceMultiPass(const double bnds[6],
                                   const int dims[3], const int subRegion[6]);
209

hrchilds's avatar
hrchilds committed
210 211 212 213 214 215 216 217 218 219
  vtkUnstructuredGrid *DiscretizeSpace(int specificZone = -1, double tol = 0.01,
                                   double minX = -10.0, double maxX = 10.0,
                                   double minY = -10.0, double maxY = 10.0,
                                   double minZ = -10.0, double maxZ = 10.0);

  //
  // A discretize method that returns the entire spatial bounding
  // box, meshed adaptively, though discontinuously to a specified
  // tolerance (smallest edge length)
  //
hrchilds's avatar
hrchilds committed
220 221
  vtkUnstructuredGrid *DiscretizeSpace3(int specificZone = -1, int rank=0, int nprocs=1,
                                       double discTol = 0.01, double flatTol = 0.01,
hrchilds's avatar
hrchilds committed
222 223 224
                                       double minX = -10.0, double maxX = 10.0,
                                       double minY = -10.0, double maxY = 10.0,
                                       double minZ = -10.0, double maxZ = 10.0);
225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293

  // Description:
  // Return the actual size of the data in kilobytes. This number
  // is valid only after the pipeline has updated. The memory size
  // returned is guaranteed to be greater than or equal to the
  // memory required to represent the data (e.g., extra space in
  // arrays, etc. are not included in the return value). THIS METHOD
  // IS THREAD SAFE.
  unsigned long GetActualMemorySize();

  // Description:
  // Shallow and Deep copy.
  void ShallowCopy(vtkDataObject *src);  
  void DeepCopy(vtkDataObject *src);

  typedef enum {
    QUADRIC_G,
    SPHERE_G,
    SPHERE_PR,
    PLANE_GEN,
    PLANE_X,
    PLANE_Y,
    PLANE_Z,
    PLANE_PN,
    PLANE_3PT,
    PLANE_PTLINE,
    PLANE_2LINE,
    CYLINDER_GEN,
    CYLINDER_PNLR,
    CYLINDER_2PT,
    BOX
  } BoundaryType;

  typedef enum {
    INNER, // replace '=' with '<' in boundary equation
    OUTER, // replace '=' with '>' in boundary equation
    ON,    // use '=' in boundary equation
    UNION,
    INTERSECT,
    DIFF,
    XFORM
  } RegionOp;

  //
  // Methods unique to CSG Description
  //

  // Add an analytic boundary
  vtkIdType AddBoundary(BoundaryType type, int numcoeffs, const double *coeffs);

  // Get a boundary description
  void GetBoundary(vtkIdType id, int *type, int *numcoeffs, double **coeffs) const;

  // Construct and add a region from a boundary
  vtkIdType AddRegion(vtkIdType bndId, RegionOp op);

  // Construct and add a region as a binary op of other regions
  vtkIdType AddRegion(vtkIdType regIdLeft, vtkIdType regIdRight, RegionOp op);

  // Construct and add a region as transform of another region 
  vtkIdType AddRegion(vtkIdType regId, const double *xform);

  // Get a region description
  void GetRegion(vtkIdType id, vtkIdType *id1, vtkIdType *id2,
                 RegionOp *op, double **xform) const;

  // Add a cell (really just a complete region expression)
  vtkIdType AddCell(vtkIdType regId);
  vtkIdType GetCellRegionId(vtkIdType cellId) const;
hrchilds's avatar
hrchilds committed
294 295 296 297 298 299 300 301 302 303 304 305 306

  //
  // Silo convenience functions
  //
  void AddBoundaries(int nbounds, const int *const typeflags, int lcoeffs,
                const double *const coeffs);
  void AddBoundaries(int nbounds, const int *const typeflags, int lcoeffs,
                const float *const coeffs);
  void AddRegions(int nregions, const int *const lids, const int *const rids,
                  const int *const typeflags,
                  int lxforms, const double *const xforms);
  void AddZones(int nzones, const int *const zoneIds);

307 308
  bool operator==(const vtkCSGGrid &) const;

309 310 311 312
protected:
  vtkCSGGrid();
  ~vtkCSGGrid();

313
    bool EvaluateRegionBits(int region, FixedLengthBitField<64> &bits);
314

315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334
  //
  // We put this in the protected part of the interface because
  // we want to force clients of the class to use the AddBroundary/
  // GetBoundary and AddRegion/GetRegion methods.
  //
  vtkSetObjectMacro(Boundaries,vtkImplicitFunctionCollection);
  vtkGetObjectMacro(Boundaries,vtkImplicitFunctionCollection);

  vtkSetObjectMacro(Regions,vtkImplicitFunctionCollection);
  vtkGetObjectMacro(Regions,vtkImplicitFunctionCollection);

  vtkSetObjectMacro(CellRegionIds,vtkIdTypeArray);
  vtkGetObjectMacro(CellRegionIds,vtkIdTypeArray);

  vtkImplicitFunctionCollection *Boundaries; // Just the equations, f(x,y) or f(x,y,z) = 0
                                             // In other words, just the analytic boundaries.
  vtkImplicitFunctionCollection *Regions;    // Regions created by picking sides ('<' or '>') of
                                             // boundaries and binary expressions of other regions 
  vtkIdTypeArray *CellRegionIds;             // Indices into Regions of the "completed" regions

335 336 337
  // These are storage of the binary partition tree unstructured grid
  // and bitfield for the boundary tags for the multipass algorithm.
  vtkUnstructuredGrid *multipassProcessedGrid;
338
  std::vector<FixedLengthBitField<64> > *multipassTags;
339 340


hrchilds's avatar
hrchilds committed
341 342 343 344 345 346 347

  int numBoundaries;
  double *gridBoundaries;
  int numRegions;
  int *leftIds, *rightIds, *regTypeFlags;
  int numZones;
  int *gridZones;
348
private:
hrchilds's avatar
hrchilds committed
349 350


hrchilds's avatar
hrchilds committed
351 352 353 354 355 356 357
// ****************************************************************************
//  Modifications:
//
//    Hank Childs, Fri Jun  9 12:54:36 PDT 2006
//    Re-order arguments to constructor to match declaration order (addresses
//    compiler warning).
//
hrchilds's avatar
hrchilds committed
358 359 360
//    Mark C. Miller, Thu Mar 22 19:09:43 PST 2007
//    Made AddCutZones return bool indicating if it indeed added anything.
//
hrchilds's avatar
hrchilds committed
361 362
// ****************************************************************************

hrchilds's avatar
hrchilds committed
363 364 365 366 367 368 369 370 371 372 373 374
class Box
{
public:
    typedef enum {
        LT_ZERO = -1,
        EQ_ZERO =  0,
        GT_ZERO = +1
    } FuncState;

    Box(double x, double X,
        double y, double Y,
        double z, double Z,
375
        const std::vector<int>& _zids,
hrchilds's avatar
hrchilds committed
376 377
        double g000, double g001, double g010, double g011,
        double g100, double g101, double g110, double g111)
hrchilds's avatar
hrchilds committed
378
        : x0(x),y0(y),z0(z),x1(X),y1(Y),z1(Z),
hrchilds's avatar
hrchilds committed
379
          f000(g000), f001(g001), f010(g010), f011(g011),
hrchilds's avatar
hrchilds committed
380
          f100(g100), f101(g101), f110(g110), f111(g111), zids(_zids) {};
hrchilds's avatar
hrchilds committed
381

hrchilds's avatar
hrchilds committed
382 383 384
    FuncState EvalBoxStateOfBoundary(const double *const a, double tol) const;

    bool IsFlatEnough2(const double *const a, int bndId, double tol);
385
    bool CanBeCut2(const double *const a, std::map<int,int>, double tol);
hrchilds's avatar
hrchilds committed
386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429

    static FuncState ValState2(double val)
        { return val > 0.0 ? GT_ZERO :  LT_ZERO; };
    static FuncState ValState3(double val)
        { return val > 0.0 ? GT_ZERO : val < 0.0 ? LT_ZERO : EQ_ZERO; };

    static bool SameState2(FuncState s1, FuncState s2)
    {
        if (s1 == LT_ZERO || s1 == EQ_ZERO)
        {
            if (s2 == LT_ZERO || s2 == EQ_ZERO)
                return true;
            return false;
        }
        else
        {
            if (s2 == GT_ZERO || s2 == EQ_ZERO)
                return true;
            return false;
        }
    }
    static bool SameState3(FuncState s1, FuncState s2)
        { return s1 == s2; } ;

    double Resolution() const
    {
        double xres = x1 - x0; 
        double yres = y1 - y0; 
        double zres = z1 - z0; 

        if (xres < yres)
        {
            if (xres < zres)
                return xres;
            return zres;
        }
        else
        {
            if (yres < zres)
                return yres;
            return zres;
        }
    };

430
    std::vector<Box*> Subdivide() const
hrchilds's avatar
hrchilds committed
431
    {
432
        std::vector<Box*> retval;
hrchilds's avatar
hrchilds committed
433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472
        double halfx = (x0 + x1) / 2.0;
        double halfy = (y0 + y1) / 2.0;
        double halfz = (z0 + z1) / 2.0;

        Box* box000 = new Box(x0, halfx, y0, halfy, z0, halfz, zids,
                                 f000, DBL_MAX, DBL_MAX, DBL_MAX,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box001 = new Box(halfx, x1, y0, halfy, z0, halfz, zids,
                              DBL_MAX,    f001, DBL_MAX, DBL_MAX,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box010 = new Box(x0, halfx, halfy, y1, z0, halfz, zids,
                              DBL_MAX, DBL_MAX,    f010, DBL_MAX,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box011 = new Box(halfx, x1, halfy, y1, z0, halfz, zids,
                              DBL_MAX, DBL_MAX, DBL_MAX,    f011,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box100 = new Box(x0, halfx, y0, halfy, halfz, z1, zids,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX,
                                 f100, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box101 = new Box(halfx, x1, y0, halfy, halfz, z1, zids,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX,
                              DBL_MAX,    f101, DBL_MAX, DBL_MAX);
        Box* box110 = new Box(x0, halfx, halfy, y1, halfz, z1, zids,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX,
                              DBL_MAX, DBL_MAX,    f110, DBL_MAX);
        Box* box111 = new Box(halfx, x1, halfy, y1, halfz, z1, zids,
                              DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX,
                              DBL_MAX, DBL_MAX, DBL_MAX,    f111);

        retval.push_back(box000);
        retval.push_back(box001);
        retval.push_back(box010);
        retval.push_back(box011);
        retval.push_back(box100);
        retval.push_back(box101);
        retval.push_back(box110);
        retval.push_back(box111);

        return retval;
    }
473
    std::vector<Box*> SubdivideX() const
hrchilds's avatar
hrchilds committed
474
    {
475
        std::vector<Box*> retval;
hrchilds's avatar
hrchilds committed
476 477 478 479 480 481 482 483 484 485 486 487 488 489
        double halfx = (x0 + x1) / 2.0;

        Box* box0 = new Box(x0, halfx, y0, y1, z0, z1, zids,
                            f000, DBL_MAX, f010, DBL_MAX,
                            f100, DBL_MAX, f110, DBL_MAX);
        Box* box1 = new Box(halfx, x1, y0, y1, z0, z1, zids,
                            DBL_MAX, f001, DBL_MAX, f011,
                            DBL_MAX, f101, DBL_MAX, f111);

        retval.push_back(box0);
        retval.push_back(box1);

        return retval;
    }
490
    std::vector<Box*> SubdivideY() const
hrchilds's avatar
hrchilds committed
491
    {
492
        std::vector<Box*> retval;
hrchilds's avatar
hrchilds committed
493 494 495 496 497 498 499 500 501 502 503 504 505 506
        double halfy = (y0 + y1) / 2.0;

        Box* box0 = new Box(x0, x1, y0, halfy, z0, z1, zids,
                            f000, f001, DBL_MAX, DBL_MAX,
                            f100, f101, DBL_MAX, DBL_MAX);
        Box* box1 = new Box(x0, x1, halfy, y1, z0, z1, zids,
                            DBL_MAX, DBL_MAX, f010, f011,
                            DBL_MAX, DBL_MAX, f110, f111);

        retval.push_back(box0);
        retval.push_back(box1);

        return retval;
    }
507
    std::vector<Box*> SubdivideZ() const
hrchilds's avatar
hrchilds committed
508
    {
509
        std::vector<Box*> retval;
hrchilds's avatar
hrchilds committed
510 511 512 513 514 515 516 517 518 519 520 521 522 523 524
        double halfz = (z0 + z1) / 2.0;

        Box* box0 = new Box(x0, x1, y0, y1, z0, halfz, zids,
                               f000,    f001,    f010,    f011,
                            DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX);
        Box* box1 = new Box(x0, x1, y0, y1, halfz, z1, zids,
                            DBL_MAX, DBL_MAX, DBL_MAX, DBL_MAX,
                               f100,    f101,    f110,    f111);

        retval.push_back(box0);
        retval.push_back(box1);

        return retval;
    }

hrchilds's avatar
hrchilds committed
525
    double *GetPoint(int i)
hrchilds's avatar
hrchilds committed
526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593
    {
        tmp[0] = (i & 01) ? x1 : x0;
        tmp[1] = (i & 02) ? y1 : y0;
        tmp[2] = (i & 04) ? z1 : z0;
        return tmp;
    }

    double GetFunctionValue(vtkImplicitFunction *func, int i)
    {
        switch (i)
        {
            case 0:
            {
                if (f000 != DBL_MAX)
                    return f000;
                f000 = func->FunctionValue(GetPoint(i));
                return f000;
            }
            case 1:
            {
                if (f001 != DBL_MAX)
                    return f001;
                f001 = func->FunctionValue(GetPoint(i));
                return f001;
            }
            case 2:
            {
                if (f010 != DBL_MAX)
                    return f010;
                f010 = func->FunctionValue(GetPoint(i));
                return f010;
            }
            case 3:
            {
                if (f011 != DBL_MAX)
                    return f011;
                f011 = func->FunctionValue(GetPoint(i));
                return f011;
            }
            case 4:
            {
                if (f100 != DBL_MAX)
                    return f100;
                f100 = func->FunctionValue(GetPoint(i));
                return f100;
            }
            case 5:
            {
                if (f101 != DBL_MAX)
                    return f101;
                f101 = func->FunctionValue(GetPoint(i));
                return f101;
            }
            case 6:
            {
                if (f110 != DBL_MAX)
                    return f110;
                f110 = func->FunctionValue(GetPoint(i));
                return f110;
            }
            case 7:
            {
                if (f111 != DBL_MAX)
                    return f111;
                f111 = func->FunctionValue(GetPoint(i));
                return f111;
            }
        }
hrchilds's avatar
hrchilds committed
594
        return 0.0;
hrchilds's avatar
hrchilds committed
595 596 597 598
    }

    double x0,y0,z0,x1,y1,z1;
    double f000,f001,f010,f011,f100,f101,f110,f111;
hrchilds's avatar
hrchilds committed
599
    double tmp[3];
600
    std::vector<int> zids;
hrchilds's avatar
hrchilds committed
601 602
};

603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650
class coord_t {
public:
    float c[3];
    coord_t() {c[0]=0; c[1]=0; c[2]=0;};
    coord_t(const float c_[3])
    {c[0]=c_[0]; c[1]=c_[1]; c[2]=c_[2];};
    coord_t(float a0, float a1, float a2)
    {c[0]=a0; c[1]=a1; c[2]=a2; };
    coord_t& operator=(const coord_t& rhs)
    { c[0]=rhs.c[0]; c[1]=rhs.c[1]; c[2]=rhs.c[2]; return *this;};
};

struct coordcomp {
    bool operator() (const coord_t& lhs, const coord_t& rhs) const
    {
        if (lhs.c[0] < rhs.c[0])
        {
            return true;
        }
        else if (lhs.c[0] == rhs.c[0])
        {
    if (lhs.c[1] < rhs.c[1])
    {
        return true;
    }
    else if (lhs.c[1] == rhs.c[1])
    {
        if (lhs.c[2] < rhs.c[2])
        {
            return true;
        }
        else 
        {
            return false;
        }
    }
    else
    {
        return false;
    }
}
else
{
    return false;
}
}
};

651
typedef std::map<coord_t,int,coordcomp> coordmap_t;
652 653 654 655 656 657 658 659

static bool AddCutZones(vtkUnstructuredGrid *cutBox, vtkPoints *points,
                   vtkUnstructuredGrid *ugrid,
                   coordmap_t& nodemap);
static void MakeMeshZone(const Box *aBox, vtkPoints *points,
                   vtkUnstructuredGrid *ugrid,
                   coordmap_t& nodemap);
bool MakeMeshZonesByCuttingBox4(const Box *aBox,
660 661
                   const std::map<int,int>& boundaryToStateMap,
                   std::map<int,int>& boundaryToSenseMap, int zoneId,
662 663 664
                   vtkPoints *points, vtkUnstructuredGrid *ugrid,
                   coordmap_t& nodemap);
bool MakeMeshZonesByCuttingBox2(const Box *aBox,
665 666
                   const std::map<int,int>& boundaryToStateMap,
                   std::map<int,int>& boundaryToSenseMap, int zoneId,
667 668 669
                   vtkPoints *points, vtkUnstructuredGrid *ugrid,
                   coordmap_t& nodemap);
static void MakeMeshZonesByCuttingBox(const Box *aBox,
670 671
                   std::map<vtkImplicitFunction*,Box::FuncState> funcToStateMap,
                   std::vector<RegionOp> senses,
672 673
                   vtkPoints *points, vtkUnstructuredGrid *ugrid,
                   coordmap_t& nodemap);
674
void AddBoundariesForZone2(int, std::vector<int> *bnds, std::vector<int> *senses);
675
void AddBoundariesForZone(vtkImplicitFunction *func,
676 677
                           std::vector<vtkImplicitFunction*> *bnds,
                           std::vector<RegionOp> *senses);
678
int EvalBoxStateOfRegion(const Box *const curBox, int regId,
679
std::map<int,int>& boundaryToStateMap, double tol);
680 681 682 683 684 685

double tmpFloats[32];                       // temporary storage to help satisfy interface
                                     //    requirements of vtkDataSet

vtkPlanes *Universe;                       // The "universe" set (a maximally sized box)

686
std::map<vtkImplicitFunction *, vtkIdType> funcMap;
687 688 689 690 691 692

vtkImplicitFunction *GetBoundaryFunc(vtkIdType id) const;
vtkImplicitFunction *GetRegionFunc(vtkIdType id) const;

vtkCSGGrid(const vtkCSGGrid&);             // Not implemented.
void operator=(const vtkCSGGrid&);         // Not implemented.
693 694 695 696 697
};


inline vtkIdType vtkCSGGrid::GetNumberOfPoints()
{
698 699 700
vtkErrorMacro("GetNumberOfPoints() means GetNumberOfBoundaries()");
vtkErrorMacro("Use GetNumberOfBoundaries() to avoid this message");
return GetNumberOfBoundaries();
701 702 703 704
};

inline vtkIdType vtkCSGGrid::GetNumberOfBoundaries() const
{
705
return (vtkIdType) this->Boundaries->GetNumberOfItems();
706 707 708 709 710 711 712 713
};

inline vtkIdType vtkCSGGrid::GetNumberOfCells() 
{
  return CellRegionIds->GetNumberOfTuples();
}

#endif