itaps-parallel partial example
Mark Miller
miller86 at llnl.gov
Fri Dec 14 06:32:44 CST 2007
Hi All,
I have tried coding up a real, live working example of using iMesh interfaces in parallel to
do something like compute all parts in a partition are adjacent to a boundary condition.
But, don't get too excited!This didn't go the direction I first thought it would take me.
A picture of the mesh I was working from is attached. Each colored region in the picture
cooresponds to a domain. The end where z=0 has a boundary condition defined. It is
totally made up. Currently, there are no ghost zones. I hope to add that.
I started from a Silo file and wrote a converter to create a MOAB h5m file. The code that
did that -- its not relevant here but if you are interested i thought I would mention it -- is
silo_to_itaps.c. It creates an ent set for each domain (colored region) in the mesh. It
creates add'l ent sets for the nodes of each domain that are on the bc (z=0). It creates
2 tag handles to tag these sets and writes the data to a single h5m file. That file is
attached also.
The interesting code is par_imesh.c. In that file, I load the whole mesh (bogus for real
life but fine for this example) onto all procs. I get top-level ent sets and assign only those
having tag 'block' to processors. I also look for ent sets with tag 'bogus_bc' for the bc.
Some processors will have ent sets with this tag. On those processors, I ask iMesh to
find all ents adjacent to the nodes of the nodes in the bc that the given processor owns.
I have one final step to make this an interesting example and that is to report back
all those 'blocks' that are 'adjacent' to the bc as a whole. It turns out that the way I've
constructed the problem, I can basically already know that by asking which 'blocks'
have sub-ent sets that have tag 'bogus_bc'. But, I am trying to do it the hard way with
adj calls to iMesh.
I will work more over weekend but wanted to give you state of things now.
Karen, I realize this is far away from pseudocode but I honestly felt better trying to
build my way UP from iMesh than trying to work my way down from the still rather
abstract (and somewhat confusing to me) and still to be defined notions of iPart et. al.
If this is NOT at all going a direction that looks useful, let me know and I will
re-group. I admit that being a little overwhelmed earlier in the week with re-familiarizing
myself with iMesh and trying to build on top of that the iPart kinds of functionality.
I am pretty sure I can easily relax the load whole mesh onto all procs constraint before
end of weekend so it really would be realistic.
Mark
--
Mark C. Miller, Lawrence Livermore National Laboratory
email: mailto:miller86 at llnl.gov
(M/T/W) (925)-423-5901 (!!LLNL BUSINESS ONLY!!)
(Th/F) (530)-753-8511 (!!LLNL BUSINESS ONLY!!)
-------------- next part --------------
#include <mpi.h>
#include <iMesh.h>
#include <iBase.h>
#include <stdarg.h>
#include <stdio.h>
// end of list
#define EoL (void*)-1
// no list
#define NoL (0,EoL)
// supports passing variable length list of void pointers that if non-zero
// need to be free'd.
static void
ITAPSErrorCleanupHelper(int dummy, ...)
{
va_list ap;
va_start(ap, dummy);
void *ptr = va_arg(ap, void*);
while (ptr != EoL)
{
if (ptr)
free(ptr);
ptr = va_arg(ap, void*);
}
va_end(ap);
}
// ****************************************************************************
// Macro: CheckITAPSError2
//
// Purpose: Very useful macro for checking error condition of ITAPS
// implementation after a call into the implementation. This macro does work
// to get an error description, issues a VisIt warning, keeps track of how
// many times a given warning has been issued and suppresses above 5 aborts
// (via 'goto') to the end of the function from which it is called and
// ensures that all pointers allocated prior to the abort get freed.
//
//
// ****************************************************************************
#define CheckITAPSError2(IMI, ERR, ARGS, THELINE, THEFILE) \
if (ERR != 0) \
{ \
int i; \
char msg[1024]; \
char desc[256]; \
for (i = 0; i < sizeof(desc); i++) desc[i] = '\0'; \
int dummyError = ERR; \
iMesh_getDescription(IMI, desc, &dummyError, sizeof(desc)); \
snprintf(msg, sizeof(msg), "Encountered ITAPS error (%d) at line %d in file \"%s\"\n" \
"The description is...\n" \
" \"%s\"\n", ERR, THELINE, THEFILE, desc); \
fprintf(stderr, "%s\n", msg); \
ITAPSErrorCleanupHelper ARGS; \
goto funcEnd; \
}
#define CheckITAPSError(IMI, ARGS) CheckITAPSError2(IMI, iMerr, ARGS, __LINE__, __FILE__)
int main(int argc, char **argv)
{
int i, j, k;
char *imeshFileName = "./myTest.h5m";
int commSize = 1;
int commRank = 0;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &commSize);
MPI_Comm_rank(MPI_COMM_WORLD, &commRank);
iMesh_Instance iMi;
int iMerr;
iMesh_newMesh(0, &iMi, &iMerr, 0);
iBase_EntitySetHandle iMroot;
iMesh_getRootSet(iMi, &iMroot, &iMerr);
CheckITAPSError(iMi, NoL);
/* bogus: loads whole mesh onto every processor. Its just an example. */
iMesh_load(iMi,
iMroot, /*in const iBase_EntitySetHandle entity_set_handle,*/
imeshFileName, /*in const char *name,*/
0, /*in const char *options,*/
&iMerr, /*out int *err,*/
strlen(imeshFileName)+2,/*in int name_len,*/
0 /*in int options_len*/);
/* query some important tags */
iBase_TagHandle blockTag, bcTag;
iMesh_getTagHandle(iMi, "block", &blockTag, &iMerr, 7);
CheckITAPSError(iMi, NoL);
iMesh_getTagHandle(iMi, "bogus_bc", &bcTag, &iMerr, 10);
CheckITAPSError(iMi, NoL);
iBase_EntityHandle *top_sets = 0;
int top_sets_allocated = 0;
int top_sets_occupied = 0;
/* why is this returning 38 sets. I called createEntSet only 36 times */
/* must be universe and empty sets */
iMesh_getEntSets(iMi, iMroot, 1, &top_sets, &top_sets_allocated, &top_sets_occupied, &iMerr);
CheckITAPSError(iMi, NoL);
if (commRank == 0)
printf("top set count = %d\n", top_sets_occupied);
int nBlocks = 0;
iBase_EntityHandle block_sets[256];
for (k = 0; k < top_sets_occupied; k++)
{
int blockId;
iMesh_getIntData(iMi, top_sets[k], blockTag, &blockId, &iMerr);
if (iMerr == 0)
block_sets[nBlocks++] = top_sets[k];
else if (iMerr != iBase_TAG_NOT_FOUND)
CheckITAPSError(iMi, NoL);
}
if (commRank == 0)
printf("block count = %d\n", nBlocks);
int nBCsets = 0;
iBase_EntityHandle bc_sets[256];
for (k = 0; k < top_sets_occupied; k++)
{
int bcId;
iMesh_getIntData(iMi, top_sets[k], bcTag, &bcId, &iMerr);
if (iMerr == 0)
bc_sets[nBCsets++] = top_sets[k];
else if (iMerr != iBase_TAG_NOT_FOUND)
CheckITAPSError(iMi, NoL);
}
/* divide top-level blocks over processors */
int baseBlocksPerProc = nBlocks / commSize;
int remainderBlocks = nBlocks % commSize;
int myStart, myEnd;
if (commRank < remainderBlocks)
{
myStart = commRank * (baseBlocksPerProc+1);
myEnd = myStart + baseBlocksPerProc+1;
}
else
{
myStart = commRank * baseBlocksPerProc + remainderBlocks;
myEnd = myStart + baseBlocksPerProc;
}
if (nBCsets > 0)
{
for (k = 0; k < nBCsets; k++)
{
int bcId;
iMesh_getIntData(iMi, bc_sets[k], bcTag, &bcId, &iMerr);
if (myStart <= bcId && bcId < myEnd)
{
/* get hex ents adjacent to the bc nodes */
iBase_EntityHandle *ents_adj_to_bc_nodes = 0;
int adj_allocated = 0;
int adj_occupied = 0;
int *offsets = 0;
int offsets_allocated = 0;
int offsets_occupied = 0;
int *in_es = 0;
int in_es_alloc = 0;
int in_es_size = 0;
iMesh_getAdjEntities(iMi,
bc_sets[k], /*in const iBase_EntityHandle entity_set_handle,*/
iBase_VERTEX, /*in const int entity_type_requestor, */
iMesh_POINT, /*in const int entity_topology_requestor,*/
iBase_REGION, /*in const int entity_type_requested,*/
&ents_adj_to_bc_nodes, /*inout iBase_EntityHandle** adj_entity_handles,*/
&adj_allocated, /*inout int* adj_entity_handles_allocated,*/
&adj_occupied, /*out int* adj_entity_handles_size,*/
&offsets, /*inout int** offset,*/
&offsets_allocated, /*inout int* offset_allocated,*/
&offsets_occupied, /*out int* offset_size,*/
&in_es, /*inout int** in_entity_set,*/
&in_es_alloc, /*inout int* in_entity_set_allocated,*/
&in_es_size, /*out int* in_entity_set_size,*/
&iMerr);
CheckITAPSError(iMi, (0,ents_adj_to_bc_nodes,offsets,in_es,EoL));
printf("on processor %d, not %d ents adjacent to bc nodes for id %d\n",
commRank, adj_occupied, bcId);
}
}
}
#if 0
for (k = myStart; k < myEnd; k++)
{
int hex=-100;
int point=-200;
/* this returns zero for ALL_TOPOLOGIES */
iMesh_getNumOfTopo(iMi, block_sets[k], iMesh_HEXAHEDRON, &hex, &iMerr);
iMesh_getNumOfTopo(iMi, block_sets[k], iMesh_POINT, &point, &iMerr);
printf("processor %d reading block %d, has %d,%d points,hexes\n", commRank, k, point,hex);
CheckITAPSError(iMi, NoL);
}
#endif
funcEnd:
MPI_Finalize();
return 0;
}
-------------- next part --------------
#include <silo.h>
#include <iMesh.h>
#include <iBase.h>
#include <stdarg.h>
// end of list
#define EoL (void*)-1
// no list
#define NoL (0,EoL)
// supports passing variable length list of void pointers that if non-zero
// need to be free'd.
static void
ITAPSErrorCleanupHelper(int dummy, ...)
{
va_list ap;
va_start(ap, dummy);
void *ptr = va_arg(ap, void*);
while (ptr != EoL)
{
if (ptr)
free(ptr);
ptr = va_arg(ap, void*);
}
va_end(ap);
}
// ****************************************************************************
// Macro: CheckITAPSError2
//
// Purpose: Very useful macro for checking error condition of ITAPS
// implementation after a call into the implementation. This macro does work
// to get an error description, issues a VisIt warning, keeps track of how
// many times a given warning has been issued and suppresses above 5 aborts
// (via 'goto') to the end of the function from which it is called and
// ensures that all pointers allocated prior to the abort get freed.
//
//
// ****************************************************************************
#define CheckITAPSError2(IMI, ERR, ARGS, THELINE, THEFILE) \
if (ERR != 0) \
{ \
int i; \
char msg[1024]; \
char desc[256]; \
for (i = 0; i < sizeof(desc); i++) desc[i] = '\0'; \
int dummyError = ERR; \
iMesh_getDescription(IMI, desc, &dummyError, sizeof(desc)); \
snprintf(msg, sizeof(msg), "Encountered ITAPS error (%d) at line %d in file \"%s\"\n" \
"The description is...\n" \
" \"%s\"\n", ERR, THELINE, THEFILE, desc); \
fprintf(stderr, "%s\n", msg); \
ITAPSErrorCleanupHelper ARGS; \
goto funcEnd; \
}
#define CheckITAPSError(IMI, ARGS) CheckITAPSError2(IMI, iMerr, ARGS, __LINE__, __FILE__)
int main(int argc, char **argv)
{
int i, j, k;
char *siloFileName = "/usr/gapps/visit/data/multi_ucd3d.silo";
char *imeshFileName = "./myTest.h5m";
DBForceSingle(1);
DBfile *sfile = DBOpen(siloFileName, DB_UNKNOWN, DB_READ);
iMesh_Instance iMi;
int iMerr;
iMesh_newMesh(0, &iMi, &iMerr, 0);
iBase_EntitySetHandle iMroot;
iMesh_getRootSet(iMi, &iMroot, &iMerr);
CheckITAPSError(iMi, NoL);
iBase_TagHandle blockTag, bcTag;
iMesh_createTag(iMi, "block", 1, iBase_INTEGER, &blockTag, &iMerr, 7);
iMesh_createTag(iMi, "bogus_bc", 1, iBase_INTEGER, &bcTag, &iMerr, 10);
/* iterate over Silo multi-blocks */
DBmultimesh *mm = DBGetMultimesh(sfile, "mesh1");
for (k = 0; k < mm->nblocks; k++)
{
char tmpName[256];
/* ent/es creation always seems to occur in a global space and then
objects are assigned to sets */
iBase_EntityHandle blockes;
iMesh_createEntSet(iMi, 0, &blockes, &iMerr);
CheckITAPSError(iMi, NoL);
/* tag the block ent set with blockTag value of domain number */
iMesh_setEntSetIntData(iMi, blockes, blockTag, k, &iMerr);
CheckITAPSError(iMi, NoL);
/* add an es handle for this block to root */
iMesh_addEntSet(iMi, blockes, &iMroot, &iMerr);
sprintf(tmpName, "block%d/mesh1", k);
printf ("processing %s\n", tmpName);
/* get a ucd mesh from Silo */
DBucdmesh *um = DBGetUcdmesh(sfile, tmpName);
double *vtx_coords = (double *) malloc(3*um->nnodes*sizeof(double));
int *bc_node_ids = (int *) malloc(um->nnodes*sizeof(int));
int bc_nodes = 0;
j = 0;
for (i = 0; i < um->nnodes; i++)
{
vtx_coords[j++] = ((float*) um->coords[0])[i];
vtx_coords[j++] = ((float*) um->coords[1])[i];
vtx_coords[j++] = ((float*) um->coords[2])[i];
if ((um->coords[2])[i] <= 0.0)
bc_node_ids[bc_nodes++] = i;
}
/*iBase_EntityHandle *vtx_ents = (iBase_EntityHandle *) malloc(um->nnodes*sizeof(iBase_EntityHandle));*/
iBase_EntityHandle *vtx_ents = 0;
int vtx_allocated = 0;
int vtx_occupied = 0;
/* these seem to be created in some kind of a 'global' space.
How do I create vertex in an entity set */
iMesh_createVtxArr(iMi, um->nnodes, /*in const int num_verts, */
iBase_INTERLEAVED, /*in const int storage_order, */
vtx_coords, /*in const double* new_coords, */
3*um->nnodes, /*in const int new_coords_size, */
&vtx_ents, /*inout iBase_EntityHandle** new_vertex_handles,*/
&vtx_allocated, /*inout int* new_vertex_handles_allocated, */
&vtx_occupied, /*inout int* new_vertex_handles_size, */
&iMerr);
CheckITAPSError(iMi, (0,vtx_ents,EoL));
/* create refs to vtx ents for hexes */
iBase_EntityHandle *vtx_ent_refs = (iBase_EntityHandle *) malloc(um->zones->lnodelist*sizeof(iBase_EntityHandle));
for (i = 0; i < um->zones->lnodelist; i++)
vtx_ent_refs[i] = vtx_ents[um->zones->nodelist[i]];
/* create all the hexes for this block */
iBase_EntityHandle *hex_ents = 0;
int hex_allocated = 0;
int hex_occupied = 0;
int *hex_status = 0;
int hex_status_alloc = 0;
int hex_status_occ = 0;
iMesh_createEntArr(iMi,
iMesh_HEXAHEDRON, /*in const int new_entity_topology,*/
vtx_ent_refs, /*in const iBase_EntityHandle* lower_order_entity_handles,*/
um->zones->lnodelist, /*in const int lower_order_entity_handles_size,*/
&hex_ents, /*out iBase_EntityHandle** new_entity_handles,*/
&hex_allocated, /*out int* new_entity_handles_allocated,*/
&hex_occupied, /*out int* new_entity_handles_size,*/
&hex_status, /*inout int** status,*/
&hex_status_alloc, /*inout int* status_allocated,*/
&hex_status_occ, /*out int* status_size,*/
&iMerr);
CheckITAPSError(iMi, (0,vtx_ent_refs,hex_ents,hex_status,EoL));
/* now, add the hexes to the blockes */
iMesh_addEntArrToSet(iMi, hex_ents, hex_occupied, &blockes, &iMerr);
CheckITAPSError(iMi, NoL);
/* create entity set for any bc nodes we might have */
if (bc_nodes > 0)
{
iBase_EntityHandle bc_nodes_es;
iMesh_createEntSet(iMi, 0, &bc_nodes_es, &iMerr);
CheckITAPSError(iMi, NoL);
/* tag the bc_nodes_es ent set with bcTag value of block number */
iMesh_setEntSetIntData(iMi, bc_nodes_es, bcTag, k, &iMerr);
CheckITAPSError(iMi, NoL);
/* add an es handle for these bc nodes to blockes */
iMesh_addEntSet(iMi, bc_nodes_es, &blockes, &iMerr);
CheckITAPSError(iMi, NoL);
/* add the node entities to this es */
iBase_EntityHandle *bc_node_ents = (iBase_EntityHandle *) malloc(bc_nodes*sizeof(iBase_EntityHandle));
for (i = 0; i < bc_nodes; i++)
bc_node_ents[i] = vtx_ents[bc_node_ids[i]];
iMesh_addEntArrToSet(iMi, bc_node_ents, bc_nodes, &bc_nodes_es, &iMerr);
CheckITAPSError(iMi, (0,bc_node_ents,EoL));
free(bc_node_ids);
}
free(vtx_ent_refs);
free(hex_ents);
free(hex_status);
free(vtx_coords);
DBFreeUcdmesh(um);
}
iMesh_save(iMi,
iMroot, /*in const iBase_EntitySetHandle entity_set_handle, */
imeshFileName, /*in const char *name,*/
0, /*in const char *options,*/
&iMerr, /*out int *err,*/
strlen(imeshFileName)+2, /*in const int name_len, */
0 /*in int options_len*/);
CheckITAPSError(iMi, NoL);
funcEnd:
return 0;
}
-------------- next part --------------
A non-text attachment was scrubbed...
Name: myTest.h5m
Type: application/octet-stream
Size: 3371504 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/itaps-parallel/attachments/20071214/adbe0553/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: visit0000.png
Type: image/png
Size: 83041 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/itaps-parallel/attachments/20071214/adbe0553/attachment.png>
More information about the itaps-parallel
mailing list