MOAB: Mesh Oriented datABase  (version 5.3.0)
ScdInterface.hpp
Go to the documentation of this file.
00001 /** \file ScdInterface.hpp
00002  */
00003 #ifndef SCD_INTERFACE
00004 #define SCD_INTERFACE
00005 
00006 #include "moab/Interface.hpp"
00007 #include "moab/HomXform.hpp"
00008 
00009 #include <iostream>
00010 #include <vector>
00011 #include <cassert>
00012 
00013 #include "moab/win32_config.h"
00014 
00015 namespace moab
00016 {
00017 
00018 class StructuredElementSeq;
00019 class EntitySequence;
00020 class ScdVertexData;
00021 class EntitySequence;
00022 class ScdBox;
00023 class ParallelComm;
00024 
00025 /** \class ScdInterface ScdInterface.hpp "moab/ScdInterface.hpp"
00026  * \brief A structured mesh interface for MOAB-based data
00027  *
00028  * Structured mesh in MOAB is created and accessed through the ScdInterface and ScdBox classes.
00029  *
00030  * \section Construction Construction and Representation
00031  * Structured mesh can be constructed in one of two ways.  First, a rectangular block of mesh,
00032  * both vertices and edges/quads/hexes, can be created in one shot, using the construct_box method.
00033  * In this case, there are single sequences of vertices/entities.  The second method for creating
00034  * structured mesh is to create the structured blocks of vertices and elements separately.  In
00035  * this case, different blocks of elements can share blocks of vertices, and each block of
00036  * elements has its own independent parametric space.  The algorithms behind this representation
00037  * are described in T. Tautges, "MOAB-SD: Integrated structured and unstructured mesh
00038  * representation", Eng. w Comp, vol 20 no. 3.
00039  *
00040  * Structured mesh is represented in MOAB down at the element sequence level, which is something
00041  * applications don't see.  In addition, when structured mesh is created, entity sets are also
00042  * created and tagged with information about the parametric space.  In particular, the BOX_DIMS
00043  * tag is used to indicate the lower and upper corners in parametric space (this
00044  * tag is integer size 6).  Structured mesh blocks are also available through ScdBox class objects
00045  * returned by ScdInterface.  These class objects should be treated only as references to the
00046  * structured mesh blocks; that is, the structured mesh referenced by these objects is not deleted
00047  * when the ScdBox instance is destroyed.  Functions for destroying the actual mesh are available
00048  * on this class, though.
00049  *
00050  * Structured mesh blocks are returned in the form of ScdBox class objects.  Each ScdBox instance
00051  * represents a rectangular block of vertices and possibly elements (edges, quads, or hexes).  The
00052  * edge/quad/hex entity handles for a ScdBox are guaranteed to be contiguous, starting at a starting
00053  * value which is also available through the ScdBox class.  However, vertex handles may or may not
00054  * be contiguous, depending on the construction method.  The start vertex handle is also available
00055  * from the ScdBox class.
00056  *
00057  * \section Parameters Parametric Space
00058  *
00059  * Each structured box has a parametric (ijk) space, which can be queried through the ScdBox
00060  * interface. For non-periodic boxes, the edge/quad/hex parameter bounds are one less in each
00061  * dimension than that of the vertices, otherwise they are the same as the vertex parameter bounds.
00062  * In a parallel representation, boxes are locally non-periodic by default, but global ids are
00063  * assigned such that the last set of vertices in a periodic direction match those of the first set
00064  * of vertices in that direction.
00065  *
00066  * Entity handles are allocated with the i parameter varying fastest, then j, then k.
00067  *
00068  * \section Per Periodic Meshes
00069  * Boxes can be periodic in i, or j, or both i and j.  If only i or j is periodic, the corresponding
00070  * mesh is a strip or an annular cylinder; if both i and j are periodic, the corresponding mesh is
00071  * an annular torus.  A box cannot be periodic in all three parameters.  If i and/or j is periodic,
00072  * and assuming IMIN/JMIN is zero, the parameter extents in the/each periodic direction (IMAX/JMAX)
00073  * for vertices and edges/faces/hexes are the same, and the vertices on the "top" end in the
00074  * periodic direction are at parameter value IMIN/JMIN.
00075  *
00076  * \section Par Parallel Representation
00077  *
00078  * For parallel structured meshes, each local mesh (the mesh on a given process) looks like a
00079  * non-periodic structured mesh, and there are both local and global parameters of the structured
00080  * mesh.  If the mesh is periodic in a given direction, the last process in the periodic direction
00081  * has local IMAX/JMAX that is one greater than the global IMAX/JMAX.
00082  *
00083  *
00084  * In parallel, the periodicity described in the previous paragraph is "local periodicity"; there is
00085  * also the notion of global periodicity.  For serial meshes, those concepts are the same.  In
00086  * parallel, a mesh can be locally non-periodic but globally periodic in a given direction.  In that
00087  * case, the local mesh is still non-periodic, i.e. the parametric extents for edges is one fewer
00088  * than that of vertices in that direction. However, vertices are given global ids such that they
00089  * match those of the parametric minimum in that direction. Geometric positions of the vertices at
00090  * the high end should still be greater than the ones just below.
00091  *
00092  * \section Adjs Adjacent Entities
00093  * This interface supports parametric access to intermediate-dimension entities, e.g. adjacent faces
00094  * and edges in a 3d mesh.  In this case, a direction parameter is added, to identify the parametric
00095  * direction of the entities being requested.  For example, to retrieve the faces adjacent to a hex
00096  * with parameters ijk, in the i parametric direction, you would use the parameters ijk0.  These
00097  * intermediate entities are not stored in a structured representation, but their parametric
00098  * positions can be evaluated based on their adjacencies to higher-dimensional entities.  Thanks to
00099  * Milad Fatenejad for the thinking behind this.
00100  *
00101  * \section Evaluation Evaluation
00102  * The ScdBox class provides functions for evaluating the mesh based on the ijk parameter space.
00103  * These functions are inlined where possible, for efficiency.
00104  */
00105 
00106 //! struct for keeping parallel data in one place
00107 class ScdParData
00108 {
00109   public:
00110     ScdParData() : partMethod( NOPART ), pComm( NULL )
00111     {
00112         gDims[0] = gDims[1] = gDims[2] = gDims[3] = gDims[4] = gDims[5] = 0;
00113         gPeriodic[0] = gPeriodic[1] = gPeriodic[2] = 0;
00114         pDims[0] = pDims[1] = pDims[2] = 0;
00115     }
00116 
00117     //! Partition method enumeration; these strategies are described in comments for
00118     //! compute_partition_alljorkori, compute_partition_alljkbal, compute_partition_sqij,
00119     //! compute_partition_sqjk, and compute_partition_sqijk
00120     enum PartitionMethod
00121     {
00122         ALLJORKORI = 0,
00123         ALLJKBAL,
00124         SQIJ,
00125         SQJK,
00126         SQIJK,
00127         TRIVIAL,
00128         RCBZOLTAN,
00129         NOPART
00130     };
00131 
00132     //! Partition method names
00133     static MOAB_EXPORT const char* PartitionMethodNames[NOPART + 1];
00134 
00135     //! partition method used to partition global parametric space
00136     int partMethod;
00137 
00138     //! lower and upper corners of global box
00139     int gDims[6];
00140 
00141     //! is globally periodic in i or j or k
00142     int gPeriodic[3];
00143 
00144     //! number of procs in each direction
00145     int pDims[3];
00146 
00147     //! parallel communicator object for this par scd mesh
00148     ParallelComm* pComm;
00149 };
00150 
00151 class MOAB_EXPORT ScdInterface
00152 {
00153   public:
00154     friend class ScdBox;
00155 
00156     //! Constructor
00157     /** Constructor; if find_boxes is true, this will search for entity sets marked as
00158      * structured blocks, based on the BOX_DIMS tag.  Structured mesh blocks will be stored
00159      * in this interface class for future retrieval.  Structured mesh blocks created through
00160      * this interface will also be stored here.
00161      * \param impl MOAB instance
00162      * \param find_boxes If true, search all the entity sets, caching the structured mesh blocks
00163      */
00164     ScdInterface( Interface* impl, bool find_boxes = false );
00165 
00166     // Destructor
00167     ~ScdInterface();
00168 
00169     //! Return the MOAB Interface instance *
00170     Interface* impl() const;
00171 
00172     //! Construct new structured mesh box, including both vertices and elements
00173     /** Parameter range
00174      * for vertex box is [low-high], for elements is [low-high).  Construct quads by passing
00175      * in low[2] == high[2], and edges by passing in low[1] == high[1] and low[2] == high[2].
00176      * The result is passed back in a ScdBox*, which is a *reference* to the box of structured mesh.
00177      * That is, the actual mesh is retained in MOAB when the ScdBox is destroyed.  To actually
00178      * destroy the mesh, call the destroy_mesh function on the ScdBox object first, before
00179      * destroying it. \param low Lower corner in parameter space \param high Higher corner in
00180      * parameter space \param coords Coordinates of vertices, interleaved (xyzxyz...); if NULL,
00181      * coords are set to parametric values \param num_coords Number of coordinate values \param
00182      * new_box Reference to box of structured mesh \param lperiodic[3] If lperiodic[s] != 0,
00183      * direction s is locally periodic \param par_data If non-NULL, this will get stored on the
00184      * ScdBox once created, contains info about global parallel nature of ScdBox across procs \param
00185      * assign_global_ids If true, assigns 1-based global ids to vertices using GLOBAL_ID_TAG_NAME
00186      * \param resolve_shared_ents If != -1, resolves shared entities up to and including dimension
00187      * equal to value
00188      */
00189     ErrorCode construct_box( HomCoord low, HomCoord high, const double* const coords, unsigned int num_coords,
00190                              ScdBox*& new_box, int* const lperiodic = NULL, ScdParData* const par_data = NULL,
00191                              bool assign_global_ids = false, int resolve_shared_ents = -1 );
00192 
00193     //! Create a structured sequence of vertices, quads, or hexes
00194     /** Starting handle for the sequence is available from the returned ScdBox.
00195      * If creating a structured quad or hex box, subsequent calls must be made to ScdBox::add_vbox,
00196      * until all the vertices have been filled in for the box.
00197      * \param low Lower corner of structured box
00198      * \param high Higher corner of structured box
00199      * \param type EntityType, one of MBVERTEX, MBEDGE, MBQUAD, MBHEX
00200      * \param starting_id Requested start id of entities
00201      * \param new_box Reference to the newly created box of entities
00202      * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
00203      * (s=[0,1,2])
00204      */
00205     ErrorCode create_scd_sequence( const HomCoord& low, const HomCoord& high, EntityType type, int starting_id,
00206                                    ScdBox*& new_box, int* is_periodic = NULL );
00207 
00208     //! Return all the structured mesh blocks in this MOAB instance, as ScdBox objects
00209     /** Return the structured blocks in this MOAB instance.  If these were not searched for
00210      * at instantiation time, then the search is done now.
00211      * \param boxes Vector of ScdBox objects representing structured mesh blocks
00212      */
00213     ErrorCode find_boxes( std::vector< ScdBox* >& boxes );
00214 
00215     //! Return all the structured mesh blocks in this MOAB instance, as entity set handles
00216     /** Return the structured blocks in this MOAB instance.  If these were not searched for
00217      * at instantiation time, then the search is done now.
00218      * \param boxes Range of entity set objects representing structured mesh blocks
00219      */
00220     ErrorCode find_boxes( Range& boxes );
00221 
00222     //! Return all the structured mesh blocks known by ScdInterface (does not search)
00223     /** Return the structured blocks in this ScdInterface instance.  Does not search for new boxes,
00224      * just returns the contents of the list.
00225      * \param boxes Structured boxes
00226      */
00227     ErrorCode get_boxes( std::vector< ScdBox* >& boxes );
00228 
00229     //! Return the tag marking the lower and upper corners of boxes
00230     /**
00231      * \param create_if_missing If the tag does not yet exist, create it
00232      */
00233     Tag box_dims_tag( bool create_if_missing = true );
00234 
00235     //! Return the tag marking the global lower and upper corners of boxes
00236     /**
00237      * \param create_if_missing If the tag does not yet exist, create it
00238      */
00239     Tag global_box_dims_tag( bool create_if_missing = true );
00240 
00241     //! Return the tag marking the partitioning method used to partition the box in parallel
00242     /**
00243      * \param create_if_missing If the tag does not yet exist, create it
00244      */
00245     Tag part_method_tag( bool create_if_missing = true );
00246 
00247     //! Return the tag marking whether box is periodic in i and j
00248     /**
00249      * \param create_if_missing If the tag does not yet exist, create it
00250      */
00251     Tag box_periodic_tag( bool create_if_missing = true );
00252 
00253     //! Return the tag marking the ScdBox for a set
00254     /**
00255      * \param create_if_missing If the tag does not yet exist, create it
00256      */
00257     Tag box_set_tag( bool create_if_missing = true );
00258 
00259     //! Return the ScdBox corresponding to the entity set passed in
00260     /** If the entity isn't a structured box set, NULL is returned.
00261      * \param eh Entity whose box is being queried
00262      */
00263     ScdBox* get_scd_box( EntityHandle eh );
00264 
00265     //! Compute a partition of structured parameter space
00266     /** Compute a partition of structured parameter space, based on data in the ScdParData
00267      * passed in.  Results are passed back in arguments, which application can set back into
00268      * par_data argument if they so desire.
00269      * \param np Number of processors
00270      * \param nr Rank of this processor
00271      * \param par_data ScdParData object that contains input global parameter space, desired
00272      *           partitioning method, and information about global periodicity.
00273      * \param ldims Local parameters for grid
00274      * \param lperiodic Whether or not a given dimension is locally periodic
00275      * \param pdims Number of procs in i, j, k directions
00276      */
00277     static ErrorCode compute_partition( int np, int nr, const ScdParData& par_data, int* ldims, int* lperiodic = NULL,
00278                                         int* pdims = NULL );
00279 
00280     //! Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
00281     //! all 3 params
00282     /** Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
00283      * all 3 params \param np (in) Total # procs \param nr Processor from which neighbor is
00284      * requested \param spd (in) ScdParData containing part method, gdims, gperiodic data \param
00285      * dijk(*) (in) Direction being queried, = +/-1 or 0 \param pto (out) Processor holding the
00286      * destination part \param rdims(6) (out) Parametric min/max of destination part \param
00287      * facedims(6) (out) Parametric min/max of interface between pfrom and pto; if at the max in a
00288      * periodic direction, set to global min of that direction \param across_bdy(3) (out) If
00289      * across_bdy[i] is -1(1), interface with pto is across periodic lower(upper) bdy in parameter
00290      * i, 0 otherwise
00291      */
00292     static ErrorCode get_neighbor( int np, int nr, const ScdParData& spd, const int* const dijk, int& pto, int* rdims,
00293                                    int* facedims, int* across_bdy );
00294 
00295     //! Tag vertices with sharing data for parallel representations
00296     /** Given the ParallelComm object to use, tag the vertices shared with other processors
00297      */
00298     ErrorCode tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth );
00299 
00300     //! Tag vertices with sharing data for parallel representations
00301     /** Given the ParallelComm object to use, tag the vertices shared with other processors
00302      */
00303     ErrorCode tag_shared_vertices( ParallelComm* pcomm, ScdBox* box );
00304 
00305   protected:
00306     //! Remove the box from the list on ScdInterface
00307     ErrorCode remove_box( ScdBox* box );
00308 
00309     //! Add the box to the list on ScdInterface
00310     ErrorCode add_box( ScdBox* box );
00311 
00312   private:
00313     //! Create an entity set for a box, and tag with the parameters
00314     /** \param low Lower corner parameters for this box
00315      * \param high Upper corner parameters for this box
00316      * \param scd_set Entity set created
00317      * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
00318      * (s=[0,1,2])
00319      */
00320     ErrorCode create_box_set( const HomCoord& low, const HomCoord& high, EntityHandle& scd_set,
00321                               int* is_periodic = NULL );
00322 
00323     //! Compute a partition of structured parameter space
00324     /** Partitions the structured parametric space by partitioning j, k, or i only.
00325      * If j is greater than #procs, partition that, else k, else i.
00326      * For description of arguments, see ScdInterface::compute_partition.
00327      */
00328     inline static ErrorCode compute_partition_alljorkori( int np, int nr, const int gijk[6], const int* const gperiodic,
00329                                                           int* lijk, int* lperiodic, int* pijk );
00330 
00331     //! Compute a partition of structured parameter space
00332     /** Partitions the structured parametric space by partitioning j, and possibly k,
00333      * seeking square regions of jk space
00334      * For description of arguments, see ScdInterface::compute_partition.
00335      */
00336     inline static ErrorCode compute_partition_alljkbal( int np, int nr, const int gijk[6], const int* const gperiodic,
00337                                                         int* lijk, int* lperiodic, int* pijk );
00338 
00339     //! Compute a partition of structured parameter space
00340     /** Partitions the structured parametric space by seeking square ij partitions
00341      * For description of arguments, see ScdInterface::compute_partition.
00342      */
00343     inline static ErrorCode compute_partition_sqij( int np, int nr, const int gijk[6], const int* const gperiodic,
00344                                                     int* lijk, int* lperiodic, int* pijk );
00345 
00346     //! Compute a partition of structured parameter space
00347     /** Partitions the structured parametric space by seeking square jk partitions
00348      * For description of arguments, see ScdInterface::compute_partition.
00349      */
00350     inline static ErrorCode compute_partition_sqjk( int np, int nr, const int gijk[6], const int* const gperiodic,
00351                                                     int* lijk, int* lperiodic, int* pijk );
00352 
00353     //! Compute a partition of structured parameter space
00354     /** Partitions the structured parametric space by seeking square ijk partitions
00355      * For description of arguments, see ScdInterface::compute_partition.
00356      */
00357     inline static ErrorCode compute_partition_sqijk( int np, int nr, const int gijk[6], const int* const gperiodic,
00358                                                      int* lijk, int* lperiodic, int* pijk );
00359 
00360     //! Get vertices shared with other processors
00361     /** Shared vertices returned as indices into each proc's handle space
00362      * \param box Box used to get parametric space info
00363      * \param procs Procs this proc shares vertices with
00364      * \param offsets Offsets into indices list for each proc
00365      * \param shared_indices local/remote indices of shared vertices
00366      */
00367     static ErrorCode get_shared_vertices( ParallelComm* pcomm, ScdBox* box, std::vector< int >& procs,
00368                                           std::vector< int >& offsets, std::vector< int >& shared_indices );
00369 
00370     static ErrorCode get_indices( const int* const ldims, const int* const rdims, const int* const across_bdy,
00371                                   int* face_dims, std::vector< int >& shared_indices );
00372 
00373     static ErrorCode get_neighbor_alljorkori( int np, int pfrom, const int* const gdims, const int* const gperiodic,
00374                                               const int* const dijk, int& pto, int* rdims, int* facedims,
00375                                               int* across_bdy );
00376 
00377     static ErrorCode get_neighbor_alljkbal( int np, int pfrom, const int* const gdims, const int* const gperiodic,
00378                                             const int* const dijk, int& pto, int* rdims, int* facedims,
00379                                             int* across_bdy );
00380 
00381     static ErrorCode get_neighbor_sqij( int np, int pfrom, const int* const gdims, const int* const gperiodic,
00382                                         const int* const dijk, int& pto, int* rdims, int* facedims, int* across_bdy );
00383 
00384     static ErrorCode get_neighbor_sqjk( int np, int pfrom, const int* const gdims, const int* const gperiodic,
00385                                         const int* const dijk, int& pto, int* rdims, int* facedims, int* across_bdy );
00386 
00387     static ErrorCode get_neighbor_sqijk( int np, int pfrom, const int* const gdims, const int* const gperiodic,
00388                                          const int* const dijk, int& pto, int* rdims, int* facedims, int* across_bdy );
00389 
00390     static int gtol( const int* gijk, int i, int j, int k );
00391 
00392     //! assign global ids to vertices in this box
00393     ErrorCode assign_global_ids( ScdBox* box );
00394 
00395     //! interface instance
00396     Interface* mbImpl;
00397 
00398     //! whether we've searched the database for boxes yet
00399     bool searchedBoxes;
00400 
00401     //! structured mesh blocks; stored as ScdBox objects, can get sets from those
00402     std::vector< ScdBox* > scdBoxes;
00403 
00404     //! tag representing whether box is periodic in i and j
00405     Tag boxPeriodicTag;
00406 
00407     //! tag representing box lower and upper corners
00408     Tag boxDimsTag;
00409 
00410     //! tag representing global lower and upper corners
00411     Tag globalBoxDimsTag;
00412 
00413     //! tag representing partition method
00414     Tag partMethodTag;
00415 
00416     //! tag pointing from set to ScdBox
00417     Tag boxSetTag;
00418 };
00419 
00420 class MOAB_EXPORT ScdBox
00421 {
00422     friend class ScdInterface;
00423 
00424   public:
00425     //! Destructor
00426     ~ScdBox();
00427 
00428     //! Return the ScdInterface responsible for this box
00429     inline ScdInterface* sc_impl() const;
00430 
00431     //! Add a vertex box to this box
00432     /* Add a vertex box to the element sequence referenced by this box.  The passed in vbox must
00433      * be a vertex box, with parametric extents no larger than that of this box.  This vbox is
00434      * oriented to this box by matching parameters from1-from3 in vbox to to1-to3 in this box.
00435      * If bb_input is true, only the part of the vertex sequence between bb_min and bb_max is
00436      * referenced \param vbox The vertex box being added to this box \param from1 1st reference
00437      * point on vbox \param to1 1st reference point on this box \param from2 2nd reference point on
00438      * vbox \param to2 2nd reference point on this box \param from3 3rd reference point on vbox
00439      * \param to3 3rd reference point on this box
00440      * \param bb_input If true, subsequent parameters list extents of vbox referenced
00441      * \param bb_min Lower corner of rectangle referenced
00442      * \param bb_max Upper corner of rectangle referenced
00443      */
00444     ErrorCode add_vbox( ScdBox* vbox, HomCoord from1, HomCoord to1, HomCoord from2, HomCoord to2, HomCoord from3,
00445                         HomCoord to3, bool bb_input = false, const HomCoord& bb_min = HomCoord::getUnitv( 0 ),
00446                         const HomCoord& bb_max = HomCoord::getUnitv( 0 ) );
00447     // const HomCoord &bb_min = HomCoord::unitv[0],
00448     // const HomCoord &bb_max = HomCoord::unitv[0]);
00449 
00450     //! Return whether this box has all its vertices defined
00451     /** Tests whether vertex boxs added with add_vbox have completely defined the vertex parametric
00452      * space for this box.
00453      *
00454      */
00455     bool boundary_complete() const;
00456 
00457     //! Return highest topological dimension of box
00458     inline int box_dimension() const;
00459 
00460     //! Starting vertex handle for this box
00461     inline EntityHandle start_vertex() const;
00462 
00463     //! Starting entity handle for this box
00464     /** If this is a vertex box, the start vertex handle is returned.
00465      */
00466     inline EntityHandle start_element() const;
00467 
00468     //! Return the number of elements in the box
00469     /* Number of elements is (boxSize[0]-1)(boxSize[1]-1)(boxSize[2]-1)
00470      */
00471     inline int num_elements() const;
00472 
00473     //! Return the number of vertices in the box
00474     /* Number of vertices is boxSize[0] * boxSize[1] * boxSize[2]
00475      */
00476     inline int num_vertices() const;
00477 
00478     //! Return the parametric coordinates for this box
00479     /**
00480      * \return IJK parameters of lower and upper corners
00481      */
00482     inline const int* box_dims() const;
00483 
00484     //! Return the lower corner parametric coordinates for this box
00485     inline HomCoord box_min() const;
00486 
00487     //! Return the upper corner parametric coordinates for this box
00488     inline HomCoord box_max() const;
00489 
00490     //! Return the parameter extents for this box
00491     inline HomCoord box_size() const;
00492 
00493     //! Return the parametric extents for this box
00494     /**
00495      * \param ijk IJK extents of this box
00496      */
00497     inline void box_size( int* ijk ) const;
00498 
00499     //! Return the parametric extents for this box
00500     /**
00501      * \param i I extent of this box
00502      * \param j J extent of this box
00503      * \param k K extent of this box
00504      */
00505     void box_size( int& i, int& j, int& k ) const;
00506 
00507     //! Get the element at the specified coordinates
00508     /**
00509      * \param ijk Parametric coordinates being evaluated
00510      */
00511     EntityHandle get_element( const HomCoord& ijk ) const;
00512 
00513     //! Get the element at the specified coordinates
00514     /**
00515      * \param i Parametric coordinates being evaluated
00516      * \param j Parametric coordinates being evaluated
00517      * \param k Parametric coordinates being evaluated
00518      */
00519     EntityHandle get_element( int i, int j = 0, int k = 0 ) const;
00520 
00521     //! Get the vertex at the specified coordinates
00522     /**
00523      * \param ijk Parametric coordinates being evaluated
00524      */
00525     EntityHandle get_vertex( const HomCoord& ijk ) const;
00526 
00527     //! Get the vertex at the specified coordinates
00528     /**
00529      * \param i Parametric coordinates being evaluated
00530      * \param j Parametric coordinates being evaluated
00531      * \param k Parametric coordinates being evaluated
00532      */
00533     EntityHandle get_vertex( int i, int j = 0, int k = 0 ) const;
00534 
00535     //! Get parametric coordinates of the specified entity
00536     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00537      * in this ScdBox.
00538      * \param ent Entity being queried
00539      * \param i Parametric coordinates returned
00540      * \param j Parametric coordinates returned
00541      * \param k Parametric coordinates returned
00542      * \param dir Parametric coordinate direction returned (in case of getting adjacent
00543      *            edges (2d, 3d) or faces (3d); not modified otherwise
00544      */
00545     ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const;
00546 
00547     //! Get parametric coordinates of the specified entity, intermediate entities not allowed (no
00548     //! dir parameter)
00549     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00550      * in this ScdBox, or MB_FAILURE if the entity is an intermediate-dimension entity.
00551      * \param ent Entity being queried
00552      * \param i Parametric coordinates returned
00553      * \param j Parametric coordinates returned
00554      * \param k Parametric coordinates returned
00555      */
00556     ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k ) const;
00557 
00558     //! Get parametric coordinates of the specified entity
00559     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00560      * in this ScdBox.
00561      * \param ent Entity being queried
00562      * \param ijkd Parametric coordinates returned (including direction, in case of
00563      *            getting adjacent edges (2d, 3d) or faces (3d))
00564      */
00565     ErrorCode get_params( EntityHandle ent, HomCoord& ijkd ) const;
00566 
00567     /** \brief Get the adjacent edge or face at a parametric location
00568      * This function gets the left (i=0), front (j=0), or bottom (k=0) edge or face for a parametric
00569      * element. Left, front, or bottom is indicated by dir = 0, 1, or 2, resp.  All edges and faces
00570      * in a structured mesh block can be accessed using these parameters. \param dim Dimension of
00571      * adjacent entity being requested \param i Parametric coordinates of cell being evaluated
00572      * \param j Parametric coordinates of cell being evaluated
00573      * \param k Parametric coordinates of cell being evaluated
00574      * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d)
00575      * \param ent Entity returned from this function
00576      * \param create_if_missing If true, creates the entity if it doesn't already exist
00577      */
00578     ErrorCode get_adj_edge_or_face( int dim, int i, int j, int k, int dir, EntityHandle& ent,
00579                                     bool create_if_missing = true ) const;
00580 
00581     //! Return whether the box contains the parameters passed in
00582     /**
00583      * \param i Parametric coordinates being evaluated
00584      * \param j Parametric coordinates being evaluated
00585      * \param k Parametric coordinates being evaluated
00586      */
00587     bool contains( int i, int j, int k ) const;
00588 
00589     //! Return whether the box contains the parameters passed in
00590     /**
00591      * \param i Parametric coordinates being evaluated
00592      * \param j Parametric coordinates being evaluated
00593      * \param k Parametric coordinates being evaluated
00594      */
00595     bool contains( const HomCoord& ijk ) const;
00596 
00597     //! Set/Get the entity set representing the box
00598     void box_set( EntityHandle this_set );
00599     EntityHandle box_set();
00600 
00601     //! Get coordinate arrays for vertex coordinates for a structured block
00602     /** Returns error if there isn't a single vertex sequence associated with this structured block
00603      * \param xc X coordinate array pointer returned
00604      * \param yc Y coordinate array pointer returned
00605      * \param zc Z coordinate array pointer returned
00606      */
00607     ErrorCode get_coordinate_arrays( double*& xc, double*& yc, double*& zc );
00608 
00609     //! Get read-only coordinate arrays for vertex coordinates for a structured block
00610     /** Returns error if there isn't a single vertex sequence associated with this structured block
00611      * \param xc X coordinate array pointer returned
00612      * \param yc Y coordinate array pointer returned
00613      * \param zc Z coordinate array pointer returned
00614      */
00615     ErrorCode get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const;
00616 
00617     //! Return whether box is locally periodic in i
00618     /** Return whether box is locally periodic in i
00619      * \return True if box is locally periodic in i direction
00620      */
00621     bool locally_periodic_i() const;
00622 
00623     //! Return whether box is locally periodic in j
00624     /** Return whether box is locally periodic in j
00625      * \return True if box is locally periodic in j direction
00626      */
00627     bool locally_periodic_j() const;
00628 
00629     //! Return whether box is locally periodic in k
00630     /** Return whether box is locally periodic in k
00631      * \return True if box is locally periodic in k direction
00632      */
00633     bool locally_periodic_k() const;
00634 
00635     //! Set local periodicity
00636     /**
00637      * \param lperiodic Vector of ijk periodicities to set this box to
00638      */
00639     void locally_periodic( bool lperiodic[3] );
00640 
00641     //! Get local periodicity
00642     /**
00643      * \return Vector of ijk periodicities for this box
00644      */
00645     const int* locally_periodic() const;
00646 
00647     //! Return parallel data
00648     /** Return parallel data, if there is any
00649      * \return par_data Parallel data set on this box
00650      */
00651     ScdParData& par_data()
00652     {
00653         return parData;
00654     }
00655 
00656     //! Return parallel data
00657     /** Return parallel data, if there is any
00658      * \return par_data Parallel data set on this box
00659      */
00660     const ScdParData& par_data() const
00661     {
00662         return parData;
00663     }
00664 
00665     //! set parallel data
00666     /** Set parallel data for this box
00667      * \param par_data Parallel data to be set on this box
00668      */
00669     void par_data( const ScdParData& par_datap )
00670     {
00671         parData = par_datap;
00672     }
00673 
00674   private:
00675     //! Constructor
00676     /** Create a structured box instance; this constructor is private because it should only be
00677      * called from ScdInterface, a friend class.  This constructor takes two sequences, one of which
00678      * can be NULL.  If both sequences come in non-NULL, the first should be a VertexSequence*
00679      * corresponding to a structured vertex sequence and the second should be a
00680      * StructuredElementSeq*.  If the 2nd is NULL, the first can be either of those types.  The
00681      * other members of this class are taken from the sequences (e.g. parametric space) or the box
00682      * set argument.  Tags on the box set should be set from the caller. \param sc_impl A
00683      * ScdInterface instance \param box_set Entity set representing this rectangle \param seq1 An
00684      * EntitySequence (see ScdBox description) \param seq2 An EntitySequence (see ScdBox
00685      * description), or NULL
00686      */
00687     ScdBox( ScdInterface* sc_impl, EntityHandle box_set, EntitySequence* seq1, EntitySequence* seq2 = NULL );
00688 
00689     //! function to get vertex handle directly from sequence
00690     /** \param i Parameter being queried
00691      * \param j Parameter being queried
00692      * \param k Parameter being queried
00693      */
00694     EntityHandle get_vertex_from_seq( int i, int j, int k ) const;
00695 
00696     //! set the vertex sequence
00697     ErrorCode vert_dat( ScdVertexData* vert_dat );
00698 
00699     //! get the vertex sequence
00700     ScdVertexData* vert_dat() const;
00701 
00702     //! set the element sequence
00703     ErrorCode elem_seq( EntitySequence* elem_seq );
00704 
00705     //! get the element sequence
00706     StructuredElementSeq* elem_seq() const;
00707 
00708     //! Set the starting vertex handle for this box
00709     void start_vertex( EntityHandle startv );
00710 
00711     //! Set the starting entity handle for this box
00712     void start_element( EntityHandle starte );
00713 
00714     //! interface instance
00715     ScdInterface* scImpl;
00716 
00717     //! entity set representing this box
00718     EntityHandle boxSet;
00719 
00720     //! vertex sequence this box represents, if there's only one, otherwise they're
00721     //! retrieved from the element sequence
00722     ScdVertexData* vertDat;
00723 
00724     //! element sequence this box represents
00725     StructuredElementSeq* elemSeq;
00726 
00727     //! starting vertex handle for this box
00728     EntityHandle startVertex;
00729 
00730     //! starting element handle for this box
00731     EntityHandle startElem;
00732 
00733     //! lower and upper corners
00734     int boxDims[6];
00735 
00736     //! is locally periodic in i or j or k
00737     int locallyPeriodic[3];
00738 
00739     //! parallel data associated with this box, if any
00740     ScdParData parData;
00741 
00742     //! parameter extents
00743     HomCoord boxSize;
00744 
00745     //! convenience parameters, (boxSize[1]-1)*(boxSize[0]-1) and boxSize[0]-1
00746     int boxSizeIJ;
00747     int boxSizeIJM1;
00748     int boxSizeIM1;
00749 };
00750 
00751 inline ErrorCode ScdInterface::compute_partition( int np, int nr, const ScdParData& par_data, int* ldims,
00752                                                   int* lperiodic, int* pdims )
00753 {
00754     ErrorCode rval = MB_SUCCESS;
00755     switch( par_data.partMethod )
00756     {
00757         case ScdParData::ALLJORKORI:
00758         case -1:
00759             rval = compute_partition_alljorkori( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00760             break;
00761         case ScdParData::ALLJKBAL:
00762             rval = compute_partition_alljkbal( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00763             break;
00764         case ScdParData::SQIJ:
00765             rval = compute_partition_sqij( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00766             break;
00767         case ScdParData::SQJK:
00768             rval = compute_partition_sqjk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00769             break;
00770         case ScdParData::SQIJK:
00771             rval = compute_partition_sqijk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00772             break;
00773         default:
00774             rval = MB_FAILURE;
00775             break;
00776     }
00777 
00778     return rval;
00779 }
00780 
00781 inline ErrorCode ScdInterface::compute_partition_alljorkori( int np, int nr, const int gijk[6],
00782                                                              const int* const gperiodic, int* ldims, int* lperiodic,
00783                                                              int* pijk )
00784 {
00785     // partition *the elements* over the parametric space; 1d partition for now, in the j, k, or i
00786     // parameters
00787     int tmp_lp[3], tmp_pijk[3];
00788     if( !lperiodic ) lperiodic = tmp_lp;
00789     if( !pijk ) pijk = tmp_pijk;
00790 
00791     for( int i = 0; i < 3; i++ )
00792         lperiodic[i] = gperiodic[i];
00793 
00794     if( np == 1 )
00795     {
00796         if( ldims )
00797         {
00798             ldims[0] = gijk[0];
00799             ldims[3] = gijk[3];
00800             ldims[1] = gijk[1];
00801             ldims[4] = gijk[4];
00802             ldims[2] = gijk[2];
00803             ldims[5] = gijk[5];
00804         }
00805         pijk[0] = pijk[1] = pijk[2] = 1;
00806     }
00807     else
00808     {
00809         if( gijk[4] - gijk[1] > np )
00810         {
00811             // partition j over procs
00812             int dj    = ( gijk[4] - gijk[1] ) / np;
00813             int extra = ( gijk[4] - gijk[1] ) % np;
00814             ldims[1]  = gijk[1] + nr * dj + std::min( nr, extra );
00815             ldims[4]  = ldims[1] + dj + ( nr < extra ? 1 : 0 );
00816 
00817             if( gperiodic[1] && np > 1 )
00818             {
00819                 lperiodic[1] = 0;
00820                 ldims[4]++;
00821             }
00822 
00823             ldims[2] = gijk[2];
00824             ldims[5] = gijk[5];
00825             ldims[0] = gijk[0];
00826             ldims[3] = gijk[3];
00827             pijk[0] = pijk[2] = 1;
00828             pijk[1]           = np;
00829         }
00830         else if( gijk[5] - gijk[2] > np )
00831         {
00832             // partition k over procs
00833             int dk    = ( gijk[5] - gijk[2] ) / np;
00834             int extra = ( gijk[5] - gijk[2] ) % np;
00835             ldims[2]  = gijk[2] + nr * dk + std::min( nr, extra );
00836             ldims[5]  = ldims[2] + dk + ( nr < extra ? 1 : 0 );
00837 
00838             ldims[1] = gijk[1];
00839             ldims[4] = gijk[4];
00840             ldims[0] = gijk[0];
00841             ldims[3] = gijk[3];
00842             pijk[0] = pijk[1] = 1;
00843             pijk[2]           = np;
00844         }
00845         else if( gijk[3] - gijk[0] > np )
00846         {
00847             // partition i over procs
00848             int di    = ( gijk[3] - gijk[0] ) / np;
00849             int extra = ( gijk[3] - gijk[0] ) % np;
00850             ldims[0]  = gijk[0] + nr * di + std::min( nr, extra );
00851             ldims[3]  = ldims[0] + di + ( nr < extra ? 1 : 0 );
00852 
00853             if( gperiodic[0] && np > 1 )
00854             {
00855                 lperiodic[0] = 0;
00856                 ldims[3]++;
00857             }
00858 
00859             ldims[2] = gijk[2];
00860             ldims[5] = gijk[5];
00861             ldims[1] = gijk[1];
00862             ldims[4] = gijk[4];
00863 
00864             pijk[1] = pijk[2] = 1;
00865             pijk[0]           = np;
00866         }
00867         else
00868         {
00869             // Couldn't find a suitable partition...
00870             return MB_FAILURE;
00871         }
00872     }
00873 
00874     return MB_SUCCESS;
00875 }
00876 
00877 inline ErrorCode ScdInterface::compute_partition_alljkbal( int np, int nr, const int gijk[6],
00878                                                            const int* const gperiodic, int* ldims, int* lperiodic,
00879                                                            int* pijk )
00880 {
00881     int tmp_lp[3], tmp_pijk[3];
00882     if( !lperiodic ) lperiodic = tmp_lp;
00883     if( !pijk ) pijk = tmp_pijk;
00884 
00885     for( int i = 0; i < 3; i++ )
00886         lperiodic[i] = gperiodic[i];
00887 
00888     if( np == 1 )
00889     {
00890         if( ldims )
00891         {
00892             ldims[0] = gijk[0];
00893             ldims[3] = gijk[3];
00894             ldims[1] = gijk[1];
00895             ldims[4] = gijk[4];
00896             ldims[2] = gijk[2];
00897             ldims[5] = gijk[5];
00898         }
00899         pijk[0] = pijk[1] = pijk[2] = 1;
00900     }
00901     else
00902     {
00903         // improved, possibly 2-d partition
00904         std::vector< double > kfactors;
00905         kfactors.push_back( 1 );
00906         int K = gijk[5] - gijk[2];
00907         for( int i = 2; i < K; i++ )
00908             if( !( K % i ) && !( np % i ) ) kfactors.push_back( i );
00909         kfactors.push_back( K );
00910 
00911         // compute the ideal nj and nk
00912         int J          = gijk[4] - gijk[1];
00913         double njideal = sqrt( ( (double)( np * J ) ) / ( (double)K ) );
00914         double nkideal = ( njideal * K ) / J;
00915 
00916         int nk, nj;
00917         if( nkideal < 1.0 )
00918         {
00919             nk = 1;
00920             nj = np;
00921         }
00922         else
00923         {
00924             std::vector< double >::iterator vit = std::lower_bound( kfactors.begin(), kfactors.end(), nkideal );
00925             if( vit == kfactors.begin() )
00926                 nk = 1;
00927             else
00928                 nk = (int)*( --vit );
00929             nj = np / nk;
00930         }
00931 
00932         int dk = K / nk;
00933         int dj = J / nj;
00934 
00935         ldims[2] = gijk[2] + ( nr % nk ) * dk;
00936         ldims[5] = ldims[2] + dk;
00937 
00938         int extra = J % nj;
00939 
00940         ldims[1] = gijk[1] + ( nr / nk ) * dj + std::min( nr / nk, extra );
00941         ldims[4] = ldims[1] + dj + ( nr / nk < extra ? 1 : 0 );
00942 
00943         ldims[0] = gijk[0];
00944         ldims[3] = gijk[3];
00945 
00946         if( gperiodic[1] && np > 1 )
00947         {
00948             lperiodic[1] = 0;
00949             if( nr / nk == nj - 1 ) { ldims[1]++; }
00950         }
00951 
00952         pijk[0] = 1;
00953         pijk[1] = nj;
00954         pijk[2] = nk;
00955     }
00956 
00957     return MB_SUCCESS;
00958 }
00959 
00960 inline ErrorCode ScdInterface::compute_partition_sqij( int np, int nr, const int gijk[6], const int* const gperiodic,
00961                                                        int* ldims, int* lperiodic, int* pijk )
00962 {
00963     int tmp_lp[3], tmp_pijk[3];
00964     if( !lperiodic ) lperiodic = tmp_lp;
00965     if( !pijk ) pijk = tmp_pijk;
00966 
00967     // square IxJ partition
00968 
00969     for( int i = 0; i < 3; i++ )
00970         lperiodic[i] = gperiodic[i];
00971 
00972     if( np == 1 )
00973     {
00974         if( ldims )
00975         {
00976             ldims[0] = gijk[0];
00977             ldims[3] = gijk[3];
00978             ldims[1] = gijk[1];
00979             ldims[4] = gijk[4];
00980             ldims[2] = gijk[2];
00981             ldims[5] = gijk[5];
00982         }
00983         pijk[0] = pijk[1] = pijk[2] = 1;
00984     }
00985     else
00986     {
00987         std::vector< double > pfactors, ppfactors;
00988         for( int i = 2; i <= np / 2; i++ )
00989             if( !( np % i ) )
00990             {
00991                 pfactors.push_back( i );
00992                 ppfactors.push_back( ( (double)( i * i ) ) / np );
00993             }
00994         pfactors.push_back( np );
00995         ppfactors.push_back( (double)np );
00996 
00997         // ideally, Px/Py = I/J
00998         double ijratio = ( (double)( gijk[3] - gijk[0] ) ) / ( (double)( gijk[4] - gijk[1] ) );
00999 
01000         unsigned int ind                        = 0;
01001         std::vector< double >::iterator optimal = std::lower_bound( ppfactors.begin(), ppfactors.end(), ijratio );
01002         if( optimal == ppfactors.end() ) { ind = ppfactors.size() - 1; }
01003         else
01004         {
01005             ind = optimal - ppfactors.begin();
01006             if( ind && fabs( ppfactors[ind - 1] - ijratio ) < fabs( ppfactors[ind] - ijratio ) ) ind--;
01007         }
01008 
01009         // VARIABLES DESCRIBING THE MESH:
01010         // pi, pj = # procs in i and j directions
01011         // nri, nrj = my proc's position in i, j directions
01012         // I, J = # edges/elements in i, j directions
01013         // iextra, jextra = # procs having extra edge in i/j direction
01014         // top_i, top_j = if true, I'm the last proc in the i/j direction
01015         // i, j = # edges locally in i/j direction, *not* including one for iextra/jextra
01016         int pi = pfactors[ind];
01017         int pj = np / pi;
01018 
01019         int I = ( gijk[3] - gijk[0] ), J = ( gijk[4] - gijk[1] );
01020         int iextra = I % pi, jextra = J % pj, i = I / pi, j = J / pj;
01021         int nri = nr % pi, nrj = nr / pi;
01022 
01023         if( ldims )
01024         {
01025             ldims[0] = gijk[0] + i * nri + std::min( iextra, nri );
01026             ldims[3] = ldims[0] + i + ( nri < iextra ? 1 : 0 );
01027             ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
01028             ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
01029 
01030             ldims[2] = gijk[2];
01031             ldims[5] = gijk[5];
01032 
01033             if( gperiodic[0] && pi > 1 )
01034             {
01035                 lperiodic[0] = 0;
01036                 if( nri == pi - 1 ) ldims[3]++;
01037             }
01038             if( gperiodic[1] && pj > 1 )
01039             {
01040                 lperiodic[1] = 0;
01041                 if( nrj == pj - 1 ) ldims[4]++;
01042             }
01043         }
01044 
01045         pijk[0] = pi;
01046         pijk[1] = pj;
01047         pijk[2] = 1;
01048     }
01049 
01050     return MB_SUCCESS;
01051 }
01052 
01053 inline ErrorCode ScdInterface::compute_partition_sqjk( int np, int nr, const int gijk[6], const int* const gperiodic,
01054                                                        int* ldims, int* lperiodic, int* pijk )
01055 {
01056     int tmp_lp[3], tmp_pijk[3];
01057     if( !lperiodic ) lperiodic = tmp_lp;
01058     if( !pijk ) pijk = tmp_pijk;
01059 
01060     // square JxK partition
01061     for( int i = 0; i < 3; i++ )
01062         lperiodic[i] = gperiodic[i];
01063 
01064     if( np == 1 )
01065     {
01066         if( ldims )
01067         {
01068             ldims[0] = gijk[0];
01069             ldims[3] = gijk[3];
01070             ldims[1] = gijk[1];
01071             ldims[4] = gijk[4];
01072             ldims[2] = gijk[2];
01073             ldims[5] = gijk[5];
01074         }
01075         pijk[0] = pijk[1] = pijk[2] = 1;
01076     }
01077     else
01078     {
01079         std::vector< double > pfactors, ppfactors;
01080         for( int p = 2; p <= np; p++ )
01081             if( !( np % p ) )
01082             {
01083                 pfactors.push_back( p );
01084                 ppfactors.push_back( ( (double)( p * p ) ) / np );
01085             }
01086 
01087         // ideally, Pj/Pk = J/K
01088         int pj, pk;
01089         if( gijk[5] == gijk[2] )
01090         {
01091             pk = 1;
01092             pj = np;
01093         }
01094         else
01095         {
01096             double jkratio = ( (double)( gijk[4] - gijk[1] ) ) / ( (double)( gijk[5] - gijk[2] ) );
01097 
01098             std::vector< double >::iterator vit = std::lower_bound( ppfactors.begin(), ppfactors.end(), jkratio );
01099             if( vit == ppfactors.end() )
01100                 --vit;
01101             else if( vit != ppfactors.begin() && fabs( *( vit - 1 ) - jkratio ) < fabs( ( *vit ) - jkratio ) )
01102                 --vit;
01103             int ind = vit - ppfactors.begin();
01104 
01105             pj = 1;
01106             if( ind >= 0 && !pfactors.empty() ) pfactors[ind];
01107             pk = np / pj;
01108         }
01109 
01110         int K = ( gijk[5] - gijk[2] ), J = ( gijk[4] - gijk[1] );
01111         int jextra = J % pj, kextra = K % pk, j = J / pj, k = K / pk;
01112         int nrj = nr % pj, nrk = nr / pj;
01113         ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
01114         ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
01115         ldims[2] = gijk[2] + k * nrk + std::min( kextra, nrk );
01116         ldims[5] = ldims[2] + k + ( nrk < kextra ? 1 : 0 );
01117 
01118         ldims[0] = gijk[0];
01119         ldims[3] = gijk[3];
01120 
01121         if( gperiodic[1] && pj > 1 )
01122         {
01123             lperiodic[1] = 0;
01124             if( nrj == pj - 1 ) ldims[4]++;
01125         }
01126 
01127         pijk[0] = 1;
01128         pijk[1] = pj;
01129         pijk[2] = pk;
01130     }
01131 
01132     return MB_SUCCESS;
01133 }
01134 
01135 inline ErrorCode ScdInterface::compute_partition_sqijk( int np, int nr, const int* const gijk,
01136                                                         const int* const gperiodic, int* ldims, int* lperiodic,
01137                                                         int* pijk )
01138 {
01139     if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE;
01140 
01141     int tmp_lp[3], tmp_pijk[3];
01142     if( !lperiodic ) lperiodic = tmp_lp;
01143     if( !pijk ) pijk = tmp_pijk;
01144 
01145     // square IxJxK partition
01146 
01147     for( int i = 0; i < 3; i++ )
01148         lperiodic[i] = gperiodic[i];
01149 
01150     if( np == 1 )
01151     {
01152         if( ldims )
01153             for( int i = 0; i < 6; i++ )
01154                 ldims[i] = gijk[i];
01155         pijk[0] = pijk[1] = pijk[2] = 1;
01156         return MB_SUCCESS;
01157     }
01158 
01159     std::vector< int > pfactors;
01160     pfactors.push_back( 1 );
01161     for( int i = 2; i <= np / 2; i++ )
01162         if( !( np % i ) ) pfactors.push_back( i );
01163     pfactors.push_back( np );
01164 
01165     // test for IJ, JK, IK
01166     int IJK[3], dIJK[3];
01167     for( int i = 0; i < 3; i++ )
01168         IJK[i] = std::max( gijk[3 + i] - gijk[i], 1 );
01169     // order IJK from lo to hi
01170     int lo = 0, hi = 0;
01171     for( int i = 1; i < 3; i++ )
01172     {
01173         if( IJK[i] < IJK[lo] ) lo = i;
01174         if( IJK[i] > IJK[hi] ) hi = i;
01175     }
01176     if( lo == hi ) hi = ( lo + 1 ) % 3;
01177     int mid = 3 - lo - hi;
01178     // search for perfect subdivision of np that balances #cells
01179     int perfa_best = -1, perfb_best = -1;
01180     double ratio = 0.0;
01181     for( int po = 0; po < (int)pfactors.size(); po++ )
01182     {
01183         for( int pi = po; pi < (int)pfactors.size() && np / ( pfactors[po] * pfactors[pi] ) >= pfactors[pi]; pi++ )
01184         {
01185             int p3 =
01186                 std::find( pfactors.begin(), pfactors.end(), np / ( pfactors[po] * pfactors[pi] ) ) - pfactors.begin();
01187             if( p3 == (int)pfactors.size() || pfactors[po] * pfactors[pi] * pfactors[p3] != np )
01188                 continue;  // po*pi should exactly factor np
01189             assert( po <= pi && pi <= p3 );
01190             // by definition, po <= pi <= p3
01191             double minl =
01192                        std::min( std::min( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] ),
01193                    maxl =
01194                        std::max( std::max( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] );
01195             if( minl / maxl > ratio )
01196             {
01197                 ratio      = minl / maxl;
01198                 perfa_best = po;
01199                 perfb_best = pi;
01200             }
01201         }
01202     }
01203     if( perfa_best == -1 || perfb_best == -1 ) return MB_FAILURE;
01204 
01205     // VARIABLES DESCRIBING THE MESH:
01206     // pijk[i] = # procs in direction i
01207     // numr[i] = my proc's position in direction i
01208     // dIJK[i] = # edges/elements in direction i
01209     // extra[i]= # procs having extra edge in direction i
01210     // top[i] = if true, I'm the last proc in direction i
01211 
01212     pijk[lo]     = pfactors[perfa_best];
01213     pijk[mid]    = pfactors[perfb_best];
01214     pijk[hi]     = ( np / ( pfactors[perfa_best] * pfactors[perfb_best] ) );
01215     int extra[3] = { 0, 0, 0 }, numr[3];
01216     for( int i = 0; i < 3; i++ )
01217     {
01218         dIJK[i]  = IJK[i] / pijk[i];
01219         extra[i] = IJK[i] % pijk[i];
01220     }
01221     numr[2] = nr / ( pijk[0] * pijk[1] );
01222     int rem = nr % ( pijk[0] * pijk[1] );
01223     numr[1] = rem / pijk[0];
01224     numr[0] = rem % pijk[0];
01225     for( int i = 0; i < 3; i++ )
01226     {
01227         extra[i]     = IJK[i] % dIJK[i];
01228         ldims[i]     = gijk[i] + numr[i] * dIJK[i] + std::min( extra[i], numr[i] );
01229         ldims[3 + i] = ldims[i] + dIJK[i] + ( numr[i] < extra[i] ? 1 : 0 );
01230     }
01231 
01232     return MB_SUCCESS;
01233 }
01234 
01235 inline int ScdInterface::gtol( const int* gijk, int i, int j, int k )
01236 {
01237     return ( ( k - gijk[2] ) * ( gijk[3] - gijk[0] + 1 ) * ( gijk[4] - gijk[1] + 1 ) +
01238              ( j - gijk[1] ) * ( gijk[3] - gijk[0] + 1 ) + i - gijk[0] );
01239 }
01240 
01241 inline ErrorCode ScdInterface::get_indices( const int* const ldims, const int* const rdims, const int* const across_bdy,
01242                                             int* face_dims, std::vector< int >& shared_indices )
01243 {
01244     // check for going across periodic bdy and face_dims not in my ldims (I'll always be on top in
01245     // that case)...
01246     if( across_bdy[0] > 0 && face_dims[0] != ldims[3] )
01247         face_dims[0] = face_dims[3] = ldims[3];
01248     else if( across_bdy[0] < 0 && face_dims[0] != ldims[0] )
01249         face_dims[0] = face_dims[3] = ldims[0];
01250     if( across_bdy[1] > 0 && face_dims[1] != ldims[4] )
01251         face_dims[1] = face_dims[4] = ldims[4];
01252     else if( across_bdy[1] < 0 && face_dims[1] != ldims[1] )
01253         face_dims[0] = face_dims[3] = ldims[1];
01254 
01255     for( int k = face_dims[2]; k <= face_dims[5]; k++ )
01256         for( int j = face_dims[1]; j <= face_dims[4]; j++ )
01257             for( int i = face_dims[0]; i <= face_dims[3]; i++ )
01258                 shared_indices.push_back( gtol( ldims, i, j, k ) );
01259 
01260     if( across_bdy[0] > 0 && face_dims[0] != rdims[0] )
01261         face_dims[0] = face_dims[3] = rdims[0];
01262     else if( across_bdy[0] < 0 && face_dims[0] != rdims[3] )
01263         face_dims[0] = face_dims[3] = rdims[3];
01264     if( across_bdy[1] > 0 && face_dims[1] != rdims[1] )
01265         face_dims[1] = face_dims[4] = rdims[1];
01266     else if( across_bdy[1] < 0 && face_dims[1] != rdims[4] )
01267         face_dims[0] = face_dims[3] = rdims[4];
01268 
01269     for( int k = face_dims[2]; k <= face_dims[5]; k++ )
01270         for( int j = face_dims[1]; j <= face_dims[4]; j++ )
01271             for( int i = face_dims[0]; i <= face_dims[3]; i++ )
01272                 shared_indices.push_back( gtol( rdims, i, j, k ) );
01273 
01274     return MB_SUCCESS;
01275 }
01276 
01277 inline ErrorCode ScdInterface::get_neighbor( int np, int pfrom, const ScdParData& spd, const int* const dijk, int& pto,
01278                                              int* rdims, int* facedims, int* across_bdy )
01279 {
01280     if( !dijk[0] && !dijk[1] && !dijk[2] )
01281     {
01282         // not going anywhere, return
01283         pto = -1;
01284         return MB_SUCCESS;
01285     }
01286 
01287     switch( spd.partMethod )
01288     {
01289         case ScdParData::ALLJORKORI:
01290         case -1:
01291             return get_neighbor_alljorkori( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims,
01292                                             across_bdy );
01293         case ScdParData::ALLJKBAL:
01294             return get_neighbor_alljkbal( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01295         case ScdParData::SQIJ:
01296             return get_neighbor_sqij( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01297         case ScdParData::SQJK:
01298             return get_neighbor_sqjk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01299         case ScdParData::SQIJK:
01300             return get_neighbor_sqijk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01301         default:
01302             break;
01303     }
01304 
01305     return MB_FAILURE;
01306 }
01307 
01308 inline ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth )
01309 {
01310     ScdBox* box = get_scd_box( seth );
01311     if( !box )
01312     {
01313         // look for contained boxes
01314         Range tmp_range;
01315         ErrorCode rval = mbImpl->get_entities_by_type( seth, MBENTITYSET, tmp_range );
01316         if( MB_SUCCESS != rval ) return rval;
01317         for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
01318         {
01319             box = get_scd_box( *rit );
01320             if( box ) break;
01321         }
01322     }
01323 
01324     if( !box ) return MB_FAILURE;
01325 
01326     return tag_shared_vertices( pcomm, box );
01327 }
01328 
01329 inline ScdInterface* ScdBox::sc_impl() const
01330 {
01331     return scImpl;
01332 }
01333 
01334 inline EntityHandle ScdBox::start_vertex() const
01335 {
01336     return startVertex;
01337 }
01338 
01339 inline void ScdBox::start_vertex( EntityHandle startv )
01340 {
01341     startVertex = startv;
01342 }
01343 
01344 inline EntityHandle ScdBox::start_element() const
01345 {
01346     return startElem;
01347 }
01348 
01349 inline void ScdBox::start_element( EntityHandle starte )
01350 {
01351     startElem = starte;
01352 }
01353 
01354 inline int ScdBox::num_elements() const
01355 {
01356     if( !startElem ) return 0;  // not initialized yet
01357 
01358     /* for a structured mesh, total number of elements is obtained by multiplying
01359         number of elements in each direction
01360       number of elements in each direction is given by number of vertices in that direction minus 1
01361       if periodic in that direction, the last vertex is the same as first one, count one more
01362       element
01363       */
01364     int num_e_i = ( -1 == boxSize[0] || 1 == boxSize[0] ) ? 1 : boxSize[0] - 1;
01365     if( locallyPeriodic[0] ) ++num_e_i;
01366 
01367     int num_e_j = ( -1 == boxSize[1] || 1 == boxSize[1] ) ? 1 : boxSize[1] - 1;
01368     if( locallyPeriodic[1] ) ++num_e_j;
01369 
01370     int num_e_k = ( -1 == boxSize[2] || 1 == boxSize[2] ) ? 1 : boxSize[2] - 1;
01371     if( locallyPeriodic[2] ) ++num_e_k;
01372 
01373     return num_e_i * num_e_j * num_e_k;
01374 }
01375 
01376 inline int ScdBox::num_vertices() const
01377 {
01378     return boxSize[0] * ( !boxSize[1] ? 1 : boxSize[1] ) * ( !boxSize[2] ? 1 : boxSize[2] );
01379 }
01380 
01381 inline const int* ScdBox::box_dims() const
01382 {
01383     return boxDims;
01384 }
01385 
01386 inline HomCoord ScdBox::box_min() const
01387 {
01388     return HomCoord( boxDims, 3 );
01389 }
01390 
01391 inline HomCoord ScdBox::box_max() const
01392 {
01393     return HomCoord( boxDims + 3, 3 );
01394 }
01395 
01396 inline HomCoord ScdBox::box_size() const
01397 {
01398     return boxSize;
01399 }
01400 
01401 inline void ScdBox::box_size( int* ijk ) const
01402 {
01403     ijk[0] = boxSize[0];
01404     ijk[1] = boxSize[1];
01405     ijk[2] = boxSize[2];
01406 }
01407 
01408 inline void ScdBox::box_size( int& i, int& j, int& k ) const
01409 {
01410     i = boxSize[0];
01411     j = boxSize[1];
01412     k = boxSize[2];
01413 }
01414 
01415 inline EntityHandle ScdBox::get_element( int i, int j, int k ) const
01416 {
01417     return ( !startElem
01418                  ? 0
01419                  : startElem + ( k - boxDims[2] ) * boxSizeIJM1 + ( j - boxDims[1] ) * boxSizeIM1 + i - boxDims[0] );
01420 }
01421 
01422 inline EntityHandle ScdBox::get_element( const HomCoord& ijk ) const
01423 {
01424     return get_element( ijk[0], ijk[1], ijk[2] );
01425 }
01426 
01427 inline EntityHandle ScdBox::get_vertex( int i, int j, int k ) const
01428 {
01429     return ( vertDat
01430                  ? startVertex + ( boxDims[2] == -1 && boxDims[5] == -1 ? 0 : ( k - boxDims[2] ) ) * boxSizeIJ +
01431                        ( boxDims[1] == -1 && boxDims[4] == -1 ? 0 : ( j - boxDims[1] ) ) * boxSize[0] + i - boxDims[0]
01432                  : get_vertex_from_seq( i, j, k ) );
01433 }
01434 
01435 inline EntityHandle ScdBox::get_vertex( const HomCoord& ijk ) const
01436 {
01437     return get_vertex( ijk[0], ijk[1], ijk[2] );
01438 }
01439 
01440 inline bool ScdBox::contains( const HomCoord& ijk ) const
01441 {
01442     return ( ijk >= HomCoord( boxDims, 3 ) && ijk <= HomCoord( boxDims + 3, 3 ) );
01443 }
01444 
01445 inline bool ScdBox::contains( int i, int j, int k ) const
01446 {
01447     return contains( HomCoord( i, j, k ) );
01448 }
01449 
01450 inline void ScdBox::box_set( EntityHandle this_set )
01451 {
01452     boxSet = this_set;
01453 }
01454 
01455 inline EntityHandle ScdBox::box_set()
01456 {
01457     return boxSet;
01458 }
01459 
01460 inline ScdVertexData* ScdBox::vert_dat() const
01461 {
01462     return vertDat;
01463 }
01464 
01465 inline StructuredElementSeq* ScdBox::elem_seq() const
01466 {
01467     return elemSeq;
01468 }
01469 
01470 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const
01471 {
01472     HomCoord hc;
01473     ErrorCode rval = get_params( ent, hc );
01474     if( MB_SUCCESS == rval )
01475     {
01476         i   = hc[0];
01477         j   = hc[1];
01478         k   = hc[2];
01479         dir = hc[3];
01480     }
01481 
01482     return rval;
01483 }
01484 
01485 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k ) const
01486 {
01487     HomCoord hc;
01488     ErrorCode rval = get_params( ent, hc );
01489     if( MB_SUCCESS == rval )
01490     {
01491         i = hc[0];
01492         j = hc[1];
01493         k = hc[2];
01494     }
01495 
01496     return rval;
01497 }
01498 
01499 inline bool ScdBox::locally_periodic_i() const
01500 {
01501     return locallyPeriodic[0];
01502 }
01503 
01504 inline bool ScdBox::locally_periodic_j() const
01505 {
01506     return locallyPeriodic[1];
01507 }
01508 
01509 inline bool ScdBox::locally_periodic_k() const
01510 {
01511     return locallyPeriodic[2];
01512 }
01513 
01514 inline void ScdBox::locally_periodic( bool lperiodic[3] )
01515 {
01516     for( int i = 0; i < 3; i++ )
01517         locallyPeriodic[i] = lperiodic[i];
01518 }
01519 
01520 inline const int* ScdBox::locally_periodic() const
01521 {
01522     return locallyPeriodic;
01523 }
01524 
01525 std::ostream& operator<<( std::ostream& str, const ScdParData& pd );
01526 
01527 }  // namespace moab
01528 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines