 |
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
|
Mesh Oriented datABase
(version 5.5.1)
An array-based unstructured mesh library
|
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
107 class ScdParData
108 {
109 public:
110 ScdParData() : partMethod( NOPART ), pComm( NULL )
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
120 enum PartitionMethod
121 {
122 ALLJORKORI = 0,
123 ALLJKBAL,
124 SQIJ,
125 SQJK,
126 SQIJK,
127 TRIVIAL,
128 RCBZOLTAN,
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
136 int partMethod;
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
148 ParallelComm* pComm;
149 };
150
151 class MOAB_EXPORT ScdInterface
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
482 Interface* mbImpl;
483
484 //! whether we've searched the database for boxes yet
485 bool searchedBoxes;
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
491 Tag boxPeriodicTag;
492
493 //! tag representing box lower and upper corners
494 Tag boxDimsTag;
495
496 //! tag representing global lower and upper corners
497 Tag globalBoxDimsTag;
498
499 //! tag representing partition method
500 Tag partMethodTag;
501
502 //! tag pointing from set to ScdBox
503 Tag boxSetTag;
504 };
505
506 class MOAB_EXPORT ScdBox
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 */
749 ScdParData& par_data()
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
813 ScdInterface* scImpl;
814
815 //! entity set representing this box
816 EntityHandle boxSet;
817
818 //! vertex sequence this box represents, if there's only one, otherwise they're
819 //! retrieved from the element sequence
820 ScdVertexData* vertDat;
821
822 //! element sequence this box represents
823 StructuredElementSeq* elemSeq;
824
825 //! starting vertex handle for this box
826 EntityHandle startVertex;
827
828 //! starting element handle for this box
829 EntityHandle startElem;
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
838 ScdParData parData;
839
840 //! parameter extents
841 HomCoord boxSize;
842
843 //! convenience parameters, (boxSize[1]-1)*(boxSize[0]-1) and boxSize[0]-1
844 int boxSizeIJ;
845 int boxSizeIJM1;
846 int boxSizeIM1;
847 };
848
849 inline ErrorCode ScdInterface::compute_partition( int np,
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 {
859 case ScdParData::ALLJORKORI:
860 case -1:
861 rval = compute_partition_alljorkori( np, nr, par_data.gDims, par_data.gPeriodic, ldims, lperiodic, pdims );
862 break;
863 case ScdParData::ALLJKBAL:
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
883 inline ErrorCode ScdInterface::compute_partition_alljorkori( int np,
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
983 inline ErrorCode ScdInterface::compute_partition_alljkbal( int np,
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
1073 inline ErrorCode ScdInterface::compute_partition_sqij( int np,
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
1174 inline ErrorCode ScdInterface::compute_partition_sqjk( int np,
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
1261 inline ErrorCode ScdInterface::compute_partition_sqijk( int np,
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
1410 inline ErrorCode ScdInterface::get_neighbor( int np,
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 {
1428 case ScdParData::ALLJORKORI:
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
1447 inline ErrorCode ScdInterface::tag_shared_vertices( ParallelComm* pcomm, EntityHandle seth )
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
1468 inline ScdInterface* ScdBox::sc_impl() const
1469 {
1470 return scImpl;
1471 }
1472
1473 inline EntityHandle ScdBox::start_vertex() const
1474 {
1475 return startVertex;
1476 }
1477
1478 inline void ScdBox::start_vertex( EntityHandle startv )
1479 {
1480 startVertex = startv;
1481 }
1482
1483 inline EntityHandle ScdBox::start_element() const
1484 {
1485 return startElem;
1486 }
1487
1488 inline void ScdBox::start_element( EntityHandle starte )
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
1525 inline HomCoord ScdBox::box_min() const
1526 {
1527 return HomCoord( boxDims, 3 );
1528 }
1529
1530 inline HomCoord ScdBox::box_max() const
1531 {
1532 return HomCoord( boxDims + 3, 3 );
1533 }
1534
1535 inline HomCoord ScdBox::box_size() const
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
1594 inline EntityHandle ScdBox::box_set()
1595 {
1596 return boxSet;
1597 }
1598
1599 inline ScdVertexData* ScdBox::vert_dat() const
1600 {
1601 return vertDat;
1602 }
1603
1604 inline StructuredElementSeq* ScdBox::elem_seq() const
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
1.9.1.