MOAB: Mesh Oriented datABase
(version 5.3.1)
|
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