Mesh Oriented datABase  (version 5.5.1)
An array-based unstructured mesh library
ScdInterface.hpp
Go to the documentation of this file.
1 /** \file ScdInterface.hpp
2  */
3 #ifndef SCD_INTERFACE
4 #define SCD_INTERFACE
5 
6 #include "moab/Interface.hpp"
7 #include "moab/HomXform.hpp"
8 
9 #include <iostream>
10 #include <vector>
11 #include <cassert>
12 
13 #include "moab/win32_config.h"
14 
15 namespace moab
16 {
17 
18 class StructuredElementSeq;
19 class EntitySequence;
20 class ScdVertexData;
21 class EntitySequence;
22 class ScdBox;
23 class ParallelComm;
24 
25 /** \class ScdInterface ScdInterface.hpp "moab/ScdInterface.hpp"
26  * \brief A structured mesh interface for MOAB-based data
27  *
28  * Structured mesh in MOAB is created and accessed through the ScdInterface and ScdBox classes.
29  *
30  * \section Construction Construction and Representation
31  * Structured mesh can be constructed in one of two ways. First, a rectangular block of mesh,
32  * both vertices and edges/quads/hexes, can be created in one shot, using the construct_box method.
33  * In this case, there are single sequences of vertices/entities. The second method for creating
34  * structured mesh is to create the structured blocks of vertices and elements separately. In
35  * this case, different blocks of elements can share blocks of vertices, and each block of
36  * elements has its own independent parametric space. The algorithms behind this representation
37  * are described in T. Tautges, "MOAB-SD: Integrated structured and unstructured mesh
38  * representation", Eng. w Comp, vol 20 no. 3.
39  *
40  * Structured mesh is represented in MOAB down at the element sequence level, which is something
41  * applications don't see. In addition, when structured mesh is created, entity sets are also
42  * created and tagged with information about the parametric space. In particular, the BOX_DIMS
43  * tag is used to indicate the lower and upper corners in parametric space (this
44  * tag is integer size 6). Structured mesh blocks are also available through ScdBox class objects
45  * returned by ScdInterface. These class objects should be treated only as references to the
46  * structured mesh blocks; that is, the structured mesh referenced by these objects is not deleted
47  * when the ScdBox instance is destroyed. Functions for destroying the actual mesh are available
48  * on this class, though.
49  *
50  * Structured mesh blocks are returned in the form of ScdBox class objects. Each ScdBox instance
51  * represents a rectangular block of vertices and possibly elements (edges, quads, or hexes). The
52  * edge/quad/hex entity handles for a ScdBox are guaranteed to be contiguous, starting at a starting
53  * value which is also available through the ScdBox class. However, vertex handles may or may not
54  * be contiguous, depending on the construction method. The start vertex handle is also available
55  * from the ScdBox class.
56  *
57  * \section Parameters Parametric Space
58  *
59  * Each structured box has a parametric (ijk) space, which can be queried through the ScdBox
60  * interface. For non-periodic boxes, the edge/quad/hex parameter bounds are one less in each
61  * dimension than that of the vertices, otherwise they are the same as the vertex parameter bounds.
62  * In a parallel representation, boxes are locally non-periodic by default, but global ids are
63  * assigned such that the last set of vertices in a periodic direction match those of the first set
64  * of vertices in that direction.
65  *
66  * Entity handles are allocated with the i parameter varying fastest, then j, then k.
67  *
68  * \section Per Periodic Meshes
69  * Boxes can be periodic in i, or j, or both i and j. If only i or j is periodic, the corresponding
70  * mesh is a strip or an annular cylinder; if both i and j are periodic, the corresponding mesh is
71  * an annular torus. A box cannot be periodic in all three parameters. If i and/or j is periodic,
72  * and assuming IMIN/JMIN is zero, the parameter extents in the/each periodic direction (IMAX/JMAX)
73  * for vertices and edges/faces/hexes are the same, and the vertices on the "top" end in the
74  * periodic direction are at parameter value IMIN/JMIN.
75  *
76  * \section Par Parallel Representation
77  *
78  * For parallel structured meshes, each local mesh (the mesh on a given process) looks like a
79  * non-periodic structured mesh, and there are both local and global parameters of the structured
80  * mesh. If the mesh is periodic in a given direction, the last process in the periodic direction
81  * has local IMAX/JMAX that is one greater than the global IMAX/JMAX.
82  *
83  *
84  * In parallel, the periodicity described in the previous paragraph is "local periodicity"; there is
85  * also the notion of global periodicity. For serial meshes, those concepts are the same. In
86  * parallel, a mesh can be locally non-periodic but globally periodic in a given direction. In that
87  * case, the local mesh is still non-periodic, i.e. the parametric extents for edges is one fewer
88  * than that of vertices in that direction. However, vertices are given global ids such that they
89  * match those of the parametric minimum in that direction. Geometric positions of the vertices at
90  * the high end should still be greater than the ones just below.
91  *
92  * \section Adjs Adjacent Entities
93  * This interface supports parametric access to intermediate-dimension entities, e.g. adjacent faces
94  * and edges in a 3d mesh. In this case, a direction parameter is added, to identify the parametric
95  * direction of the entities being requested. For example, to retrieve the faces adjacent to a hex
96  * with parameters ijk, in the i parametric direction, you would use the parameters ijk0. These
97  * intermediate entities are not stored in a structured representation, but their parametric
98  * positions can be evaluated based on their adjacencies to higher-dimensional entities. Thanks to
99  * Milad Fatenejad for the thinking behind this.
100  *
101  * \section Evaluation Evaluation
102  * The ScdBox class provides functions for evaluating the mesh based on the ijk parameter space.
103  * These functions are inlined where possible, for efficiency.
104  */
105 
106 //! struct for keeping parallel data in one place
108 {
109  public:
111  {
112  gDims[0] = gDims[1] = gDims[2] = gDims[3] = gDims[4] = gDims[5] = 0;
113  gPeriodic[0] = gPeriodic[1] = gPeriodic[2] = 0;
114  pDims[0] = pDims[1] = pDims[2] = 0;
115  }
116 
117  //! Partition method enumeration; these strategies are described in comments for
118  //! compute_partition_alljorkori, compute_partition_alljkbal, compute_partition_sqij,
119  //! compute_partition_sqjk, and compute_partition_sqijk
121  {
129  NOPART
130  };
131 
132  //! Partition method names
133  static MOAB_EXPORT const char* PartitionMethodNames[NOPART + 1];
134 
135  //! partition method used to partition global parametric space
137 
138  //! lower and upper corners of global box
139  int gDims[6];
140 
141  //! is globally periodic in i or j or k
142  int gPeriodic[3];
143 
144  //! number of procs in each direction
145  int pDims[3];
146 
147  //! parallel communicator object for this par scd mesh
149 };
150 
152 {
153  public:
154  friend class ScdBox;
155 
156  //! Constructor
157  /** Constructor; if find_boxes is true, this will search for entity sets marked as
158  * structured blocks, based on the BOX_DIMS tag. Structured mesh blocks will be stored
159  * in this interface class for future retrieval. Structured mesh blocks created through
160  * this interface will also be stored here.
161  * \param impl MOAB instance
162  * \param find_boxes If true, search all the entity sets, caching the structured mesh blocks
163  */
164  ScdInterface( Interface* impl, bool find_boxes = false );
165 
166  // Destructor
167  ~ScdInterface();
168 
169  //! Return the MOAB Interface instance *
170  Interface* impl() const;
171 
172  //! Construct new structured mesh box, including both vertices and elements
173  /** Parameter range
174  * for vertex box is [low-high], for elements is [low-high). Construct quads by passing
175  * in low[2] == high[2], and edges by passing in low[1] == high[1] and low[2] == high[2].
176  * The result is passed back in a ScdBox*, which is a *reference* to the box of structured mesh.
177  * That is, the actual mesh is retained in MOAB when the ScdBox is destroyed. To actually
178  * destroy the mesh, call the destroy_mesh function on the ScdBox object first, before
179  * destroying it. \param low Lower corner in parameter space \param high Higher corner in
180  * parameter space \param coords Coordinates of vertices, interleaved (xyzxyz...); if NULL,
181  * coords are set to parametric values \param num_coords Number of coordinate values \param
182  * new_box Reference to box of structured mesh \param lperiodic[3] If lperiodic[s] != 0,
183  * direction s is locally periodic \param par_data If non-NULL, this will get stored on the
184  * ScdBox once created, contains info about global parallel nature of ScdBox across procs \param
185  * assign_global_ids If true, assigns 1-based global ids to vertices using GLOBAL_ID_TAG_NAME
186  * \param resolve_shared_ents If != -1, resolves shared entities up to and including dimension
187  * equal to value
188  */
189  ErrorCode construct_box( HomCoord low,
190  HomCoord high,
191  const double* const coords,
192  unsigned int num_coords,
193  ScdBox*& new_box,
194  int* const lperiodic = NULL,
195  ScdParData* const par_data = NULL,
196  bool assign_global_ids = false,
197  int resolve_shared_ents = -1 );
198 
199  //! Create a structured sequence of vertices, quads, or hexes
200  /** Starting handle for the sequence is available from the returned ScdBox.
201  * If creating a structured quad or hex box, subsequent calls must be made to ScdBox::add_vbox,
202  * until all the vertices have been filled in for the box.
203  * \param low Lower corner of structured box
204  * \param high Higher corner of structured box
205  * \param type EntityType, one of MBVERTEX, MBEDGE, MBQUAD, MBHEX
206  * \param starting_id Requested start id of entities
207  * \param new_box Reference to the newly created box of entities
208  * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
209  * (s=[0,1,2])
210  */
211  ErrorCode create_scd_sequence( const HomCoord& low,
212  const HomCoord& high,
213  EntityType type,
214  int starting_id,
215  ScdBox*& new_box,
216  int* is_periodic = NULL );
217 
218  //! Return all the structured mesh blocks in this MOAB instance, as ScdBox objects
219  /** Return the structured blocks in this MOAB instance. If these were not searched for
220  * at instantiation time, then the search is done now.
221  * \param boxes Vector of ScdBox objects representing structured mesh blocks
222  */
223  ErrorCode find_boxes( std::vector< ScdBox* >& boxes );
224 
225  //! Return all the structured mesh blocks in this MOAB instance, as entity set handles
226  /** Return the structured blocks in this MOAB instance. If these were not searched for
227  * at instantiation time, then the search is done now.
228  * \param boxes Range of entity set objects representing structured mesh blocks
229  */
230  ErrorCode find_boxes( Range& boxes );
231 
232  //! Return all the structured mesh blocks known by ScdInterface (does not search)
233  /** Return the structured blocks in this ScdInterface instance. Does not search for new boxes,
234  * just returns the contents of the list.
235  * \param boxes Structured boxes
236  */
237  ErrorCode get_boxes( std::vector< ScdBox* >& boxes );
238 
239  //! Return the tag marking the lower and upper corners of boxes
240  /**
241  * \param create_if_missing If the tag does not yet exist, create it
242  */
243  Tag box_dims_tag( bool create_if_missing = true );
244 
245  //! Return the tag marking the global lower and upper corners of boxes
246  /**
247  * \param create_if_missing If the tag does not yet exist, create it
248  */
249  Tag global_box_dims_tag( bool create_if_missing = true );
250 
251  //! Return the tag marking the partitioning method used to partition the box in parallel
252  /**
253  * \param create_if_missing If the tag does not yet exist, create it
254  */
255  Tag part_method_tag( bool create_if_missing = true );
256 
257  //! Return the tag marking whether box is periodic in i and j
258  /**
259  * \param create_if_missing If the tag does not yet exist, create it
260  */
261  Tag box_periodic_tag( bool create_if_missing = true );
262 
263  //! Return the tag marking the ScdBox for a set
264  /**
265  * \param create_if_missing If the tag does not yet exist, create it
266  */
267  Tag box_set_tag( bool create_if_missing = true );
268 
269  //! Return the ScdBox corresponding to the entity set passed in
270  /** If the entity isn't a structured box set, NULL is returned.
271  * \param eh Entity whose box is being queried
272  */
273  ScdBox* get_scd_box( EntityHandle eh );
274 
275  //! Compute a partition of structured parameter space
276  /** Compute a partition of structured parameter space, based on data in the ScdParData
277  * passed in. Results are passed back in arguments, which application can set back into
278  * par_data argument if they so desire.
279  * \param np Number of processors
280  * \param nr Rank of this processor
281  * \param par_data ScdParData object that contains input global parameter space, desired
282  * partitioning method, and information about global periodicity.
283  * \param ldims Local parameters for grid
284  * \param lperiodic Whether or not a given dimension is locally periodic
285  * \param pdims Number of procs in i, j, k directions
286  */
287  static ErrorCode compute_partition( int np,
288  int nr,
289  const ScdParData& par_data,
290  int* ldims,
291  int* lperiodic = NULL,
292  int* pdims = NULL );
293 
294  //! Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
295  //! all 3 params
296  /** Get information about the neighbor in the dijk[] direction, where dijk can be -1 or 1 for
297  * all 3 params \param np (in) Total # procs \param nr Processor from which neighbor is
298  * requested \param spd (in) ScdParData containing part method, gdims, gperiodic data \param
299  * dijk(*) (in) Direction being queried, = +/-1 or 0 \param pto (out) Processor holding the
300  * destination part \param rdims(6) (out) Parametric min/max of destination part \param
301  * facedims(6) (out) Parametric min/max of interface between pfrom and pto; if at the max in a
302  * periodic direction, set to global min of that direction \param across_bdy(3) (out) If
303  * across_bdy[i] is -1(1), interface with pto is across periodic lower(upper) bdy in parameter
304  * i, 0 otherwise
305  */
306  static ErrorCode get_neighbor( int np,
307  int nr,
308  const ScdParData& spd,
309  const int* const dijk,
310  int& pto,
311  int* rdims,
312  int* facedims,
313  int* across_bdy );
314 
315  //! Tag vertices with sharing data for parallel representations
316  /** Given the ParallelComm object to use, tag the vertices shared with other processors
317  */
318  ErrorCode tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth );
319 
320  //! Tag vertices with sharing data for parallel representations
321  /** Given the ParallelComm object to use, tag the vertices shared with other processors
322  */
323  ErrorCode tag_shared_vertices( ParallelComm* pcomm, ScdBox* box );
324 
325  protected:
326  //! Remove the box from the list on ScdInterface
327  ErrorCode remove_box( ScdBox* box );
328 
329  //! Add the box to the list on ScdInterface
330  ErrorCode add_box( ScdBox* box );
331 
332  private:
333  //! Create an entity set for a box, and tag with the parameters
334  /** \param low Lower corner parameters for this box
335  * \param high Upper corner parameters for this box
336  * \param scd_set Entity set created
337  * \param is_periodic[3] If is_periodic[s] is non-zero, mesh should be periodic in direction s
338  * (s=[0,1,2])
339  */
340  ErrorCode create_box_set( const HomCoord& low,
341  const HomCoord& high,
342  EntityHandle& scd_set,
343  int* is_periodic = NULL );
344 
345  //! Compute a partition of structured parameter space
346  /** Partitions the structured parametric space by partitioning j, k, or i only.
347  * If j is greater than #procs, partition that, else k, else i.
348  * For description of arguments, see ScdInterface::compute_partition.
349  */
350  inline static ErrorCode compute_partition_alljorkori( int np,
351  int nr,
352  const int gijk[6],
353  const int* const gperiodic,
354  int* lijk,
355  int* lperiodic,
356  int* pijk );
357 
358  //! Compute a partition of structured parameter space
359  /** Partitions the structured parametric space by partitioning j, and possibly k,
360  * seeking square regions of jk space
361  * For description of arguments, see ScdInterface::compute_partition.
362  */
363  inline static ErrorCode compute_partition_alljkbal( int np,
364  int nr,
365  const int gijk[6],
366  const int* const gperiodic,
367  int* lijk,
368  int* lperiodic,
369  int* pijk );
370 
371  //! Compute a partition of structured parameter space
372  /** Partitions the structured parametric space by seeking square ij partitions
373  * For description of arguments, see ScdInterface::compute_partition.
374  */
375  inline static ErrorCode compute_partition_sqij( int np,
376  int nr,
377  const int gijk[6],
378  const int* const gperiodic,
379  int* lijk,
380  int* lperiodic,
381  int* pijk );
382 
383  //! Compute a partition of structured parameter space
384  /** Partitions the structured parametric space by seeking square jk partitions
385  * For description of arguments, see ScdInterface::compute_partition.
386  */
387  inline static ErrorCode compute_partition_sqjk( int np,
388  int nr,
389  const int gijk[6],
390  const int* const gperiodic,
391  int* lijk,
392  int* lperiodic,
393  int* pijk );
394 
395  //! Compute a partition of structured parameter space
396  /** Partitions the structured parametric space by seeking square ijk partitions
397  * For description of arguments, see ScdInterface::compute_partition.
398  */
399  inline static ErrorCode compute_partition_sqijk( int np,
400  int nr,
401  const int gijk[6],
402  const int* const gperiodic,
403  int* lijk,
404  int* lperiodic,
405  int* pijk );
406 
407  //! Get vertices shared with other processors
408  /** Shared vertices returned as indices into each proc's handle space
409  * \param box Box used to get parametric space info
410  * \param procs Procs this proc shares vertices with
411  * \param offsets Offsets into indices list for each proc
412  * \param shared_indices local/remote indices of shared vertices
413  */
414  static ErrorCode get_shared_vertices( ParallelComm* pcomm,
415  ScdBox* box,
416  std::vector< int >& procs,
417  std::vector< int >& offsets,
418  std::vector< int >& shared_indices );
419 
420  static ErrorCode get_indices( const int* const ldims,
421  const int* const rdims,
422  const int* const across_bdy,
423  int* face_dims,
424  std::vector< int >& shared_indices );
425 
426  static ErrorCode get_neighbor_alljorkori( int np,
427  int pfrom,
428  const int* const gdims,
429  const int* const gperiodic,
430  const int* const dijk,
431  int& pto,
432  int* rdims,
433  int* facedims,
434  int* across_bdy );
435 
436  static ErrorCode get_neighbor_alljkbal( int np,
437  int pfrom,
438  const int* const gdims,
439  const int* const gperiodic,
440  const int* const dijk,
441  int& pto,
442  int* rdims,
443  int* facedims,
444  int* across_bdy );
445 
446  static ErrorCode get_neighbor_sqij( int np,
447  int pfrom,
448  const int* const gdims,
449  const int* const gperiodic,
450  const int* const dijk,
451  int& pto,
452  int* rdims,
453  int* facedims,
454  int* across_bdy );
455 
456  static ErrorCode get_neighbor_sqjk( int np,
457  int pfrom,
458  const int* const gdims,
459  const int* const gperiodic,
460  const int* const dijk,
461  int& pto,
462  int* rdims,
463  int* facedims,
464  int* across_bdy );
465 
466  static ErrorCode get_neighbor_sqijk( int np,
467  int pfrom,
468  const int* const gdims,
469  const int* const gperiodic,
470  const int* const dijk,
471  int& pto,
472  int* rdims,
473  int* facedims,
474  int* across_bdy );
475 
476  static int gtol( const int* gijk, int i, int j, int k );
477 
478  //! assign global ids to vertices in this box
479  ErrorCode assign_global_ids( ScdBox* box );
480 
481  //! interface instance
483 
484  //! whether we've searched the database for boxes yet
486 
487  //! structured mesh blocks; stored as ScdBox objects, can get sets from those
488  std::vector< ScdBox* > scdBoxes;
489 
490  //! tag representing whether box is periodic in i and j
492 
493  //! tag representing box lower and upper corners
495 
496  //! tag representing global lower and upper corners
498 
499  //! tag representing partition method
501 
502  //! tag pointing from set to ScdBox
504 };
505 
507 {
508  friend class ScdInterface;
509 
510  public:
511  //! Destructor
512  ~ScdBox();
513 
514  //! Return the ScdInterface responsible for this box
515  inline ScdInterface* sc_impl() const;
516 
517  //! Add a vertex box to this box
518  /* Add a vertex box to the element sequence referenced by this box. The passed in vbox must
519  * be a vertex box, with parametric extents no larger than that of this box. This vbox is
520  * oriented to this box by matching parameters from1-from3 in vbox to to1-to3 in this box.
521  * If bb_input is true, only the part of the vertex sequence between bb_min and bb_max is
522  * referenced \param vbox The vertex box being added to this box \param from1 1st reference
523  * point on vbox \param to1 1st reference point on this box \param from2 2nd reference point on
524  * vbox \param to2 2nd reference point on this box \param from3 3rd reference point on vbox
525  * \param to3 3rd reference point on this box
526  * \param bb_input If true, subsequent parameters list extents of vbox referenced
527  * \param bb_min Lower corner of rectangle referenced
528  * \param bb_max Upper corner of rectangle referenced
529  */
530  ErrorCode add_vbox( ScdBox* vbox,
531  HomCoord from1,
532  HomCoord to1,
533  HomCoord from2,
534  HomCoord to2,
535  HomCoord from3,
536  HomCoord to3,
537  bool bb_input = false,
538  const HomCoord& bb_min = HomCoord::getUnitv( 0 ),
539  const HomCoord& bb_max = HomCoord::getUnitv( 0 ) );
540  // const HomCoord &bb_min = HomCoord::unitv[0],
541  // const HomCoord &bb_max = HomCoord::unitv[0]);
542 
543  //! Return whether this box has all its vertices defined
544  /** Tests whether vertex boxs added with add_vbox have completely defined the vertex parametric
545  * space for this box.
546  *
547  */
548  bool boundary_complete() const;
549 
550  //! Return highest topological dimension of box
551  inline int box_dimension() const;
552 
553  //! Starting vertex handle for this box
554  inline EntityHandle start_vertex() const;
555 
556  //! Starting entity handle for this box
557  /** If this is a vertex box, the start vertex handle is returned.
558  */
559  inline EntityHandle start_element() const;
560 
561  //! Return the number of elements in the box
562  /* Number of elements is (boxSize[0]-1)(boxSize[1]-1)(boxSize[2]-1)
563  */
564  inline int num_elements() const;
565 
566  //! Return the number of vertices in the box
567  /* Number of vertices is boxSize[0] * boxSize[1] * boxSize[2]
568  */
569  inline int num_vertices() const;
570 
571  //! Return the parametric coordinates for this box
572  /**
573  * \return IJK parameters of lower and upper corners
574  */
575  inline const int* box_dims() const;
576 
577  //! Return the lower corner parametric coordinates for this box
578  inline HomCoord box_min() const;
579 
580  //! Return the upper corner parametric coordinates for this box
581  inline HomCoord box_max() const;
582 
583  //! Return the parameter extents for this box
584  inline HomCoord box_size() const;
585 
586  //! Return the parametric extents for this box
587  /**
588  * \param ijk IJK extents of this box
589  */
590  inline void box_size( int* ijk ) const;
591 
592  //! Return the parametric extents for this box
593  /**
594  * \param i I extent of this box
595  * \param j J extent of this box
596  * \param k K extent of this box
597  */
598  void box_size( int& i, int& j, int& k ) const;
599 
600  //! Get the element at the specified coordinates
601  /**
602  * \param ijk Parametric coordinates being evaluated
603  */
604  EntityHandle get_element( const HomCoord& ijk ) const;
605 
606  //! Get the element at the specified coordinates
607  /**
608  * \param i Parametric coordinates being evaluated
609  * \param j Parametric coordinates being evaluated
610  * \param k Parametric coordinates being evaluated
611  */
612  EntityHandle get_element( int i, int j = 0, int k = 0 ) const;
613 
614  //! Get the vertex at the specified coordinates
615  /**
616  * \param ijk Parametric coordinates being evaluated
617  */
618  EntityHandle get_vertex( const HomCoord& ijk ) const;
619 
620  //! Get the vertex at the specified coordinates
621  /**
622  * \param i Parametric coordinates being evaluated
623  * \param j Parametric coordinates being evaluated
624  * \param k Parametric coordinates being evaluated
625  */
626  EntityHandle get_vertex( int i, int j = 0, int k = 0 ) const;
627 
628  //! Get parametric coordinates of the specified entity
629  /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
630  * in this ScdBox.
631  * \param ent Entity being queried
632  * \param i Parametric coordinates returned
633  * \param j Parametric coordinates returned
634  * \param k Parametric coordinates returned
635  * \param dir Parametric coordinate direction returned (in case of getting adjacent
636  * edges (2d, 3d) or faces (3d); not modified otherwise
637  */
638  ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const;
639 
640  //! Get parametric coordinates of the specified entity, intermediate entities not allowed (no
641  //! dir parameter)
642  /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
643  * in this ScdBox, or MB_FAILURE if the entity is an intermediate-dimension entity.
644  * \param ent Entity being queried
645  * \param i Parametric coordinates returned
646  * \param j Parametric coordinates returned
647  * \param k Parametric coordinates returned
648  */
649  ErrorCode get_params( EntityHandle ent, int& i, int& j, int& k ) const;
650 
651  //! Get parametric coordinates of the specified entity
652  /** This function returns MB_ENTITY_NOT_FOUND if the entity is not
653  * in this ScdBox.
654  * \param ent Entity being queried
655  * \param ijkd Parametric coordinates returned (including direction, in case of
656  * getting adjacent edges (2d, 3d) or faces (3d))
657  */
658  ErrorCode get_params( EntityHandle ent, HomCoord& ijkd ) const;
659 
660  /** \brief Get the adjacent edge or face at a parametric location
661  * This function gets the left (i=0), front (j=0), or bottom (k=0) edge or face for a parametric
662  * element. Left, front, or bottom is indicated by dir = 0, 1, or 2, resp. All edges and faces
663  * in a structured mesh block can be accessed using these parameters. \param dim Dimension of
664  * adjacent entity being requested \param i Parametric coordinates of cell being evaluated
665  * \param j Parametric coordinates of cell being evaluated
666  * \param k Parametric coordinates of cell being evaluated
667  * \param dir Direction (0, 1, or 2), for getting adjacent edges (2d, 3d) or faces (3d)
668  * \param ent Entity returned from this function
669  * \param create_if_missing If true, creates the entity if it doesn't already exist
670  */
671  ErrorCode get_adj_edge_or_face( int dim,
672  int i,
673  int j,
674  int k,
675  int dir,
676  EntityHandle& ent,
677  bool create_if_missing = true ) const;
678 
679  //! Return whether the box contains the parameters passed in
680  /**
681  * \param i Parametric coordinates being evaluated
682  * \param j Parametric coordinates being evaluated
683  * \param k Parametric coordinates being evaluated
684  */
685  bool contains( int i, int j, int k ) const;
686 
687  //! Return whether the box contains the parameters passed in
688  /**
689  * \param i Parametric coordinates being evaluated
690  * \param j Parametric coordinates being evaluated
691  * \param k Parametric coordinates being evaluated
692  */
693  bool contains( const HomCoord& ijk ) const;
694 
695  //! Set/Get the entity set representing the box
696  void box_set( EntityHandle this_set );
697  EntityHandle box_set();
698 
699  //! Get coordinate arrays for vertex coordinates for a structured block
700  /** Returns error if there isn't a single vertex sequence associated with this structured block
701  * \param xc X coordinate array pointer returned
702  * \param yc Y coordinate array pointer returned
703  * \param zc Z coordinate array pointer returned
704  */
705  ErrorCode get_coordinate_arrays( double*& xc, double*& yc, double*& zc );
706 
707  //! Get read-only coordinate arrays for vertex coordinates for a structured block
708  /** Returns error if there isn't a single vertex sequence associated with this structured block
709  * \param xc X coordinate array pointer returned
710  * \param yc Y coordinate array pointer returned
711  * \param zc Z coordinate array pointer returned
712  */
713  ErrorCode get_coordinate_arrays( const double*& xc, const double*& yc, const double*& zc ) const;
714 
715  //! Return whether box is locally periodic in i
716  /** Return whether box is locally periodic in i
717  * \return True if box is locally periodic in i direction
718  */
719  bool locally_periodic_i() const;
720 
721  //! Return whether box is locally periodic in j
722  /** Return whether box is locally periodic in j
723  * \return True if box is locally periodic in j direction
724  */
725  bool locally_periodic_j() const;
726 
727  //! Return whether box is locally periodic in k
728  /** Return whether box is locally periodic in k
729  * \return True if box is locally periodic in k direction
730  */
731  bool locally_periodic_k() const;
732 
733  //! Set local periodicity
734  /**
735  * \param lperiodic Vector of ijk periodicities to set this box to
736  */
737  void locally_periodic( bool lperiodic[3] );
738 
739  //! Get local periodicity
740  /**
741  * \return Vector of ijk periodicities for this box
742  */
743  const int* locally_periodic() const;
744 
745  //! Return parallel data
746  /** Return parallel data, if there is any
747  * \return par_data Parallel data set on this box
748  */
750  {
751  return parData;
752  }
753 
754  //! Return parallel data
755  /** Return parallel data, if there is any
756  * \return par_data Parallel data set on this box
757  */
758  const ScdParData& par_data() const
759  {
760  return parData;
761  }
762 
763  //! set parallel data
764  /** Set parallel data for this box
765  * \param par_data Parallel data to be set on this box
766  */
767  void par_data( const ScdParData& par_datap )
768  {
769  parData = par_datap;
770  }
771 
772  private:
773  //! Constructor
774  /** Create a structured box instance; this constructor is private because it should only be
775  * called from ScdInterface, a friend class. This constructor takes two sequences, one of which
776  * can be NULL. If both sequences come in non-NULL, the first should be a VertexSequence*
777  * corresponding to a structured vertex sequence and the second should be a
778  * StructuredElementSeq*. If the 2nd is NULL, the first can be either of those types. The
779  * other members of this class are taken from the sequences (e.g. parametric space) or the box
780  * set argument. Tags on the box set should be set from the caller. \param sc_impl A
781  * ScdInterface instance \param box_set Entity set representing this rectangle \param seq1 An
782  * EntitySequence (see ScdBox description) \param seq2 An EntitySequence (see ScdBox
783  * description), or NULL
784  */
785  ScdBox( ScdInterface* sc_impl, EntityHandle box_set, EntitySequence* seq1, EntitySequence* seq2 = NULL );
786 
787  //! function to get vertex handle directly from sequence
788  /** \param i Parameter being queried
789  * \param j Parameter being queried
790  * \param k Parameter being queried
791  */
792  EntityHandle get_vertex_from_seq( int i, int j, int k ) const;
793 
794  //! set the vertex sequence
795  ErrorCode vert_dat( ScdVertexData* vert_dat );
796 
797  //! get the vertex sequence
798  ScdVertexData* vert_dat() const;
799 
800  //! set the element sequence
801  ErrorCode elem_seq( EntitySequence* elem_seq );
802 
803  //! get the element sequence
804  StructuredElementSeq* elem_seq() const;
805 
806  //! Set the starting vertex handle for this box
807  void start_vertex( EntityHandle startv );
808 
809  //! Set the starting entity handle for this box
810  void start_element( EntityHandle starte );
811 
812  //! interface instance
814 
815  //! entity set representing this box
817 
818  //! vertex sequence this box represents, if there's only one, otherwise they're
819  //! retrieved from the element sequence
821 
822  //! element sequence this box represents
824 
825  //! starting vertex handle for this box
827 
828  //! starting element handle for this box
830 
831  //! lower and upper corners
832  int boxDims[6];
833 
834  //! is locally periodic in i or j or k
835  int locallyPeriodic[3];
836 
837  //! parallel data associated with this box, if any
839 
840  //! parameter extents
842 
843  //! convenience parameters, (boxSize[1]-1)*(boxSize[0]-1) and boxSize[0]-1
847 };
848 
850  int nr,
851  const ScdParData& par_data,
852  int* ldims,
853  int* lperiodic,
854  int* pdims )
855 {
856  ErrorCode rval = MB_SUCCESS;
857  switch( par_data.partMethod )
858  {
860  case -1:
861  rval = compute_partition_alljorkori( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
862  break;
864  rval = compute_partition_alljkbal( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
865  break;
866  case ScdParData::SQIJ:
867  rval = compute_partition_sqij( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
868  break;
869  case ScdParData::SQJK:
870  rval = compute_partition_sqjk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
871  break;
872  case ScdParData::SQIJK:
873  rval = compute_partition_sqijk( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
874  break;
875  default:
876  rval = MB_FAILURE;
877  break;
878  }
879 
880  return rval;
881 }
882 
884  int nr,
885  const int gijk[6],
886  const int* const gperiodic,
887  int* ldims,
888  int* lperiodic,
889  int* pijk )
890 {
891  // partition *the elements* over the parametric space; 1d partition for now, in the j, k, or i
892  // parameters
893  int tmp_lp[3], tmp_pijk[3];
894  if( !lperiodic ) lperiodic = tmp_lp;
895  if( !pijk ) pijk = tmp_pijk;
896 
897  for( int i = 0; i < 3; i++ )
898  lperiodic[i] = gperiodic[i];
899 
900  if( np == 1 )
901  {
902  if( ldims )
903  {
904  ldims[0] = gijk[0];
905  ldims[3] = gijk[3];
906  ldims[1] = gijk[1];
907  ldims[4] = gijk[4];
908  ldims[2] = gijk[2];
909  ldims[5] = gijk[5];
910  }
911  pijk[0] = pijk[1] = pijk[2] = 1;
912  }
913  else
914  {
915  if( gijk[4] - gijk[1] > np )
916  {
917  // partition j over procs
918  int dj = ( gijk[4] - gijk[1] ) / np;
919  int extra = ( gijk[4] - gijk[1] ) % np;
920  ldims[1] = gijk[1] + nr * dj + std::min( nr, extra );
921  ldims[4] = ldims[1] + dj + ( nr < extra ? 1 : 0 );
922 
923  if( gperiodic[1] && np > 1 )
924  {
925  lperiodic[1] = 0;
926  ldims[4]++;
927  }
928 
929  ldims[2] = gijk[2];
930  ldims[5] = gijk[5];
931  ldims[0] = gijk[0];
932  ldims[3] = gijk[3];
933  pijk[0] = pijk[2] = 1;
934  pijk[1] = np;
935  }
936  else if( gijk[5] - gijk[2] > np )
937  {
938  // partition k over procs
939  int dk = ( gijk[5] - gijk[2] ) / np;
940  int extra = ( gijk[5] - gijk[2] ) % np;
941  ldims[2] = gijk[2] + nr * dk + std::min( nr, extra );
942  ldims[5] = ldims[2] + dk + ( nr < extra ? 1 : 0 );
943 
944  ldims[1] = gijk[1];
945  ldims[4] = gijk[4];
946  ldims[0] = gijk[0];
947  ldims[3] = gijk[3];
948  pijk[0] = pijk[1] = 1;
949  pijk[2] = np;
950  }
951  else if( gijk[3] - gijk[0] > np )
952  {
953  // partition i over procs
954  int di = ( gijk[3] - gijk[0] ) / np;
955  int extra = ( gijk[3] - gijk[0] ) % np;
956  ldims[0] = gijk[0] + nr * di + std::min( nr, extra );
957  ldims[3] = ldims[0] + di + ( nr < extra ? 1 : 0 );
958 
959  if( gperiodic[0] && np > 1 )
960  {
961  lperiodic[0] = 0;
962  ldims[3]++;
963  }
964 
965  ldims[2] = gijk[2];
966  ldims[5] = gijk[5];
967  ldims[1] = gijk[1];
968  ldims[4] = gijk[4];
969 
970  pijk[1] = pijk[2] = 1;
971  pijk[0] = np;
972  }
973  else
974  {
975  // Couldn't find a suitable partition...
976  return MB_FAILURE;
977  }
978  }
979 
980  return MB_SUCCESS;
981 }
982 
984  int nr,
985  const int gijk[6],
986  const int* const gperiodic,
987  int* ldims,
988  int* lperiodic,
989  int* pijk )
990 {
991  int tmp_lp[3], tmp_pijk[3];
992  if( !lperiodic ) lperiodic = tmp_lp;
993  if( !pijk ) pijk = tmp_pijk;
994 
995  for( int i = 0; i < 3; i++ )
996  lperiodic[i] = gperiodic[i];
997 
998  if( np == 1 )
999  {
1000  if( ldims )
1001  {
1002  ldims[0] = gijk[0];
1003  ldims[3] = gijk[3];
1004  ldims[1] = gijk[1];
1005  ldims[4] = gijk[4];
1006  ldims[2] = gijk[2];
1007  ldims[5] = gijk[5];
1008  }
1009  pijk[0] = pijk[1] = pijk[2] = 1;
1010  }
1011  else
1012  {
1013  // improved, possibly 2-d partition
1014  std::vector< double > kfactors;
1015  kfactors.push_back( 1 );
1016  int K = gijk[5] - gijk[2];
1017  for( int i = 2; i < K; i++ )
1018  if( !( K % i ) && !( np % i ) ) kfactors.push_back( i );
1019  kfactors.push_back( K );
1020 
1021  // compute the ideal nj and nk
1022  int J = gijk[4] - gijk[1];
1023  double njideal = sqrt( ( (double)( np * J ) ) / ( (double)K ) );
1024  double nkideal = ( njideal * K ) / J;
1025 
1026  int nk, nj;
1027  if( nkideal < 1.0 )
1028  {
1029  nk = 1;
1030  nj = np;
1031  }
1032  else
1033  {
1034  std::vector< double >::iterator vit = std::lower_bound( kfactors.begin(), kfactors.end(), nkideal );
1035  if( vit == kfactors.begin() )
1036  nk = 1;
1037  else
1038  nk = (int)*( --vit );
1039  nj = np / nk;
1040  }
1041 
1042  int dk = K / nk;
1043  int dj = J / nj;
1044 
1045  ldims[2] = gijk[2] + ( nr % nk ) * dk;
1046  ldims[5] = ldims[2] + dk;
1047 
1048  int extra = J % nj;
1049 
1050  ldims[1] = gijk[1] + ( nr / nk ) * dj + std::min( nr / nk, extra );
1051  ldims[4] = ldims[1] + dj + ( nr / nk < extra ? 1 : 0 );
1052 
1053  ldims[0] = gijk[0];
1054  ldims[3] = gijk[3];
1055 
1056  if( gperiodic[1] && np > 1 )
1057  {
1058  lperiodic[1] = 0;
1059  if( nr / nk == nj - 1 )
1060  {
1061  ldims[1]++;
1062  }
1063  }
1064 
1065  pijk[0] = 1;
1066  pijk[1] = nj;
1067  pijk[2] = nk;
1068  }
1069 
1070  return MB_SUCCESS;
1071 }
1072 
1074  int nr,
1075  const int gijk[6],
1076  const int* const gperiodic,
1077  int* ldims,
1078  int* lperiodic,
1079  int* pijk )
1080 {
1081  int tmp_lp[3], tmp_pijk[3];
1082  if( !lperiodic ) lperiodic = tmp_lp;
1083  if( !pijk ) pijk = tmp_pijk;
1084 
1085  // square IxJ partition
1086 
1087  for( int i = 0; i < 3; i++ )
1088  lperiodic[i] = gperiodic[i];
1089 
1090  if( np == 1 )
1091  {
1092  if( ldims )
1093  {
1094  ldims[0] = gijk[0];
1095  ldims[3] = gijk[3];
1096  ldims[1] = gijk[1];
1097  ldims[4] = gijk[4];
1098  ldims[2] = gijk[2];
1099  ldims[5] = gijk[5];
1100  }
1101  pijk[0] = pijk[1] = pijk[2] = 1;
1102  }
1103  else
1104  {
1105  std::vector< double > pfactors, ppfactors;
1106  for( int i = 2; i <= np / 2; i++ )
1107  if( !( np % i ) )
1108  {
1109  pfactors.push_back( i );
1110  ppfactors.push_back( ( (double)( i * i ) ) / np );
1111  }
1112  pfactors.push_back( np );
1113  ppfactors.push_back( (double)np );
1114 
1115  // ideally, Px/Py = I/J
1116  double ijratio = ( (double)( gijk[3] - gijk[0] ) ) / ( (double)( gijk[4] - gijk[1] ) );
1117 
1118  unsigned int ind = 0;
1119  std::vector< double >::iterator optimal = std::lower_bound( ppfactors.begin(), ppfactors.end(), ijratio );
1120  if( optimal == ppfactors.end() )
1121  {
1122  ind = ppfactors.size() - 1;
1123  }
1124  else
1125  {
1126  ind = optimal - ppfactors.begin();
1127  if( ind && fabs( ppfactors[ind - 1] - ijratio ) < fabs( ppfactors[ind] - ijratio ) ) ind--;
1128  }
1129 
1130  // VARIABLES DESCRIBING THE MESH:
1131  // pi, pj = # procs in i and j directions
1132  // nri, nrj = my proc's position in i, j directions
1133  // I, J = # edges/elements in i, j directions
1134  // iextra, jextra = # procs having extra edge in i/j direction
1135  // top_i, top_j = if true, I'm the last proc in the i/j direction
1136  // i, j = # edges locally in i/j direction, *not* including one for iextra/jextra
1137  int pi = pfactors[ind];
1138  int pj = np / pi;
1139 
1140  int I = ( gijk[3] - gijk[0] ), J = ( gijk[4] - gijk[1] );
1141  int iextra = I % pi, jextra = J % pj, i = I / pi, j = J / pj;
1142  int nri = nr % pi, nrj = nr / pi;
1143 
1144  if( ldims )
1145  {
1146  ldims[0] = gijk[0] + i * nri + std::min( iextra, nri );
1147  ldims[3] = ldims[0] + i + ( nri < iextra ? 1 : 0 );
1148  ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
1149  ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
1150 
1151  ldims[2] = gijk[2];
1152  ldims[5] = gijk[5];
1153 
1154  if( gperiodic[0] && pi > 1 )
1155  {
1156  lperiodic[0] = 0;
1157  if( nri == pi - 1 ) ldims[3]++;
1158  }
1159  if( gperiodic[1] && pj > 1 )
1160  {
1161  lperiodic[1] = 0;
1162  if( nrj == pj - 1 ) ldims[4]++;
1163  }
1164  }
1165 
1166  pijk[0] = pi;
1167  pijk[1] = pj;
1168  pijk[2] = 1;
1169  }
1170 
1171  return MB_SUCCESS;
1172 }
1173 
1175  int nr,
1176  const int gijk[6],
1177  const int* const gperiodic,
1178  int* ldims,
1179  int* lperiodic,
1180  int* pijk )
1181 {
1182  int tmp_lp[3], tmp_pijk[3];
1183  if( !lperiodic ) lperiodic = tmp_lp;
1184  if( !pijk ) pijk = tmp_pijk;
1185 
1186  // square JxK partition
1187  for( int i = 0; i < 3; i++ )
1188  lperiodic[i] = gperiodic[i];
1189 
1190  if( np == 1 )
1191  {
1192  if( ldims )
1193  {
1194  ldims[0] = gijk[0];
1195  ldims[3] = gijk[3];
1196  ldims[1] = gijk[1];
1197  ldims[4] = gijk[4];
1198  ldims[2] = gijk[2];
1199  ldims[5] = gijk[5];
1200  }
1201  pijk[0] = pijk[1] = pijk[2] = 1;
1202  }
1203  else
1204  {
1205  std::vector< double > pfactors, ppfactors;
1206  for( int p = 2; p <= np; p++ )
1207  if( !( np % p ) )
1208  {
1209  pfactors.push_back( p );
1210  ppfactors.push_back( ( (double)( p * p ) ) / np );
1211  }
1212 
1213  // ideally, Pj/Pk = J/K
1214  int pj, pk;
1215  if( gijk[5] == gijk[2] )
1216  {
1217  pk = 1;
1218  pj = np;
1219  }
1220  else
1221  {
1222  double jkratio = ( (double)( gijk[4] - gijk[1] ) ) / ( (double)( gijk[5] - gijk[2] ) );
1223 
1224  std::vector< double >::iterator vit = std::lower_bound( ppfactors.begin(), ppfactors.end(), jkratio );
1225  if( vit == ppfactors.end() )
1226  --vit;
1227  else if( vit != ppfactors.begin() && fabs( *( vit - 1 ) - jkratio ) < fabs( ( *vit ) - jkratio ) )
1228  --vit;
1229  int ind = vit - ppfactors.begin();
1230 
1231  pj = 1;
1232  if( ind >= 0 && !pfactors.empty() ) pj = pfactors[ind];
1233  pk = np / pj;
1234  }
1235 
1236  int K = ( gijk[5] - gijk[2] ), J = ( gijk[4] - gijk[1] );
1237  int jextra = J % pj, kextra = K % pk, j = J / pj, k = K / pk;
1238  int nrj = nr % pj, nrk = nr / pj;
1239  ldims[1] = gijk[1] + j * nrj + std::min( jextra, nrj );
1240  ldims[4] = ldims[1] + j + ( nrj < jextra ? 1 : 0 );
1241  ldims[2] = gijk[2] + k * nrk + std::min( kextra, nrk );
1242  ldims[5] = ldims[2] + k + ( nrk < kextra ? 1 : 0 );
1243 
1244  ldims[0] = gijk[0];
1245  ldims[3] = gijk[3];
1246 
1247  if( gperiodic[1] && pj > 1 )
1248  {
1249  lperiodic[1] = 0;
1250  if( nrj == pj - 1 ) ldims[4]++;
1251  }
1252 
1253  pijk[0] = 1;
1254  pijk[1] = pj;
1255  pijk[2] = pk;
1256  }
1257 
1258  return MB_SUCCESS;
1259 }
1260 
1262  int nr,
1263  const int* const gijk,
1264  const int* const gperiodic,
1265  int* ldims,
1266  int* lperiodic,
1267  int* pijk )
1268 {
1269  if( gperiodic[0] || gperiodic[1] || gperiodic[2] ) return MB_FAILURE;
1270 
1271  int tmp_lp[3], tmp_pijk[3];
1272  if( !lperiodic ) lperiodic = tmp_lp;
1273  if( !pijk ) pijk = tmp_pijk;
1274 
1275  // square IxJxK partition
1276 
1277  for( int i = 0; i < 3; i++ )
1278  lperiodic[i] = gperiodic[i];
1279 
1280  if( np == 1 )
1281  {
1282  if( ldims )
1283  for( int i = 0; i < 6; i++ )
1284  ldims[i] = gijk[i];
1285  pijk[0] = pijk[1] = pijk[2] = 1;
1286  return MB_SUCCESS;
1287  }
1288 
1289  std::vector< int > pfactors;
1290  pfactors.push_back( 1 );
1291  for( int i = 2; i <= np / 2; i++ )
1292  if( !( np % i ) ) pfactors.push_back( i );
1293  pfactors.push_back( np );
1294 
1295  // test for IJ, JK, IK
1296  int IJK[3], dIJK[3];
1297  for( int i = 0; i < 3; i++ )
1298  IJK[i] = std::max( gijk[3 + i] - gijk[i], 1 );
1299  // order IJK from lo to hi
1300  int lo = 0, hi = 0;
1301  for( int i = 1; i < 3; i++ )
1302  {
1303  if( IJK[i] < IJK[lo] ) lo = i;
1304  if( IJK[i] > IJK[hi] ) hi = i;
1305  }
1306  if( lo == hi ) hi = ( lo + 1 ) % 3;
1307  int mid = 3 - lo - hi;
1308  // search for perfect subdivision of np that balances #cells
1309  int perfa_best = -1, perfb_best = -1;
1310  double ratio = 0.0;
1311  for( int po = 0; po < (int)pfactors.size(); po++ )
1312  {
1313  for( int pi = po; pi < (int)pfactors.size() && np / ( pfactors[po] * pfactors[pi] ) >= pfactors[pi]; pi++ )
1314  {
1315  int p3 =
1316  std::find( pfactors.begin(), pfactors.end(), np / ( pfactors[po] * pfactors[pi] ) ) - pfactors.begin();
1317  if( p3 == (int)pfactors.size() || pfactors[po] * pfactors[pi] * pfactors[p3] != np )
1318  continue; // po*pi should exactly factor np
1319  assert( po <= pi && pi <= p3 );
1320  // by definition, po <= pi <= p3
1321  double minl =
1322  std::min( std::min( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] ),
1323  maxl =
1324  std::max( std::max( IJK[lo] / pfactors[po], IJK[mid] / pfactors[pi] ), IJK[hi] / pfactors[p3] );
1325  if( minl / maxl > ratio )
1326  {
1327  ratio = minl / maxl;
1328  perfa_best = po;
1329  perfb_best = pi;
1330  }
1331  }
1332  }
1333  if( perfa_best == -1 || perfb_best == -1 ) return MB_FAILURE;
1334 
1335  // VARIABLES DESCRIBING THE MESH:
1336  // pijk[i] = # procs in direction i
1337  // numr[i] = my proc's position in direction i
1338  // dIJK[i] = # edges/elements in direction i
1339  // extra[i]= # procs having extra edge in direction i
1340  // top[i] = if true, I'm the last proc in direction i
1341 
1342  pijk[lo] = pfactors[perfa_best];
1343  pijk[mid] = pfactors[perfb_best];
1344  pijk[hi] = ( np / ( pfactors[perfa_best] * pfactors[perfb_best] ) );
1345  int extra[3] = { 0, 0, 0 }, numr[3];
1346  for( int i = 0; i < 3; i++ )
1347  {
1348  dIJK[i] = IJK[i] / pijk[i];
1349  extra[i] = IJK[i] % pijk[i];
1350  }
1351  numr[2] = nr / ( pijk[0] * pijk[1] );
1352  int rem = nr % ( pijk[0] * pijk[1] );
1353  numr[1] = rem / pijk[0];
1354  numr[0] = rem % pijk[0];
1355  for( int i = 0; i < 3; i++ )
1356  {
1357  extra[i] = IJK[i] % dIJK[i];
1358  ldims[i] = gijk[i] + numr[i] * dIJK[i] + std::min( extra[i], numr[i] );
1359  ldims[3 + i] = ldims[i] + dIJK[i] + ( numr[i] < extra[i] ? 1 : 0 );
1360  }
1361 
1362  return MB_SUCCESS;
1363 }
1364 
1365 inline int ScdInterface::gtol( const int* gijk, int i, int j, int k )
1366 {
1367  return ( ( k - gijk[2] ) * ( gijk[3] - gijk[0] + 1 ) * ( gijk[4] - gijk[1] + 1 ) +
1368  ( j - gijk[1] ) * ( gijk[3] - gijk[0] + 1 ) + i - gijk[0] );
1369 }
1370 
1371 inline ErrorCode ScdInterface::get_indices( const int* const ldims,
1372  const int* const rdims,
1373  const int* const across_bdy,
1374  int* face_dims,
1375  std::vector< int >& shared_indices )
1376 {
1377  // check for going across periodic bdy and face_dims not in my ldims (I'll always be on top in
1378  // that case)...
1379  if( across_bdy[0] > 0 && face_dims[0] != ldims[3] )
1380  face_dims[0] = face_dims[3] = ldims[3];
1381  else if( across_bdy[0] < 0 && face_dims[0] != ldims[0] )
1382  face_dims[0] = face_dims[3] = ldims[0];
1383  if( across_bdy[1] > 0 && face_dims[1] != ldims[4] )
1384  face_dims[1] = face_dims[4] = ldims[4];
1385  else if( across_bdy[1] < 0 && face_dims[1] != ldims[1] )
1386  face_dims[0] = face_dims[3] = ldims[1];
1387 
1388  for( int k = face_dims[2]; k <= face_dims[5]; k++ )
1389  for( int j = face_dims[1]; j <= face_dims[4]; j++ )
1390  for( int i = face_dims[0]; i <= face_dims[3]; i++ )
1391  shared_indices.push_back( gtol( ldims, i, j, k ) );
1392 
1393  if( across_bdy[0] > 0 && face_dims[0] != rdims[0] )
1394  face_dims[0] = face_dims[3] = rdims[0];
1395  else if( across_bdy[0] < 0 && face_dims[0] != rdims[3] )
1396  face_dims[0] = face_dims[3] = rdims[3];
1397  if( across_bdy[1] > 0 && face_dims[1] != rdims[1] )
1398  face_dims[1] = face_dims[4] = rdims[1];
1399  else if( across_bdy[1] < 0 && face_dims[1] != rdims[4] )
1400  face_dims[0] = face_dims[3] = rdims[4];
1401 
1402  for( int k = face_dims[2]; k <= face_dims[5]; k++ )
1403  for( int j = face_dims[1]; j <= face_dims[4]; j++ )
1404  for( int i = face_dims[0]; i <= face_dims[3]; i++ )
1405  shared_indices.push_back( gtol( rdims, i, j, k ) );
1406 
1407  return MB_SUCCESS;
1408 }
1409 
1411  int pfrom,
1412  const ScdParData& spd,
1413  const int* const dijk,
1414  int& pto,
1415  int* rdims,
1416  int* facedims,
1417  int* across_bdy )
1418 {
1419  if( !dijk[0] && !dijk[1] && !dijk[2] )
1420  {
1421  // not going anywhere, return
1422  pto = -1;
1423  return MB_SUCCESS;
1424  }
1425 
1426  switch( spd.partMethod )
1427  {
1429  case -1:
1430  return get_neighbor_alljorkori( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims,
1431  across_bdy );
1432  case ScdParData::ALLJKBAL:
1433  return get_neighbor_alljkbal( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
1434  case ScdParData::SQIJ:
1435  return get_neighbor_sqij( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
1436  case ScdParData::SQJK:
1437  return get_neighbor_sqjk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
1438  case ScdParData::SQIJK:
1439  return get_neighbor_sqijk( np, pfrom, spd.gDims, spd.gPeriodic, dijk, pto, rdims, facedims, across_bdy );
1440  default:
1441  break;
1442  }
1443 
1444  return MB_FAILURE;
1445 }
1446 
1448 {
1449  ScdBox* box = get_scd_box( seth );
1450  if( !box )
1451  {
1452  // look for contained boxes
1453  Range tmp_range;
1454  ErrorCode rval = mbImpl->get_entities_by_type( seth, MBENTITYSET, tmp_range );
1455  if( MB_SUCCESS != rval ) return rval;
1456  for( Range::iterator rit = tmp_range.begin(); rit != tmp_range.end(); ++rit )
1457  {
1458  box = get_scd_box( *rit );
1459  if( box ) break;
1460  }
1461  }
1462 
1463  if( !box ) return MB_FAILURE;
1464 
1465  return tag_shared_vertices( pcomm, box );
1466 }
1467 
1469 {
1470  return scImpl;
1471 }
1472 
1474 {
1475  return startVertex;
1476 }
1477 
1478 inline void ScdBox::start_vertex( EntityHandle startv )
1479 {
1480  startVertex = startv;
1481 }
1482 
1484 {
1485  return startElem;
1486 }
1487 
1489 {
1490  startElem = starte;
1491 }
1492 
1493 inline int ScdBox::num_elements() const
1494 {
1495  if( !startElem ) return 0; // not initialized yet
1496 
1497  /* for a structured mesh, total number of elements is obtained by multiplying
1498  number of elements in each direction
1499  number of elements in each direction is given by number of vertices in that direction minus 1
1500  if periodic in that direction, the last vertex is the same as first one, count one more
1501  element
1502  */
1503  int num_e_i = ( -1 == boxSize[0] || 1 == boxSize[0] ) ? 1 : boxSize[0] - 1;
1504  if( locallyPeriodic[0] ) ++num_e_i;
1505 
1506  int num_e_j = ( -1 == boxSize[1] || 1 == boxSize[1] ) ? 1 : boxSize[1] - 1;
1507  if( locallyPeriodic[1] ) ++num_e_j;
1508 
1509  int num_e_k = ( -1 == boxSize[2] || 1 == boxSize[2] ) ? 1 : boxSize[2] - 1;
1510  if( locallyPeriodic[2] ) ++num_e_k;
1511 
1512  return num_e_i * num_e_j * num_e_k;
1513 }
1514 
1515 inline int ScdBox::num_vertices() const
1516 {
1517  return boxSize[0] * ( !boxSize[1] ? 1 : boxSize[1] ) * ( !boxSize[2] ? 1 : boxSize[2] );
1518 }
1519 
1520 inline const int* ScdBox::box_dims() const
1521 {
1522  return boxDims;
1523 }
1524 
1526 {
1527  return HomCoord( boxDims, 3 );
1528 }
1529 
1531 {
1532  return HomCoord( boxDims + 3, 3 );
1533 }
1534 
1536 {
1537  return boxSize;
1538 }
1539 
1540 inline void ScdBox::box_size( int* ijk ) const
1541 {
1542  ijk[0] = boxSize[0];
1543  ijk[1] = boxSize[1];
1544  ijk[2] = boxSize[2];
1545 }
1546 
1547 inline void ScdBox::box_size( int& i, int& j, int& k ) const
1548 {
1549  i = boxSize[0];
1550  j = boxSize[1];
1551  k = boxSize[2];
1552 }
1553 
1554 inline EntityHandle ScdBox::get_element( int i, int j, int k ) const
1555 {
1556  return ( !startElem
1557  ? 0
1558  : startElem + ( k - boxDims[2] ) * boxSizeIJM1 + ( j - boxDims[1] ) * boxSizeIM1 + i - boxDims[0] );
1559 }
1560 
1561 inline EntityHandle ScdBox::get_element( const HomCoord& ijk ) const
1562 {
1563  return get_element( ijk[0], ijk[1], ijk[2] );
1564 }
1565 
1566 inline EntityHandle ScdBox::get_vertex( int i, int j, int k ) const
1567 {
1568  return ( vertDat
1569  ? startVertex + ( boxDims[2] == -1 && boxDims[5] == -1 ? 0 : ( k - boxDims[2] ) ) * boxSizeIJ +
1570  ( boxDims[1] == -1 && boxDims[4] == -1 ? 0 : ( j - boxDims[1] ) ) * boxSize[0] + i - boxDims[0]
1571  : get_vertex_from_seq( i, j, k ) );
1572 }
1573 
1574 inline EntityHandle ScdBox::get_vertex( const HomCoord& ijk ) const
1575 {
1576  return get_vertex( ijk[0], ijk[1], ijk[2] );
1577 }
1578 
1579 inline bool ScdBox::contains( const HomCoord& ijk ) const
1580 {
1581  return ( ijk >= HomCoord( boxDims, 3 ) && ijk <= HomCoord( boxDims + 3, 3 ) );
1582 }
1583 
1584 inline bool ScdBox::contains( int i, int j, int k ) const
1585 {
1586  return contains( HomCoord( i, j, k ) );
1587 }
1588 
1589 inline void ScdBox::box_set( EntityHandle this_set )
1590 {
1591  boxSet = this_set;
1592 }
1593 
1595 {
1596  return boxSet;
1597 }
1598 
1600 {
1601  return vertDat;
1602 }
1603 
1605 {
1606  return elemSeq;
1607 }
1608 
1609 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k, int& dir ) const
1610 {
1611  HomCoord hc;
1612  ErrorCode rval = get_params( ent, hc );
1613  if( MB_SUCCESS == rval )
1614  {
1615  i = hc[0];
1616  j = hc[1];
1617  k = hc[2];
1618  dir = hc[3];
1619  }
1620 
1621  return rval;
1622 }
1623 
1624 inline ErrorCode ScdBox::get_params( EntityHandle ent, int& i, int& j, int& k ) const
1625 {
1626  HomCoord hc;
1627  ErrorCode rval = get_params( ent, hc );
1628  if( MB_SUCCESS == rval )
1629  {
1630  i = hc[0];
1631  j = hc[1];
1632  k = hc[2];
1633  }
1634 
1635  return rval;
1636 }
1637 
1638 inline bool ScdBox::locally_periodic_i() const
1639 {
1640  return locallyPeriodic[0];
1641 }
1642 
1643 inline bool ScdBox::locally_periodic_j() const
1644 {
1645  return locallyPeriodic[1];
1646 }
1647 
1648 inline bool ScdBox::locally_periodic_k() const
1649 {
1650  return locallyPeriodic[2];
1651 }
1652 
1653 inline void ScdBox::locally_periodic( bool lperiodic[3] )
1654 {
1655  for( int i = 0; i < 3; i++ )
1656  locallyPeriodic[i] = lperiodic[i];
1657 }
1658 
1659 inline const int* ScdBox::locally_periodic() const
1660 {
1661  return locallyPeriodic;
1662 }
1663 
1664 std::ostream& operator<<( std::ostream& str, const ScdParData& pd );
1665 
1666 } // namespace moab
1667 #endif