MOAB: Mesh Oriented datABase  (version 5.4.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,
00190                              HomCoord high,
00191                              const double* const coords,
00192                              unsigned int num_coords,
00193                              ScdBox*& new_box,
00194                              int* const lperiodic       = NULL,
00195                              ScdParData* const par_data = NULL,
00196                              bool assign_global_ids     = false,
00197                              int resolve_shared_ents    = -1 );
00198 
00199     //! Create a structured sequence of vertices, quads, or hexes
00200     /** Starting handle for the sequence is available from the returned ScdBox.
00201      * If creating a structured quad or hex box, subsequent calls must be made to ScdBox::add_vbox,
00202      * until all the vertices have been filled in for the box.
00203      * \param low Lower corner of structured box
00204      * \param high Higher corner of structured box
00205      * \param type EntityType, one of MBVERTEX, MBEDGE, MBQUAD, MBHEX
00206      * \param starting_id Requested start id of entities
00207      * \param new_box Reference to the newly created box of entities
00208      * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
00209      * (s=[0,1,2])
00210      */
00211     ErrorCode create_scd_sequence( const HomCoord& low,
00212                                    const HomCoord& high,
00213                                    EntityType type,
00214                                    int starting_id,
00215                                    ScdBox*& new_box,
00216                                    int* is_periodic = NULL );
00217 
00218     //! Return all the structured mesh blocks in this MOAB instance, as ScdBox objects
00219     /** Return the structured blocks in this MOAB instance.  If these were not searched for
00220      * at instantiation time, then the search is done now.
00221      * \param boxes Vector of ScdBox objects representing structured mesh blocks
00222      */
00223     ErrorCode find_boxes( std::vector< ScdBox* >& boxes );
00224 
00225     //! Return all the structured mesh blocks in this MOAB instance, as entity set handles
00226     /** Return the structured blocks in this MOAB instance.  If these were not searched for
00227      * at instantiation time, then the search is done now.
00228      * \param boxes Range of entity set objects representing structured mesh blocks
00229      */
00230     ErrorCode find_boxes( Range& boxes );
00231 
00232     //! Return all the structured mesh blocks known by ScdInterface (does not search)
00233     /** Return the structured blocks in this ScdInterface instance.  Does not search for new boxes,
00234      * just returns the contents of the list.
00235      * \param boxes Structured boxes
00236      */
00237     ErrorCode get_boxes( std::vector< ScdBox* >& boxes );
00238 
00239     //! Return the tag marking the lower and upper corners of boxes
00240     /**
00241      * \param create_if_missing If the tag does not yet exist, create it
00242      */
00243     Tag box_dims_tag( bool create_if_missing = true );
00244 
00245     //! Return the tag marking the global lower and upper corners of boxes
00246     /**
00247      * \param create_if_missing If the tag does not yet exist, create it
00248      */
00249     Tag global_box_dims_tag( bool create_if_missing = true );
00250 
00251     //! Return the tag marking the partitioning method used to partition the box in parallel
00252     /**
00253      * \param create_if_missing If the tag does not yet exist, create it
00254      */
00255     Tag part_method_tag( bool create_if_missing = true );
00256 
00257     //! Return the tag marking whether box is periodic in i and j
00258     /**
00259      * \param create_if_missing If the tag does not yet exist, create it
00260      */
00261     Tag box_periodic_tag( bool create_if_missing = true );
00262 
00263     //! Return the tag marking the ScdBox for a set
00264     /**
00265      * \param create_if_missing If the tag does not yet exist, create it
00266      */
00267     Tag box_set_tag( bool create_if_missing = true );
00268 
00269     //! Return the ScdBox corresponding to the entity set passed in
00270     /** If the entity isn't a structured box set, NULL is returned.
00271      * \param eh Entity whose box is being queried
00272      */
00273     ScdBox* get_scd_box( EntityHandle eh );
00274 
00275     //! Compute a partition of structured parameter space
00276     /** Compute a partition of structured parameter space, based on data in the ScdParData
00277      * passed in.  Results are passed back in arguments, which application can set back into
00278      * par_data argument if they so desire.
00279      * \param np Number of processors
00280      * \param nr Rank of this processor
00281      * \param par_data ScdParData object that contains input global parameter space, desired
00282      *           partitioning method, and information about global periodicity.
00283      * \param ldims Local parameters for grid
00284      * \param lperiodic Whether or not a given dimension is locally periodic
00285      * \param pdims Number of procs in i, j, k directions
00286      */
00287     static ErrorCode compute_partition( int np,
00288                                         int nr,
00289                                         const ScdParData& par_data,
00290                                         int* ldims,
00291                                         int* lperiodic = NULL,
00292                                         int* pdims     = NULL );
00293 
00294     //! Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
00295     //! all 3 params
00296     /** Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
00297      * all 3 params \param np (in) Total # procs \param nr Processor from which neighbor is
00298      * requested \param spd (in) ScdParData containing part method, gdims, gperiodic data \param
00299      * dijk(*) (in) Direction being queried, = +/-1 or 0 \param pto (out) Processor holding the
00300      * destination part \param rdims(6) (out) Parametric min/max of destination part \param
00301      * facedims(6) (out) Parametric min/max of interface between pfrom and pto; if at the max in a
00302      * periodic direction, set to global min of that direction \param across_bdy(3) (out) If
00303      * across_bdy[i] is -1(1), interface with pto is across periodic lower(upper) bdy in parameter
00304      * i, 0 otherwise
00305      */
00306     static ErrorCode get_neighbor( int np,
00307                                    int nr,
00308                                    const ScdParData& spd,
00309                                    const int* const dijk,
00310                                    int& pto,
00311                                    int* rdims,
00312                                    int* facedims,
00313                                    int* across_bdy );
00314 
00315     //! Tag vertices with sharing data for parallel representations
00316     /** Given the ParallelComm object to use, tag the vertices shared with other processors
00317      */
00318     ErrorCode tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth );
00319 
00320     //! Tag vertices with sharing data for parallel representations
00321     /** Given the ParallelComm object to use, tag the vertices shared with other processors
00322      */
00323     ErrorCode tag_shared_vertices( ParallelComm* pcomm, ScdBox* box );
00324 
00325   protected:
00326     //! Remove the box from the list on ScdInterface
00327     ErrorCode remove_box( ScdBox* box );
00328 
00329     //! Add the box to the list on ScdInterface
00330     ErrorCode add_box( ScdBox* box );
00331 
00332   private:
00333     //! Create an entity set for a box, and tag with the parameters
00334     /** \param low Lower corner parameters for this box
00335      * \param high Upper corner parameters for this box
00336      * \param scd_set Entity set created
00337      * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
00338      * (s=[0,1,2])
00339      */
00340     ErrorCode create_box_set( const HomCoord& low,
00341                               const HomCoord& high,
00342                               EntityHandle& scd_set,
00343                               int* is_periodic = NULL );
00344 
00345     //! Compute a partition of structured parameter space
00346     /** Partitions the structured parametric space by partitioning j, k, or i only.
00347      * If j is greater than #procs, partition that, else k, else i.
00348      * For description of arguments, see ScdInterface::compute_partition.
00349      */
00350     inline static ErrorCode compute_partition_alljorkori( int np,
00351                                                           int nr,
00352                                                           const int gijk[6],
00353                                                           const int* const gperiodic,
00354                                                           int* lijk,
00355                                                           int* lperiodic,
00356                                                           int* pijk );
00357 
00358     //! Compute a partition of structured parameter space
00359     /** Partitions the structured parametric space by partitioning j, and possibly k,
00360      * seeking square regions of jk space
00361      * For description of arguments, see ScdInterface::compute_partition.
00362      */
00363     inline static ErrorCode compute_partition_alljkbal( int np,
00364                                                         int nr,
00365                                                         const int gijk[6],
00366                                                         const int* const gperiodic,
00367                                                         int* lijk,
00368                                                         int* lperiodic,
00369                                                         int* pijk );
00370 
00371     //! Compute a partition of structured parameter space
00372     /** Partitions the structured parametric space by seeking square ij partitions
00373      * For description of arguments, see ScdInterface::compute_partition.
00374      */
00375     inline static ErrorCode compute_partition_sqij( int np,
00376                                                     int nr,
00377                                                     const int gijk[6],
00378                                                     const int* const gperiodic,
00379                                                     int* lijk,
00380                                                     int* lperiodic,
00381                                                     int* pijk );
00382 
00383     //! Compute a partition of structured parameter space
00384     /** Partitions the structured parametric space by seeking square jk partitions
00385      * For description of arguments, see ScdInterface::compute_partition.
00386      */
00387     inline static ErrorCode compute_partition_sqjk( int np,
00388                                                     int nr,
00389                                                     const int gijk[6],
00390                                                     const int* const gperiodic,
00391                                                     int* lijk,
00392                                                     int* lperiodic,
00393                                                     int* pijk );
00394 
00395     //! Compute a partition of structured parameter space
00396     /** Partitions the structured parametric space by seeking square ijk partitions
00397      * For description of arguments, see ScdInterface::compute_partition.
00398      */
00399     inline static ErrorCode compute_partition_sqijk( int np,
00400                                                      int nr,
00401                                                      const int gijk[6],
00402                                                      const int* const gperiodic,
00403                                                      int* lijk,
00404                                                      int* lperiodic,
00405                                                      int* pijk );
00406 
00407     //! Get vertices shared with other processors
00408     /** Shared vertices returned as indices into each proc's handle space
00409      * \param box Box used to get parametric space info
00410      * \param procs Procs this proc shares vertices with
00411      * \param offsets Offsets into indices list for each proc
00412      * \param shared_indices local/remote indices of shared vertices
00413      */
00414     static ErrorCode get_shared_vertices( ParallelComm* pcomm,
00415                                           ScdBox* box,
00416                                           std::vector< int >& procs,
00417                                           std::vector< int >& offsets,
00418                                           std::vector< int >& shared_indices );
00419 
00420     static ErrorCode get_indices( const int* const ldims,
00421                                   const int* const rdims,
00422                                   const int* const across_bdy,
00423                                   int* face_dims,
00424                                   std::vector< int >& shared_indices );
00425 
00426     static ErrorCode get_neighbor_alljorkori( int np,
00427                                               int pfrom,
00428                                               const int* const gdims,
00429                                               const int* const gperiodic,
00430                                               const int* const dijk,
00431                                               int& pto,
00432                                               int* rdims,
00433                                               int* facedims,
00434                                               int* across_bdy );
00435 
00436     static ErrorCode get_neighbor_alljkbal( int np,
00437                                             int pfrom,
00438                                             const int* const gdims,
00439                                             const int* const gperiodic,
00440                                             const int* const dijk,
00441                                             int& pto,
00442                                             int* rdims,
00443                                             int* facedims,
00444                                             int* across_bdy );
00445 
00446     static ErrorCode get_neighbor_sqij( int np,
00447                                         int pfrom,
00448                                         const int* const gdims,
00449                                         const int* const gperiodic,
00450                                         const int* const dijk,
00451                                         int& pto,
00452                                         int* rdims,
00453                                         int* facedims,
00454                                         int* across_bdy );
00455 
00456     static ErrorCode get_neighbor_sqjk( int np,
00457                                         int pfrom,
00458                                         const int* const gdims,
00459                                         const int* const gperiodic,
00460                                         const int* const dijk,
00461                                         int& pto,
00462                                         int* rdims,
00463                                         int* facedims,
00464                                         int* across_bdy );
00465 
00466     static ErrorCode get_neighbor_sqijk( int np,
00467                                          int pfrom,
00468                                          const int* const gdims,
00469                                          const int* const gperiodic,
00470                                          const int* const dijk,
00471                                          int& pto,
00472                                          int* rdims,
00473                                          int* facedims,
00474                                          int* across_bdy );
00475 
00476     static int gtol( const int* gijk, int i, int j, int k );
00477 
00478     //! assign global ids to vertices in this box
00479     ErrorCode assign_global_ids( ScdBox* box );
00480 
00481     //! interface instance
00482     Interface* mbImpl;
00483 
00484     //! whether we've searched the database for boxes yet
00485     bool searchedBoxes;
00486 
00487     //! structured mesh blocks; stored as ScdBox objects, can get sets from those
00488     std::vector< ScdBox* > scdBoxes;
00489 
00490     //! tag representing whether box is periodic in i and j
00491     Tag boxPeriodicTag;
00492 
00493     //! tag representing box lower and upper corners
00494     Tag boxDimsTag;
00495 
00496     //! tag representing global lower and upper corners
00497     Tag globalBoxDimsTag;
00498 
00499     //! tag representing partition method
00500     Tag partMethodTag;
00501 
00502     //! tag pointing from set to ScdBox
00503     Tag boxSetTag;
00504 };
00505 
00506 class MOAB_EXPORT ScdBox
00507 {
00508     friend class ScdInterface;
00509 
00510   public:
00511     //! Destructor
00512     ~ScdBox();
00513 
00514     //! Return the ScdInterface responsible for this box
00515     inline ScdInterface* sc_impl() const;
00516 
00517     //! Add a vertex box to this box
00518     /* Add a vertex box to the element sequence referenced by this box.  The passed in vbox must
00519      * be a vertex box, with parametric extents no larger than that of this box.  This vbox is
00520      * oriented to this box by matching parameters from1-from3 in vbox to to1-to3 in this box.
00521      * If bb_input is true, only the part of the vertex sequence between bb_min and bb_max is
00522      * referenced \param vbox The vertex box being added to this box \param from1 1st reference
00523      * point on vbox \param to1 1st reference point on this box \param from2 2nd reference point on
00524      * vbox \param to2 2nd reference point on this box \param from3 3rd reference point on vbox
00525      * \param to3 3rd reference point on this box
00526      * \param bb_input If true, subsequent parameters list extents of vbox referenced
00527      * \param bb_min Lower corner of rectangle referenced
00528      * \param bb_max Upper corner of rectangle referenced
00529      */
00530     ErrorCode add_vbox( ScdBox* vbox,
00531                         HomCoord from1,
00532                         HomCoord to1,
00533                         HomCoord from2,
00534                         HomCoord to2,
00535                         HomCoord from3,
00536                         HomCoord to3,
00537                         bool bb_input          = false,
00538                         const HomCoord& bb_min = HomCoord::getUnitv( 0 ),
00539                         const HomCoord& bb_max = HomCoord::getUnitv( 0 ) );
00540     // const HomCoord &bb_min = HomCoord::unitv[0],
00541     // const HomCoord &bb_max = HomCoord::unitv[0]);
00542 
00543     //! Return whether this box has all its vertices defined
00544     /** Tests whether vertex boxs added with add_vbox have completely defined the vertex parametric
00545      * space for this box.
00546      *
00547      */
00548     bool boundary_complete() const;
00549 
00550     //! Return highest topological dimension of box
00551     inline int box_dimension() const;
00552 
00553     //! Starting vertex handle for this box
00554     inline EntityHandle start_vertex() const;
00555 
00556     //! Starting entity handle for this box
00557     /** If this is a vertex box, the start vertex handle is returned.
00558      */
00559     inline EntityHandle start_element() const;
00560 
00561     //! Return the number of elements in the box
00562     /* Number of elements is (boxSize[0]-1)(boxSize[1]-1)(boxSize[2]-1)
00563      */
00564     inline int num_elements() const;
00565 
00566     //! Return the number of vertices in the box
00567     /* Number of vertices is boxSize[0] * boxSize[1] * boxSize[2]
00568      */
00569     inline int num_vertices() const;
00570 
00571     //! Return the parametric coordinates for this box
00572     /**
00573      * \return IJK parameters of lower and upper corners
00574      */
00575     inline const int* box_dims() const;
00576 
00577     //! Return the lower corner parametric coordinates for this box
00578     inline HomCoord box_min() const;
00579 
00580     //! Return the upper corner parametric coordinates for this box
00581     inline HomCoord box_max() const;
00582 
00583     //! Return the parameter extents for this box
00584     inline HomCoord box_size() const;
00585 
00586     //! Return the parametric extents for this box
00587     /**
00588      * \param ijk IJK extents of this box
00589      */
00590     inline void box_size( int* ijk ) const;
00591 
00592     //! Return the parametric extents for this box
00593     /**
00594      * \param i I extent of this box
00595      * \param j J extent of this box
00596      * \param k K extent of this box
00597      */
00598     void box_size( int& i, int& j, int& k ) const;
00599 
00600     //! Get the element at the specified coordinates
00601     /**
00602      * \param ijk Parametric coordinates being evaluated
00603      */
00604     EntityHandle get_element( const HomCoord& ijk ) const;
00605 
00606     //! Get the element at the specified coordinates
00607     /**
00608      * \param i Parametric coordinates being evaluated
00609      * \param j Parametric coordinates being evaluated
00610      * \param k Parametric coordinates being evaluated
00611      */
00612     EntityHandle get_element( int i, int j = 0, int k = 0 ) const;
00613 
00614     //! Get the vertex at the specified coordinates
00615     /**
00616      * \param ijk Parametric coordinates being evaluated
00617      */
00618     EntityHandle get_vertex( const HomCoord& ijk ) const;
00619 
00620     //! Get the vertex at the specified coordinates
00621     /**
00622      * \param i Parametric coordinates being evaluated
00623      * \param j Parametric coordinates being evaluated
00624      * \param k Parametric coordinates being evaluated
00625      */
00626     EntityHandle get_vertex( int i, int j = 0, int k = 0 ) const;
00627 
00628     //! Get parametric coordinates of the specified entity
00629     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00630      * in this ScdBox.
00631      * \param ent Entity being queried
00632      * \param i Parametric coordinates returned
00633      * \param j Parametric coordinates returned
00634      * \param k Parametric coordinates returned
00635      * \param dir Parametric coordinate direction returned (in case of getting adjacent
00636      *            edges (2d, 3d) or faces (3d); not modified otherwise
00637      */
00638     ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const;
00639 
00640     //! Get parametric coordinates of the specified entity, intermediate entities not allowed (no
00641     //! dir parameter)
00642     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00643      * in this ScdBox, or MB_FAILURE if the entity is an intermediate-dimension entity.
00644      * \param ent Entity being queried
00645      * \param i Parametric coordinates returned
00646      * \param j Parametric coordinates returned
00647      * \param k Parametric coordinates returned
00648      */
00649     ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k ) const;
00650 
00651     //! Get parametric coordinates of the specified entity
00652     /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
00653      * in this ScdBox.
00654      * \param ent Entity being queried
00655      * \param ijkd Parametric coordinates returned (including direction, in case of
00656      *            getting adjacent edges (2d, 3d) or faces (3d))
00657      */
00658     ErrorCode get_params( EntityHandle ent, HomCoord& ijkd ) const;
00659 
00660     /** \brief Get the adjacent edge or face at a parametric location
00661      * This function gets the left (i=0), front (j=0), or bottom (k=0) edge or face for a parametric
00662      * element. Left, front, or bottom is indicated by dir = 0, 1, or 2, resp.  All edges and faces
00663      * in a structured mesh block can be accessed using these parameters. \param dim Dimension of
00664      * adjacent entity being requested \param i Parametric coordinates of cell being evaluated
00665      * \param j Parametric coordinates of cell being evaluated
00666      * \param k Parametric coordinates of cell being evaluated
00667      * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d)
00668      * \param ent Entity returned from this function
00669      * \param create_if_missing If true, creates the entity if it doesn't already exist
00670      */
00671     ErrorCode get_adj_edge_or_face( int dim,
00672                                     int i,
00673                                     int j,
00674                                     int k,
00675                                     int dir,
00676                                     EntityHandle& ent,
00677                                     bool create_if_missing = true ) const;
00678 
00679     //! Return whether the box contains the parameters passed in
00680     /**
00681      * \param i Parametric coordinates being evaluated
00682      * \param j Parametric coordinates being evaluated
00683      * \param k Parametric coordinates being evaluated
00684      */
00685     bool contains( int i, int j, int k ) const;
00686 
00687     //! Return whether the box contains the parameters passed in
00688     /**
00689      * \param i Parametric coordinates being evaluated
00690      * \param j Parametric coordinates being evaluated
00691      * \param k Parametric coordinates being evaluated
00692      */
00693     bool contains( const HomCoord& ijk ) const;
00694 
00695     //! Set/Get the entity set representing the box
00696     void box_set( EntityHandle this_set );
00697     EntityHandle box_set();
00698 
00699     //! Get coordinate arrays for vertex coordinates for a structured block
00700     /** Returns error if there isn't a single vertex sequence associated with this structured block
00701      * \param xc X coordinate array pointer returned
00702      * \param yc Y coordinate array pointer returned
00703      * \param zc Z coordinate array pointer returned
00704      */
00705     ErrorCode get_coordinate_arrays( double*& xc, double*& yc, double*& zc );
00706 
00707     //! Get read-only coordinate arrays for vertex coordinates for a structured block
00708     /** Returns error if there isn't a single vertex sequence associated with this structured block
00709      * \param xc X coordinate array pointer returned
00710      * \param yc Y coordinate array pointer returned
00711      * \param zc Z coordinate array pointer returned
00712      */
00713     ErrorCode get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const;
00714 
00715     //! Return whether box is locally periodic in i
00716     /** Return whether box is locally periodic in i
00717      * \return True if box is locally periodic in i direction
00718      */
00719     bool locally_periodic_i() const;
00720 
00721     //! Return whether box is locally periodic in j
00722     /** Return whether box is locally periodic in j
00723      * \return True if box is locally periodic in j direction
00724      */
00725     bool locally_periodic_j() const;
00726 
00727     //! Return whether box is locally periodic in k
00728     /** Return whether box is locally periodic in k
00729      * \return True if box is locally periodic in k direction
00730      */
00731     bool locally_periodic_k() const;
00732 
00733     //! Set local periodicity
00734     /**
00735      * \param lperiodic Vector of ijk periodicities to set this box to
00736      */
00737     void locally_periodic( bool lperiodic[3] );
00738 
00739     //! Get local periodicity
00740     /**
00741      * \return Vector of ijk periodicities for this box
00742      */
00743     const int* locally_periodic() const;
00744 
00745     //! Return parallel data
00746     /** Return parallel data, if there is any
00747      * \return par_data Parallel data set on this box
00748      */
00749     ScdParData& par_data()
00750     {
00751         return parData;
00752     }
00753 
00754     //! Return parallel data
00755     /** Return parallel data, if there is any
00756      * \return par_data Parallel data set on this box
00757      */
00758     const ScdParData& par_data() const
00759     {
00760         return parData;
00761     }
00762 
00763     //! set parallel data
00764     /** Set parallel data for this box
00765      * \param par_data Parallel data to be set on this box
00766      */
00767     void par_data( const ScdParData& par_datap )
00768     {
00769         parData = par_datap;
00770     }
00771 
00772   private:
00773     //! Constructor
00774     /** Create a structured box instance; this constructor is private because it should only be
00775      * called from ScdInterface, a friend class.  This constructor takes two sequences, one of which
00776      * can be NULL.  If both sequences come in non-NULL, the first should be a VertexSequence*
00777      * corresponding to a structured vertex sequence and the second should be a
00778      * StructuredElementSeq*.  If the 2nd is NULL, the first can be either of those types.  The
00779      * other members of this class are taken from the sequences (e.g. parametric space) or the box
00780      * set argument.  Tags on the box set should be set from the caller. \param sc_impl A
00781      * ScdInterface instance \param box_set Entity set representing this rectangle \param seq1 An
00782      * EntitySequence (see ScdBox description) \param seq2 An EntitySequence (see ScdBox
00783      * description), or NULL
00784      */
00785     ScdBox( ScdInterface* sc_impl, EntityHandle box_set, EntitySequence* seq1, EntitySequence* seq2 = NULL );
00786 
00787     //! function to get vertex handle directly from sequence
00788     /** \param i Parameter being queried
00789      * \param j Parameter being queried
00790      * \param k Parameter being queried
00791      */
00792     EntityHandle get_vertex_from_seq( int i, int j, int k ) const;
00793 
00794     //! set the vertex sequence
00795     ErrorCode vert_dat( ScdVertexData* vert_dat );
00796 
00797     //! get the vertex sequence
00798     ScdVertexData* vert_dat() const;
00799 
00800     //! set the element sequence
00801     ErrorCode elem_seq( EntitySequence* elem_seq );
00802 
00803     //! get the element sequence
00804     StructuredElementSeq* elem_seq() const;
00805 
00806     //! Set the starting vertex handle for this box
00807     void start_vertex( EntityHandle startv );
00808 
00809     //! Set the starting entity handle for this box
00810     void start_element( EntityHandle starte );
00811 
00812     //! interface instance
00813     ScdInterface* scImpl;
00814 
00815     //! entity set representing this box
00816     EntityHandle boxSet;
00817 
00818     //! vertex sequence this box represents, if there's only one, otherwise they're
00819     //! retrieved from the element sequence
00820     ScdVertexData* vertDat;
00821 
00822     //! element sequence this box represents
00823     StructuredElementSeq* elemSeq;
00824 
00825     //! starting vertex handle for this box
00826     EntityHandle startVertex;
00827 
00828     //! starting element handle for this box
00829     EntityHandle startElem;
00830 
00831     //! lower and upper corners
00832     int boxDims[6];
00833 
00834     //! is locally periodic in i or j or k
00835     int locallyPeriodic[3];
00836 
00837     //! parallel data associated with this box, if any
00838     ScdParData parData;
00839 
00840     //! parameter extents
00841     HomCoord boxSize;
00842 
00843     //! convenience parameters, (boxSize[1]-1)*(boxSize[0]-1) and boxSize[0]-1
00844     int boxSizeIJ;
00845     int boxSizeIJM1;
00846     int boxSizeIM1;
00847 };
00848 
00849 inline ErrorCode ScdInterface::compute_partition( int np,
00850                                                   int nr,
00851                                                   const ScdParData& par_data,
00852                                                   int* ldims,
00853                                                   int* lperiodic,
00854                                                   int* pdims )
00855 {
00856     ErrorCode rval = MB_SUCCESS;
00857     switch( par_data.partMethod )
00858     {
00859         case ScdParData::ALLJORKORI:
00860         case -1:
00861             rval = compute_partition_alljorkori( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00862             break;
00863         case ScdParData::ALLJKBAL:
00864             rval = compute_partition_alljkbal( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00865             break;
00866         case ScdParData::SQIJ:
00867             rval = compute_partition_sqij( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00868             break;
00869         case ScdParData::SQJK:
00870             rval = compute_partition_sqjk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00871             break;
00872         case ScdParData::SQIJK:
00873             rval = compute_partition_sqijk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
00874             break;
00875         default:
00876             rval = MB_FAILURE;
00877             break;
00878     }
00879 
00880     return rval;
00881 }
00882 
00883 inline ErrorCode ScdInterface::compute_partition_alljorkori( int np,
00884                                                              int nr,
00885                                                              const int gijk[6],
00886                                                              const int* const gperiodic,
00887                                                              int* ldims,
00888                                                              int* lperiodic,
00889                                                              int* pijk )
00890 {
00891     // partition *the elements* over the parametric space; 1d partition for now, in the j, k, or i
00892     // parameters
00893     int tmp_lp[3], tmp_pijk[3];
00894     if( !lperiodic ) lperiodic = tmp_lp;
00895     if( !pijk ) pijk = tmp_pijk;
00896 
00897     for( int i = 0; i < 3; i++ )
00898         lperiodic[i] = gperiodic[i];
00899 
00900     if( np == 1 )
00901     {
00902         if( ldims )
00903         {
00904             ldims[0] = gijk[0];
00905             ldims[3] = gijk[3];
00906             ldims[1] = gijk[1];
00907             ldims[4] = gijk[4];
00908             ldims[2] = gijk[2];
00909             ldims[5] = gijk[5];
00910         }
00911         pijk[0] = pijk[1] = pijk[2] = 1;
00912     }
00913     else
00914     {
00915         if( gijk[4] - gijk[1] > np )
00916         {
00917             // partition j over procs
00918             int dj    = ( gijk[4] - gijk[1] ) / np;
00919             int extra = ( gijk[4] - gijk[1] ) % np;
00920             ldims[1]  = gijk[1] + nr * dj + std::min( nr, extra );
00921             ldims[4]  = ldims[1] + dj + ( nr < extra ? 1 : 0 );
00922 
00923             if( gperiodic[1] && np > 1 )
00924             {
00925                 lperiodic[1] = 0;
00926                 ldims[4]++;
00927             }
00928 
00929             ldims[2] = gijk[2];
00930             ldims[5] = gijk[5];
00931             ldims[0] = gijk[0];
00932             ldims[3] = gijk[3];
00933             pijk[0] = pijk[2] = 1;
00934             pijk[1]           = np;
00935         }
00936         else if( gijk[5] - gijk[2] > np )
00937         {
00938             // partition k over procs
00939             int dk    = ( gijk[5] - gijk[2] ) / np;
00940             int extra = ( gijk[5] - gijk[2] ) % np;
00941             ldims[2]  = gijk[2] + nr * dk + std::min( nr, extra );
00942             ldims[5]  = ldims[2] + dk + ( nr < extra ? 1 : 0 );
00943 
00944             ldims[1] = gijk[1];
00945             ldims[4] = gijk[4];
00946             ldims[0] = gijk[0];
00947             ldims[3] = gijk[3];
00948             pijk[0] = pijk[1] = 1;
00949             pijk[2]           = np;
00950         }
00951         else if( gijk[3] - gijk[0] > np )
00952         {
00953             // partition i over procs
00954             int di    = ( gijk[3] - gijk[0] ) / np;
00955             int extra = ( gijk[3] - gijk[0] ) % np;
00956             ldims[0]  = gijk[0] + nr * di + std::min( nr, extra );
00957             ldims[3]  = ldims[0] + di + ( nr < extra ? 1 : 0 );
00958 
00959             if( gperiodic[0] && np > 1 )
00960             {
00961                 lperiodic[0] = 0;
00962                 ldims[3]++;
00963             }
00964 
00965             ldims[2] = gijk[2];
00966             ldims[5] = gijk[5];
00967             ldims[1] = gijk[1];
00968             ldims[4] = gijk[4];
00969 
00970             pijk[1] = pijk[2] = 1;
00971             pijk[0]           = np;
00972         }
00973         else
00974         {
00975             // Couldn't find a suitable partition...
00976             return MB_FAILURE;
00977         }
00978     }
00979 
00980     return MB_SUCCESS;
00981 }
00982 
00983 inline ErrorCode ScdInterface::compute_partition_alljkbal( int np,
00984                                                            int nr,
00985                                                            const int gijk[6],
00986                                                            const int* const gperiodic,
00987                                                            int* ldims,
00988                                                            int* lperiodic,
00989                                                            int* pijk )
00990 {
00991     int tmp_lp[3], tmp_pijk[3];
00992     if( !lperiodic ) lperiodic = tmp_lp;
00993     if( !pijk ) pijk = tmp_pijk;
00994 
00995     for( int i = 0; i < 3; i++ )
00996         lperiodic[i] = gperiodic[i];
00997 
00998     if( np == 1 )
00999     {
01000         if( ldims )
01001         {
01002             ldims[0] = gijk[0];
01003             ldims[3] = gijk[3];
01004             ldims[1] = gijk[1];
01005             ldims[4] = gijk[4];
01006             ldims[2] = gijk[2];
01007             ldims[5] = gijk[5];
01008         }
01009         pijk[0] = pijk[1] = pijk[2] = 1;
01010     }
01011     else
01012     {
01013         // improved, possibly 2-d partition
01014         std::vector< double > kfactors;
01015         kfactors.push_back( 1 );
01016         int K = gijk[5] - gijk[2];
01017         for( int i = 2; i < K; i++ )
01018             if( !( K % i ) && !( np % i ) ) kfactors.push_back( i );
01019         kfactors.push_back( K );
01020 
01021         // compute the ideal nj and nk
01022         int J          = gijk[4] - gijk[1];
01023         double njideal = sqrt( ( (double)( np * J ) ) / ( (double)K ) );
01024         double nkideal = ( njideal * K ) / J;
01025 
01026         int nk, nj;
01027         if( nkideal < 1.0 )
01028         {
01029             nk = 1;
01030             nj = np;
01031         }
01032         else
01033         {
01034             std::vector< double >::iterator vit = std::lower_bound( kfactors.begin(), kfactors.end(), nkideal );
01035             if( vit == kfactors.begin() )
01036                 nk = 1;
01037             else
01038                 nk = (int)*( --vit );
01039             nj = np / nk;
01040         }
01041 
01042         int dk = K / nk;
01043         int dj = J / nj;
01044 
01045         ldims[2] = gijk[2] + ( nr % nk ) * dk;
01046         ldims[5] = ldims[2] + dk;
01047 
01048         int extra = J % nj;
01049 
01050         ldims[1] = gijk[1] + ( nr / nk ) * dj + std::min( nr / nk, extra );
01051         ldims[4] = ldims[1] + dj + ( nr / nk < extra ? 1 : 0 );
01052 
01053         ldims[0] = gijk[0];
01054         ldims[3] = gijk[3];
01055 
01056         if( gperiodic[1] && np > 1 )
01057         {
01058             lperiodic[1] = 0;
01059             if( nr / nk == nj - 1 )
01060             {
01061                 ldims[1]++;
01062             }
01063         }
01064 
01065         pijk[0] = 1;
01066         pijk[1] = nj;
01067         pijk[2] = nk;
01068     }
01069 
01070     return MB_SUCCESS;
01071 }
01072 
01073 inline ErrorCode ScdInterface::compute_partition_sqij( int np,
01074                                                        int nr,
01075                                                        const int gijk[6],
01076                                                        const int* const gperiodic,
01077                                                        int* ldims,
01078                                                        int* lperiodic,
01079                                                        int* pijk )
01080 {
01081     int tmp_lp[3], tmp_pijk[3];
01082     if( !lperiodic ) lperiodic = tmp_lp;
01083     if( !pijk ) pijk = tmp_pijk;
01084 
01085     // square IxJ partition
01086 
01087     for( int i = 0; i < 3; i++ )
01088         lperiodic[i] = gperiodic[i];
01089 
01090     if( np == 1 )
01091     {
01092         if( ldims )
01093         {
01094             ldims[0] = gijk[0];
01095             ldims[3] = gijk[3];
01096             ldims[1] = gijk[1];
01097             ldims[4] = gijk[4];
01098             ldims[2] = gijk[2];
01099             ldims[5] = gijk[5];
01100         }
01101         pijk[0] = pijk[1] = pijk[2] = 1;
01102     }
01103     else
01104     {
01105         std::vector< double > pfactors, ppfactors;
01106         for( int i = 2; i <= np / 2; i++ )
01107             if( !( np % i ) )
01108             {
01109                 pfactors.push_back( i );
01110                 ppfactors.push_back( ( (double)( i * i ) ) / np );
01111             }
01112         pfactors.push_back( np );
01113         ppfactors.push_back( (double)np );
01114 
01115         // ideally, Px/Py = I/J
01116         double ijratio = ( (double)( gijk[3] - gijk[0] ) ) / ( (double)( gijk[4] - gijk[1] ) );
01117 
01118         unsigned int ind                        = 0;
01119         std::vector< double >::iterator optimal = std::lower_bound( ppfactors.begin(), ppfactors.end(), ijratio );
01120         if( optimal == ppfactors.end() )
01121         {
01122             ind = ppfactors.size() - 1;
01123         }
01124         else
01125         {
01126             ind = optimal - ppfactors.begin();
01127             if( ind && fabs( ppfactors[ind - 1] - ijratio ) < fabs( ppfactors[ind] - ijratio ) ) ind--;
01128         }
01129 
01130         // VARIABLES DESCRIBING THE MESH:
01131         // pi, pj = # procs in i and j directions
01132         // nri, nrj = my proc's position in i, j directions
01133         // I, J = # edges/elements in i, j directions
01134         // iextra, jextra = # procs having extra edge in i/j direction
01135         // top_i, top_j = if true, I'm the last proc in the i/j direction
01136         // i, j = # edges locally in i/j direction, *not* including one for iextra/jextra
01137         int pi = pfactors[ind];
01138         int pj = np / pi;
01139 
01140         int I = ( gijk[3] - gijk[0] ), J = ( gijk[4] - gijk[1] );
01141         int iextra = I % pi, jextra = J % pj, i = I / pi, j = J / pj;
01142         int nri = nr % pi, nrj = nr / pi;
01143 
01144         if( ldims )
01145         {
01146             ldims[0] = gijk[0] + i * nri + std::min( iextra, nri );
01147             ldims[3] = ldims[0] + i + ( nri < iextra ? 1 : 0 );
01148             ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
01149             ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
01150 
01151             ldims[2] = gijk[2];
01152             ldims[5] = gijk[5];
01153 
01154             if( gperiodic[0] && pi > 1 )
01155             {
01156                 lperiodic[0] = 0;
01157                 if( nri == pi - 1 ) ldims[3]++;
01158             }
01159             if( gperiodic[1] && pj > 1 )
01160             {
01161                 lperiodic[1] = 0;
01162                 if( nrj == pj - 1 ) ldims[4]++;
01163             }
01164         }
01165 
01166         pijk[0] = pi;
01167         pijk[1] = pj;
01168         pijk[2] = 1;
01169     }
01170 
01171     return MB_SUCCESS;
01172 }
01173 
01174 inline ErrorCode ScdInterface::compute_partition_sqjk( int np,
01175                                                        int nr,
01176                                                        const int gijk[6],
01177                                                        const int* const gperiodic,
01178                                                        int* ldims,
01179                                                        int* lperiodic,
01180                                                        int* pijk )
01181 {
01182     int tmp_lp[3], tmp_pijk[3];
01183     if( !lperiodic ) lperiodic = tmp_lp;
01184     if( !pijk ) pijk = tmp_pijk;
01185 
01186     // square JxK partition
01187     for( int i = 0; i < 3; i++ )
01188         lperiodic[i] = gperiodic[i];
01189 
01190     if( np == 1 )
01191     {
01192         if( ldims )
01193         {
01194             ldims[0] = gijk[0];
01195             ldims[3] = gijk[3];
01196             ldims[1] = gijk[1];
01197             ldims[4] = gijk[4];
01198             ldims[2] = gijk[2];
01199             ldims[5] = gijk[5];
01200         }
01201         pijk[0] = pijk[1] = pijk[2] = 1;
01202     }
01203     else
01204     {
01205         std::vector< double > pfactors, ppfactors;
01206         for( int p = 2; p <= np; p++ )
01207             if( !( np % p ) )
01208             {
01209                 pfactors.push_back( p );
01210                 ppfactors.push_back( ( (double)( p * p ) ) / np );
01211             }
01212 
01213         // ideally, Pj/Pk = J/K
01214         int pj, pk;
01215         if( gijk[5] == gijk[2] )
01216         {
01217             pk = 1;
01218             pj = np;
01219         }
01220         else
01221         {
01222             double jkratio = ( (double)( gijk[4] - gijk[1] ) ) / ( (double)( gijk[5] - gijk[2] ) );
01223 
01224             std::vector< double >::iterator vit = std::lower_bound( ppfactors.begin(), ppfactors.end(), jkratio );
01225             if( vit == ppfactors.end() )
01226                 --vit;
01227             else if( vit != ppfactors.begin() && fabs( *( vit - 1 ) - jkratio ) < fabs( ( *vit ) - jkratio ) )
01228                 --vit;
01229             int ind = vit - ppfactors.begin();
01230 
01231             pj = 1;
01232             if( ind >= 0 && !pfactors.empty() ) pfactors[ind];
01233             pk = np / pj;
01234         }
01235 
01236         int K = ( gijk[5] - gijk[2] ), J = ( gijk[4] - gijk[1] );
01237         int jextra = J % pj, kextra = K % pk, j = J / pj, k = K / pk;
01238         int nrj = nr % pj, nrk = nr / pj;
01239         ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
01240         ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
01241         ldims[2] = gijk[2] + k * nrk + std::min( kextra, nrk );
01242         ldims[5] = ldims[2] + k + ( nrk < kextra ? 1 : 0 );
01243 
01244         ldims[0] = gijk[0];
01245         ldims[3] = gijk[3];
01246 
01247         if( gperiodic[1] && pj > 1 )
01248         {
01249             lperiodic[1] = 0;
01250             if( nrj == pj - 1 ) ldims[4]++;
01251         }
01252 
01253         pijk[0] = 1;
01254         pijk[1] = pj;
01255         pijk[2] = pk;
01256     }
01257 
01258     return MB_SUCCESS;
01259 }
01260 
01261 inline ErrorCode ScdInterface::compute_partition_sqijk( int np,
01262                                                         int nr,
01263                                                         const int* const gijk,
01264                                                         const int* const gperiodic,
01265                                                         int* ldims,
01266                                                         int* lperiodic,
01267                                                         int* pijk )
01268 {
01269     if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE;
01270 
01271     int tmp_lp[3], tmp_pijk[3];
01272     if( !lperiodic ) lperiodic = tmp_lp;
01273     if( !pijk ) pijk = tmp_pijk;
01274 
01275     // square IxJxK partition
01276 
01277     for( int i = 0; i < 3; i++ )
01278         lperiodic[i] = gperiodic[i];
01279 
01280     if( np == 1 )
01281     {
01282         if( ldims )
01283             for( int i = 0; i < 6; i++ )
01284                 ldims[i] = gijk[i];
01285         pijk[0] = pijk[1] = pijk[2] = 1;
01286         return MB_SUCCESS;
01287     }
01288 
01289     std::vector< int > pfactors;
01290     pfactors.push_back( 1 );
01291     for( int i = 2; i <= np / 2; i++ )
01292         if( !( np % i ) ) pfactors.push_back( i );
01293     pfactors.push_back( np );
01294 
01295     // test for IJ, JK, IK
01296     int IJK[3], dIJK[3];
01297     for( int i = 0; i < 3; i++ )
01298         IJK[i] = std::max( gijk[3 + i] - gijk[i], 1 );
01299     // order IJK from lo to hi
01300     int lo = 0, hi = 0;
01301     for( int i = 1; i < 3; i++ )
01302     {
01303         if( IJK[i] < IJK[lo] ) lo = i;
01304         if( IJK[i] > IJK[hi] ) hi = i;
01305     }
01306     if( lo == hi ) hi = ( lo + 1 ) % 3;
01307     int mid = 3 - lo - hi;
01308     // search for perfect subdivision of np that balances #cells
01309     int perfa_best = -1, perfb_best = -1;
01310     double ratio = 0.0;
01311     for( int po = 0; po < (int)pfactors.size(); po++ )
01312     {
01313         for( int pi = po; pi < (int)pfactors.size() && np / ( pfactors[po] * pfactors[pi] ) >= pfactors[pi]; pi++ )
01314         {
01315             int p3 =
01316                 std::find( pfactors.begin(), pfactors.end(), np / ( pfactors[po] * pfactors[pi] ) ) - pfactors.begin();
01317             if( p3 == (int)pfactors.size() || pfactors[po] * pfactors[pi] * pfactors[p3] != np )
01318                 continue;  // po*pi should exactly factor np
01319             assert( po <= pi && pi <= p3 );
01320             // by definition, po <= pi <= p3
01321             double minl =
01322                        std::min( std::min( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] ),
01323                    maxl =
01324                        std::max( std::max( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] );
01325             if( minl / maxl > ratio )
01326             {
01327                 ratio      = minl / maxl;
01328                 perfa_best = po;
01329                 perfb_best = pi;
01330             }
01331         }
01332     }
01333     if( perfa_best == -1 || perfb_best == -1 ) return MB_FAILURE;
01334 
01335     // VARIABLES DESCRIBING THE MESH:
01336     // pijk[i] = # procs in direction i
01337     // numr[i] = my proc's position in direction i
01338     // dIJK[i] = # edges/elements in direction i
01339     // extra[i]= # procs having extra edge in direction i
01340     // top[i] = if true, I'm the last proc in direction i
01341 
01342     pijk[lo]     = pfactors[perfa_best];
01343     pijk[mid]    = pfactors[perfb_best];
01344     pijk[hi]     = ( np / ( pfactors[perfa_best] * pfactors[perfb_best] ) );
01345     int extra[3] = { 0, 0, 0 }, numr[3];
01346     for( int i = 0; i < 3; i++ )
01347     {
01348         dIJK[i]  = IJK[i] / pijk[i];
01349         extra[i] = IJK[i] % pijk[i];
01350     }
01351     numr[2] = nr / ( pijk[0] * pijk[1] );
01352     int rem = nr % ( pijk[0] * pijk[1] );
01353     numr[1] = rem / pijk[0];
01354     numr[0] = rem % pijk[0];
01355     for( int i = 0; i < 3; i++ )
01356     {
01357         extra[i]     = IJK[i] % dIJK[i];
01358         ldims[i]     = gijk[i] + numr[i] * dIJK[i] + std::min( extra[i], numr[i] );
01359         ldims[3 + i] = ldims[i] + dIJK[i] + ( numr[i] < extra[i] ? 1 : 0 );
01360     }
01361 
01362     return MB_SUCCESS;
01363 }
01364 
01365 inline int ScdInterface::gtol( const int* gijk, int i, int j, int k )
01366 {
01367     return ( ( k - gijk[2] ) * ( gijk[3] - gijk[0] + 1 ) * ( gijk[4] - gijk[1] + 1 ) +
01368              ( j - gijk[1] ) * ( gijk[3] - gijk[0] + 1 ) + i - gijk[0] );
01369 }
01370 
01371 inline ErrorCode ScdInterface::get_indices( const int* const ldims,
01372                                             const int* const rdims,
01373                                             const int* const across_bdy,
01374                                             int* face_dims,
01375                                             std::vector< int >& shared_indices )
01376 {
01377     // check for going across periodic bdy and face_dims not in my ldims (I'll always be on top in
01378     // that case)...
01379     if( across_bdy[0] > 0 && face_dims[0] != ldims[3] )
01380         face_dims[0] = face_dims[3] = ldims[3];
01381     else if( across_bdy[0] < 0 && face_dims[0] != ldims[0] )
01382         face_dims[0] = face_dims[3] = ldims[0];
01383     if( across_bdy[1] > 0 && face_dims[1] != ldims[4] )
01384         face_dims[1] = face_dims[4] = ldims[4];
01385     else if( across_bdy[1] < 0 && face_dims[1] != ldims[1] )
01386         face_dims[0] = face_dims[3] = ldims[1];
01387 
01388     for( int k = face_dims[2]; k <= face_dims[5]; k++ )
01389         for( int j = face_dims[1]; j <= face_dims[4]; j++ )
01390             for( int i = face_dims[0]; i <= face_dims[3]; i++ )
01391                 shared_indices.push_back( gtol( ldims, i, j, k ) );
01392 
01393     if( across_bdy[0] > 0 && face_dims[0] != rdims[0] )
01394         face_dims[0] = face_dims[3] = rdims[0];
01395     else if( across_bdy[0] < 0 && face_dims[0] != rdims[3] )
01396         face_dims[0] = face_dims[3] = rdims[3];
01397     if( across_bdy[1] > 0 && face_dims[1] != rdims[1] )
01398         face_dims[1] = face_dims[4] = rdims[1];
01399     else if( across_bdy[1] < 0 && face_dims[1] != rdims[4] )
01400         face_dims[0] = face_dims[3] = rdims[4];
01401 
01402     for( int k = face_dims[2]; k <= face_dims[5]; k++ )
01403         for( int j = face_dims[1]; j <= face_dims[4]; j++ )
01404             for( int i = face_dims[0]; i <= face_dims[3]; i++ )
01405                 shared_indices.push_back( gtol( rdims, i, j, k ) );
01406 
01407     return MB_SUCCESS;
01408 }
01409 
01410 inline ErrorCode ScdInterface::get_neighbor( int np,
01411                                              int pfrom,
01412                                              const ScdParData& spd,
01413                                              const int* const dijk,
01414                                              int& pto,
01415                                              int* rdims,
01416                                              int* facedims,
01417                                              int* across_bdy )
01418 {
01419     if( !dijk[0] && !dijk[1] && !dijk[2] )
01420     {
01421         // not going anywhere, return
01422         pto = -1;
01423         return MB_SUCCESS;
01424     }
01425 
01426     switch( spd.partMethod )
01427     {
01428         case ScdParData::ALLJORKORI:
01429         case -1:
01430             return get_neighbor_alljorkori( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims,
01431                                             across_bdy );
01432         case ScdParData::ALLJKBAL:
01433             return get_neighbor_alljkbal( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01434         case ScdParData::SQIJ:
01435             return get_neighbor_sqij( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01436         case ScdParData::SQJK:
01437             return get_neighbor_sqjk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01438         case ScdParData::SQIJK:
01439             return get_neighbor_sqijk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
01440         default:
01441             break;
01442     }
01443 
01444     return MB_FAILURE;
01445 }
01446 
01447 inline ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth )
01448 {
01449     ScdBox* box = get_scd_box( seth );
01450     if( !box )
01451     {
01452         // look for contained boxes
01453         Range tmp_range;
01454         ErrorCode rval = mbImpl->get_entities_by_type( seth, MBENTITYSET, tmp_range );
01455         if( MB_SUCCESS != rval ) return rval;
01456         for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
01457         {
01458             box = get_scd_box( *rit );
01459             if( box ) break;
01460         }
01461     }
01462 
01463     if( !box ) return MB_FAILURE;
01464 
01465     return tag_shared_vertices( pcomm, box );
01466 }
01467 
01468 inline ScdInterface* ScdBox::sc_impl() const
01469 {
01470     return scImpl;
01471 }
01472 
01473 inline EntityHandle ScdBox::start_vertex() const
01474 {
01475     return startVertex;
01476 }
01477 
01478 inline void ScdBox::start_vertex( EntityHandle startv )
01479 {
01480     startVertex = startv;
01481 }
01482 
01483 inline EntityHandle ScdBox::start_element() const
01484 {
01485     return startElem;
01486 }
01487 
01488 inline void ScdBox::start_element( EntityHandle starte )
01489 {
01490     startElem = starte;
01491 }
01492 
01493 inline int ScdBox::num_elements() const
01494 {
01495     if( !startElem ) return 0;  // not initialized yet
01496 
01497     /* for a structured mesh, total number of elements is obtained by multiplying
01498         number of elements in each direction
01499       number of elements in each direction is given by number of vertices in that direction minus 1
01500       if periodic in that direction, the last vertex is the same as first one, count one more
01501       element
01502       */
01503     int num_e_i = ( -1 == boxSize[0] || 1 == boxSize[0] ) ? 1 : boxSize[0] - 1;
01504     if( locallyPeriodic[0] ) ++num_e_i;
01505 
01506     int num_e_j = ( -1 == boxSize[1] || 1 == boxSize[1] ) ? 1 : boxSize[1] - 1;
01507     if( locallyPeriodic[1] ) ++num_e_j;
01508 
01509     int num_e_k = ( -1 == boxSize[2] || 1 == boxSize[2] ) ? 1 : boxSize[2] - 1;
01510     if( locallyPeriodic[2] ) ++num_e_k;
01511 
01512     return num_e_i * num_e_j * num_e_k;
01513 }
01514 
01515 inline int ScdBox::num_vertices() const
01516 {
01517     return boxSize[0] * ( !boxSize[1] ? 1 : boxSize[1] ) * ( !boxSize[2] ? 1 : boxSize[2] );
01518 }
01519 
01520 inline const int* ScdBox::box_dims() const
01521 {
01522     return boxDims;
01523 }
01524 
01525 inline HomCoord ScdBox::box_min() const
01526 {
01527     return HomCoord( boxDims, 3 );
01528 }
01529 
01530 inline HomCoord ScdBox::box_max() const
01531 {
01532     return HomCoord( boxDims + 3, 3 );
01533 }
01534 
01535 inline HomCoord ScdBox::box_size() const
01536 {
01537     return boxSize;
01538 }
01539 
01540 inline void ScdBox::box_size( int* ijk ) const
01541 {
01542     ijk[0] = boxSize[0];
01543     ijk[1] = boxSize[1];
01544     ijk[2] = boxSize[2];
01545 }
01546 
01547 inline void ScdBox::box_size( int& i, int& j, int& k ) const
01548 {
01549     i = boxSize[0];
01550     j = boxSize[1];
01551     k = boxSize[2];
01552 }
01553 
01554 inline EntityHandle ScdBox::get_element( int i, int j, int k ) const
01555 {
01556     return ( !startElem
01557                  ? 0
01558                  : startElem + ( k - boxDims[2] ) * boxSizeIJM1 + ( j - boxDims[1] ) * boxSizeIM1 + i - boxDims[0] );
01559 }
01560 
01561 inline EntityHandle ScdBox::get_element( const HomCoord& ijk ) const
01562 {
01563     return get_element( ijk[0], ijk[1], ijk[2] );
01564 }
01565 
01566 inline EntityHandle ScdBox::get_vertex( int i, int j, int k ) const
01567 {
01568     return ( vertDat
01569                  ? startVertex + ( boxDims[2] == -1 && boxDims[5] == -1 ? 0 : ( k - boxDims[2] ) ) * boxSizeIJ +
01570                        ( boxDims[1] == -1 && boxDims[4] == -1 ? 0 : ( j - boxDims[1] ) ) * boxSize[0] + i - boxDims[0]
01571                  : get_vertex_from_seq( i, j, k ) );
01572 }
01573 
01574 inline EntityHandle ScdBox::get_vertex( const HomCoord& ijk ) const
01575 {
01576     return get_vertex( ijk[0], ijk[1], ijk[2] );
01577 }
01578 
01579 inline bool ScdBox::contains( const HomCoord& ijk ) const
01580 {
01581     return ( ijk >= HomCoord( boxDims, 3 ) && ijk <= HomCoord( boxDims + 3, 3 ) );
01582 }
01583 
01584 inline bool ScdBox::contains( int i, int j, int k ) const
01585 {
01586     return contains( HomCoord( i, j, k ) );
01587 }
01588 
01589 inline void ScdBox::box_set( EntityHandle this_set )
01590 {
01591     boxSet = this_set;
01592 }
01593 
01594 inline EntityHandle ScdBox::box_set()
01595 {
01596     return boxSet;
01597 }
01598 
01599 inline ScdVertexData* ScdBox::vert_dat() const
01600 {
01601     return vertDat;
01602 }
01603 
01604 inline StructuredElementSeq* ScdBox::elem_seq() const
01605 {
01606     return elemSeq;
01607 }
01608 
01609 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const
01610 {
01611     HomCoord hc;
01612     ErrorCode rval = get_params( ent, hc );
01613     if( MB_SUCCESS == rval )
01614     {
01615         i   = hc[0];
01616         j   = hc[1];
01617         k   = hc[2];
01618         dir = hc[3];
01619     }
01620 
01621     return rval;
01622 }
01623 
01624 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k ) const
01625 {
01626     HomCoord hc;
01627     ErrorCode rval = get_params( ent, hc );
01628     if( MB_SUCCESS == rval )
01629     {
01630         i = hc[0];
01631         j = hc[1];
01632         k = hc[2];
01633     }
01634 
01635     return rval;
01636 }
01637 
01638 inline bool ScdBox::locally_periodic_i() const
01639 {
01640     return locallyPeriodic[0];
01641 }
01642 
01643 inline bool ScdBox::locally_periodic_j() const
01644 {
01645     return locallyPeriodic[1];
01646 }
01647 
01648 inline bool ScdBox::locally_periodic_k() const
01649 {
01650     return locallyPeriodic[2];
01651 }
01652 
01653 inline void ScdBox::locally_periodic( bool lperiodic[3] )
01654 {
01655     for( int i = 0; i < 3; i++ )
01656         locallyPeriodic[i] = lperiodic[i];
01657 }
01658 
01659 inline const int* ScdBox::locally_periodic() const
01660 {
01661     return locallyPeriodic;
01662 }
01663 
01664 std::ostream& operator<<( std::ostream& str, const ScdParData& pd );
01665 
01666 }  // namespace moab
01667 #endif
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines