Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
normals.cpp
Go to the documentation of this file.
00001 /*
00002  * normals.cpp
00003  *
00004  *  Created on: Oct 30, 2015
00005  *     will take an occ model and compute normals at all vertices, using exact geometry info
00006  */
00007 
00008 #include "iGeom.h"
00009 #include "iMesh.h"
00010 #include "iRel.h"
00011 
00012 #include <cstdio>
00013 #include <cstdlib>
00014 #include <cstring>
00015 
00016 #include "MBiMesh.hpp"
00017 #include "moab/Core.hpp"
00018 
00019 #define STRINGIFY_( X ) #X
00020 #define STRINGIFY( X )  STRINGIFY_( X )
00021 #ifdef MESHDIR
00022 #ifdef HAVE_OCC
00023 #define DEFAULT_GEOM STRINGIFY( MESHDIR / irel / cyl2.stp )
00024 #define DEFAULT_MESH STRINGIFY( MESHDIR / irel / cyl2.h5m )
00025 #else
00026 #define DEFAULT_GEOM STRINGIFY( MESHDIR / irel / brick.facet )
00027 #define DEFAULT_MESH STRINGIFY( MESHDIR / irel / brick.h5m )
00028 #endif
00029 #endif
00030 
00031 #define CHECK_SIZE_C( type, array, allocated_size, size )                 \
00032     if( NULL == *( array ) || *( allocated_size ) == 0 )                  \
00033     {                                                                     \
00034         *( array )          = (type*)malloc( sizeof( type ) * ( size ) ); \
00035         *( allocated_size ) = size;                                       \
00036     }                                                                     \
00037     else if( *( allocated_size ) < ( size ) )                             \
00038     {                                                                     \
00039         printf( "   Array passed in is non-zero but too short.\n" );      \
00040     }
00041 
00042 typedef void* iRel_EntityHandle;
00043 
00044 int print_geom_info( iGeom_Instance geom, iBase_EntityHandle gent )
00045 {
00046     /* print information about this entity */
00047     int ent_type;
00048     int result;
00049     const char* type_names[] = { "Vertex", "Edge", "Face", "Region" };
00050 
00051     iGeom_getEntType( geom, gent, &ent_type, &result );
00052 
00053     if( iBase_SUCCESS != result )
00054     {
00055         printf( "Trouble getting entity adjacencies or types." );
00056         return 0;
00057     }
00058 
00059     printf( "%s 0x%lx\n", type_names[ent_type], (unsigned long)gent );
00060 
00061     return 1;
00062 }
00063 
00064 int print_mesh_info( iMesh_Instance mesh, iBase_EntityHandle ment )
00065 {
00066     /* print information about this entity */
00067 
00068     /* get adjacencies first; assume not more than 50 */
00069     iBase_EntityHandle adj_ents[500], *adj_ents_ptr = adj_ents;
00070     int ent_types[500], *ent_types_ptr              = ent_types;
00071     int adj_ents_alloc = 500, adj_ents_size, ent_types_size, ent_types_allocated = 500;
00072     int result;
00073 
00074     iBase_TagHandle* ment_tags = NULL;
00075     int ment_tags_size, ment_tags_alloc;
00076 
00077     char** tag_names;
00078     int i;
00079 
00080     const char* type_names[] = { "Vertex", "Edge", "Face", "Region" };
00081 
00082     iMesh_getEntAdj( mesh, ment, iBase_ALL_TYPES, &adj_ents_ptr, &adj_ents_alloc, &adj_ents_size, &result );
00083 
00084     if( iBase_SUCCESS != result ) return 0;
00085 
00086     /* put this ent on the end, then get types */
00087     adj_ents[adj_ents_size] = ment;
00088     iMesh_getEntArrType( mesh, adj_ents, adj_ents_size + 1, &ent_types_ptr, &ent_types_allocated, &ent_types_size,
00089                          &result );
00090     if( iBase_SUCCESS != result )
00091     {
00092         printf( "Trouble getting entity adjacencies or types." );
00093         return 0;
00094     }
00095 
00096     /* get tags on ment */
00097     iMesh_getAllTags( mesh, ment, &ment_tags, &ment_tags_alloc, &ment_tags_size, &result );
00098 
00099     printf( "Trouble getting tags on an entity or their names." );
00100 
00101     /* while we're at it, get all the tag names */
00102 
00103     tag_names = (char**)malloc( ment_tags_size * sizeof( char* ) );
00104 
00105     for( i = 0; i < ment_tags_size; i++ )
00106     {
00107         tag_names[i] = (char*)malloc( 120 * sizeof( char ) );
00108         iMesh_getTagName( mesh, ment_tags[i], tag_names[i], &result, 120 );
00109     }
00110 
00111     /* now print the information */
00112     printf( "%s %ld:\n", type_names[ent_types[ent_types_size - 1]], (long)ment );
00113     printf( "Adjacencies:" );
00114     for( i = 0; i < adj_ents_size; i++ )
00115     {
00116         if( i > 0 ) printf( ", " );
00117         printf( "%s %ld", type_names[ent_types[i]], (long)adj_ents[i] );
00118     }
00119     printf( "\nTags: \n" );
00120     for( i = 0; i < ment_tags_size; i++ )
00121     {
00122         int tag_type;
00123 
00124         printf( "%s ", tag_names[i] );
00125         iMesh_getTagType( mesh, ment_tags[i], &tag_type, &result );
00126         if( iBase_SUCCESS != result )
00127             printf( "(trouble getting type...)\n" );
00128         else
00129         {
00130             char* dum_handle     = NULL;
00131             int dum_handle_alloc = 0, dum_handle_size = 0;
00132             int int_data;
00133             double dbl_data;
00134             iBase_EntityHandle eh_data;
00135 
00136             switch( tag_type )
00137             {
00138                 case iBase_INTEGER:
00139                     iMesh_getIntData( mesh, ment, ment_tags[i], &int_data, &result );
00140                     printf( "(Int value=%d)", int_data );
00141                     break;
00142                 case iBase_DOUBLE:
00143                     iMesh_getDblData( mesh, ment, ment_tags[i], &dbl_data, &result );
00144                     printf( "(Dbl value=%f)", dbl_data );
00145                     break;
00146                 case iBase_ENTITY_HANDLE:
00147                     iMesh_getEHData( mesh, ment, ment_tags[i], &eh_data, &result );
00148                     printf( "(EH value=%ld)", (long)eh_data );
00149                     break;
00150                 case iBase_BYTES:
00151                     iMesh_getData( mesh, ment, ment_tags[i], (void**)&dum_handle, &dum_handle_alloc, &dum_handle_size,
00152                                    &result );
00153                     if( NULL != dum_handle && dum_handle_size > 0 ) printf( "(Opaque value=%c)", dum_handle[0] );
00154                     break;
00155             }
00156         }
00157 
00158         printf( "\n" );
00159     }
00160     printf( "(end tags)\n\n" );
00161     free( ment_tags );
00162 
00163     return 1;
00164 }
00165 
00166 /*!
00167   @li Load a geom and a mesh file
00168 */
00169 int load_geom_mesh( const char* geom_filename, const char* mesh_filename, iGeom_Instance geom, iMesh_Instance mesh )
00170 {
00171     /* load a geom */
00172     int result;
00173     iGeom_load( geom, geom_filename, 0, &result, strlen( geom_filename ), 0 );
00174     if( iBase_SUCCESS != result )
00175     {
00176         printf( "ERROR : can not load a geometry\n" );
00177         return 0;
00178     }
00179 
00180     /* load a mesh */
00181     iMesh_load( mesh, 0, mesh_filename, 0, &result, strlen( mesh_filename ), 0 );
00182     if( iBase_SUCCESS != result )
00183     {
00184         printf( "ERROR : can not load a mesh\n" );
00185         return 0;
00186     }
00187 
00188     return 1;
00189 }
00190 
00191 /*!
00192 
00193   @li Create relation between geom and mesh
00194 */
00195 int create_relation( iRel_Instance assoc, iGeom_Instance geom, iMesh_Instance mesh, iRel_PairHandle* pair )
00196 {
00197     int result;
00198     iBase_Instance iface1, iface2;
00199     int type1, type2;
00200     int ent_or_set1, ent_or_set2;
00201     int status1, status2;
00202 
00203     iRel_PairHandle tmp_pair;
00204     iRel_PairHandle* pair_ptr = &tmp_pair;
00205     int pairs_alloc           = 1, pairs_size;
00206 
00207     /* create an relation, entity to set */
00208     iRel_createPair( assoc, geom, iRel_ENTITY, iRel_IGEOM_IFACE, iRel_ACTIVE, mesh, iRel_SET, iRel_IMESH_IFACE,
00209                      iRel_ACTIVE, pair, &result );
00210     if( iBase_SUCCESS != result )
00211     {
00212         printf( "Couldn't create a new relation.\n" );
00213         return 0;
00214     }
00215 
00216     iRel_getPairInfo( assoc, *pair, &iface1, &ent_or_set1, &type1, &status1, &iface2, &ent_or_set2, &type2, &status2,
00217                       &result );
00218     if( iBase_SUCCESS != result )
00219     {
00220         printf( "Couldn't retrieve relation info.\n" );
00221         return 0;
00222     }
00223     if( iface1 != geom || ent_or_set1 != iRel_ENTITY || type1 != iRel_IGEOM_IFACE || iface2 != mesh ||
00224         ent_or_set2 != iRel_SET || type2 != iRel_IMESH_IFACE )
00225     {
00226         printf( "Unexpected relation info returned.\n" );
00227         return 0;
00228     }
00229 
00230     iRel_findPairs( assoc, geom, &pair_ptr, &pairs_alloc, &pairs_size, &result );
00231     if( iBase_SUCCESS != result )
00232     {
00233         printf( "Couldn't find relation pair when querying geom.\n" );
00234         return 0;
00235     }
00236     if( pairs_size != 1 || tmp_pair != *pair )
00237     {
00238         printf( "Unexpected relation pairs returned when querying geom.\n" );
00239         return 0;
00240     }
00241 
00242     iRel_findPairs( assoc, mesh, &pair_ptr, &pairs_alloc, &pairs_size, &result );
00243     if( iBase_SUCCESS != result )
00244     {
00245         printf( "Couldn't find relation pair when querying mesh.\n" );
00246         return 0;
00247     }
00248     if( pairs_size != 1 || tmp_pair != *pair )
00249     {
00250         printf( "Unexpected relation pairs returned when querying mesh.\n" );
00251         return 0;
00252     }
00253 
00254     return 1;
00255 }
00256 
00257 /*!
00258   @li Check relation between geom and mesh
00259 */
00260 int relate_geom_mesh( iRel_Instance assoc, iGeom_Instance geom, iMesh_Instance mesh, iRel_PairHandle pair )
00261 {
00262     int result;
00263 
00264     iBase_EntityHandle* gentities = NULL;
00265     int gentities_size = 0, gentities_alloc = 0;
00266 
00267     iBase_EntitySetHandle* mentity_handles = NULL;
00268     int mentity_handles_size = 0, mentity_handles_alloc = 0;
00269 
00270     const char* dim_tag_name = "GEOM_DIMENSION";
00271     iBase_TagHandle dim_tag_mesh;
00272 
00273     iBase_EntitySetHandle* mentities_vec;
00274     int mentities_vec_size = 0;
00275     int i;
00276 
00277     iBase_EntitySetHandle* out_mentities = NULL;
00278     int out_mentities_size = 0, out_mentities_alloc = 0;
00279 
00280     iBase_EntitySetHandle* out_mentities2 = NULL;
00281     int out_mentities2_size = 0, out_mentities2_alloc = 0;
00282 
00283     iBase_EntityHandle* out_gentities = NULL;
00284     int out_gentities_size = 0, out_gentities_alloc = 0;
00285 
00286     /* relate geometry entities with coresponding mesh entity sets */
00287     iGeom_getEntities( geom, NULL, iBase_VERTEX, &gentities, &gentities_alloc, &gentities_size, &result );
00288     if( iBase_SUCCESS != result )
00289     {
00290         printf( "Failed to get gentities by type in relate_geom_mesh.\n" );
00291         return 0;
00292     }
00293 
00294     iRel_inferEntArrRelations( assoc, pair, gentities, gentities_size, 0, &result );
00295     if( iBase_SUCCESS != result )
00296     {
00297         printf( "Failed to relate geom entities in relate_geom_mesh.\n" );
00298         return 0;
00299     }
00300 
00301     /* relate coresponding mesh entity sets for geometry entities */
00302     /* get 1-dimensional mesh entitysets */
00303     iMesh_getEntSets( mesh, NULL, 1, &mentity_handles, &mentity_handles_alloc, &mentity_handles_size, &result );
00304     if( iBase_SUCCESS != result )
00305     {
00306         printf( "Problem to get all entity sets.\n" );
00307         return 0;
00308     }
00309 
00310     /* get geom dimension tags for mesh entitysets */
00311     iMesh_createTag( mesh, dim_tag_name, 1, iBase_INTEGER, &dim_tag_mesh, &result, 15 );
00312     if( iBase_SUCCESS != result && result != iBase_TAG_ALREADY_EXISTS )
00313     {
00314         printf( "Couldn't create geom dim tag for mesh entities.\n" );
00315         return 0;
00316     }
00317 
00318     /* get 1-dimensional mesh entitysets */
00319     mentities_vec = (iBase_EntitySetHandle*)malloc( mentity_handles_size * sizeof( iBase_EntitySetHandle ) );
00320     for( i = 0; i < mentity_handles_size; i++ )
00321     { /* test */
00322         int dim;
00323         iMesh_getEntSetIntData( mesh, mentity_handles[i], dim_tag_mesh, &dim, &result );
00324         if( iBase_SUCCESS != result ) continue;
00325 
00326         if( dim == 1 ) mentities_vec[mentities_vec_size++] = mentity_handles[i];
00327     }
00328 
00329     iRel_inferSetArrRelations( assoc, pair, mentities_vec, mentities_vec_size, 1, &result );
00330     if( iBase_SUCCESS != result )
00331     {
00332         printf( "Failed to relate mesh entities in relate_geom_mesh.\n" );
00333         return 0;
00334     }
00335 
00336     /* relate all geometry and mesh entities */
00337     iRel_inferAllRelations( assoc, pair, &result );
00338     if( iBase_SUCCESS != result )
00339     {
00340         printf( "Failed to relate all geom and mesh entities in relate_geom_mesh.\n" );
00341         return 0;
00342     }
00343 
00344     /* reset geom entities list and get all geom entities (prev
00345        only vertices) */
00346     free( gentities );
00347     gentities       = NULL;
00348     gentities_alloc = 0;
00349     iGeom_getEntities( geom, NULL, iBase_ALL_TYPES, &gentities, &gentities_alloc, &gentities_size, &result );
00350     if( iBase_SUCCESS != result )
00351     {
00352         printf( "Failed to get gentities by type in relate_geom_mesh.\n" );
00353         return 0;
00354     }
00355 
00356     /* get related mesh entity sets for geometry entities */
00357     iRel_getEntArrSetArrRelation( assoc, pair, gentities, gentities_size, 0, &out_mentities, &out_mentities_alloc,
00358                                   &out_mentities_size, &result );
00359     if( iBase_SUCCESS != result )
00360     {
00361         printf( "Failed to get geom entities in relate_geom_mesh.\n" );
00362         return 0;
00363     }
00364 
00365     if( out_mentities_size != gentities_size )
00366     {
00367         printf( "Number of input geom entities and output mesh entity sets should be same\n" );
00368         return 0;
00369     }
00370 
00371     /* now try deleting this relation */
00372     iRel_rmvEntArrRelation( assoc, pair, gentities, gentities_size, 0, &result );
00373     if( iBase_SUCCESS != result )
00374     {
00375         printf( "Failed to remove relation in relate_geom_mesh.\n" );
00376         return 0;
00377     }
00378 
00379     iRel_getEntArrSetArrRelation( assoc, pair, gentities, gentities_size, 0, &out_mentities2, &out_mentities2_alloc,
00380                                   &out_mentities2_size, &result );
00381     if( iBase_SUCCESS == result )
00382     {
00383         printf( "Shouldn't have gotten mesh sets in relate_geom_mesh.\n" );
00384         return 0;
00385     }
00386 
00387     /* restore the relation, since we need it later */
00388     iRel_setEntArrSetArrRelation( assoc, pair, gentities, gentities_size, out_mentities, out_mentities_size, &result );
00389     if( iBase_SUCCESS != result )
00390     {
00391         printf( "Failed to restore relation in relate_geom_mesh.\n" );
00392         return 0;
00393     }
00394 
00395     /* get related geometry entities for mesh entity sets */
00396     iRel_getSetArrEntArrRelation( assoc, pair, out_mentities, out_mentities_size, 1, &out_gentities,
00397                                   &out_gentities_alloc, &out_gentities_size, &result );
00398     if( iBase_SUCCESS != result )
00399     {
00400         printf( "Failed to get mesh entities in relate_geom_mesh.\n" );
00401         return 0;
00402     }
00403 
00404     if( out_mentities_size != out_gentities_size )
00405     {
00406         printf( "Number of input mesh entity sets and output geom entities should be same\n" );
00407         return 0;
00408     }
00409     free( mentity_handles );
00410     mentity_handles = NULL;
00411     free( gentities );
00412     gentities = NULL;
00413     free( mentity_handles );
00414     mentity_handles = NULL;
00415     free( out_mentities );
00416     out_mentities = NULL;
00417     free( mentities_vec );
00418     mentities_vec = NULL;
00419     free( out_gentities );
00420     out_gentities = NULL;
00421     return 1;
00422 }
00423 
00424 /*!
00425   @li compute normals onto the given geometry
00426 */
00427 int compute_normals( iRel_Instance assoc, iGeom_Instance geom, iMesh_Instance mesh, iRel_PairHandle pair )
00428 {
00429     int result;
00430     int i;
00431 
00432     iBase_EntityHandle* gentities = NULL;
00433     int gentities_size = 0, gentities_alloc = 0;
00434 
00435     iBase_EntitySetHandle* out_mentities = NULL;
00436     int out_mentities_size, out_mentities_alloc = 0;
00437 
00438     /* get all the face geo entities, and find relation to some mesh entity */
00439     iGeom_getEntities( geom, NULL, iBase_FACE, &gentities, &gentities_alloc, &gentities_size, &result );
00440     if( iBase_SUCCESS != result )
00441     {
00442         printf( "Problem getting all geom entities.\n" );
00443         return 0;
00444     }
00445     for( i = 0; i < gentities_size; i++ )
00446     {
00447         print_geom_info( geom, gentities[i] );
00448     }
00449 
00450     iRel_getEntArrSetArrRelation( assoc, pair, gentities, gentities_size, 0, &out_mentities, &out_mentities_alloc,
00451                                   &out_mentities_size, &result );
00452     /* might not all be */
00453     if( iBase_SUCCESS != result )
00454     {
00455         printf( "Failed to get mesh entities related to geom entities in compute_normals.\n" );
00456         return 0;
00457     }
00458 
00459     /* check that they're all non-null */
00460     if( out_mentities_size != gentities_size )
00461     {
00462         printf( "Number of mesh & related geom entities don't match.\n" );
00463         return 0;
00464     }
00465 
00466     /* check to make sure they're mesh sets; how to do that? */
00467     for( i = 0; i < out_mentities_size; i++ )
00468     {
00469         int is_list;
00470         iMesh_isList( mesh, (iBase_EntitySetHandle)out_mentities[i], &is_list, &result );
00471         if( iBase_SUCCESS != result )
00472         {
00473             printf( "Entity set returned from classification wasn't valid.\n" );
00474             return 0;
00475         }
00476     }
00477 
00478     moab::Interface* mb = reinterpret_cast< MBiMesh* >( mesh )->mbImpl;
00479     if( !mb ) return 0;
00480     moab::ErrorCode rval = MB_SUCCESS;
00481 
00482     moab::Tag idtag;
00483     rval = mb->tag_get_handle( "GLOBAL_ID", idtag );
00484     if( MB_SUCCESS != rval ) return 0;
00485 
00486     for( i = 0; i < gentities_size; i++ )
00487     {
00488 
00489         EntityHandle meshSet = (EntityHandle)out_mentities[i];
00490         Range faces;
00491         rval = mb->get_entities_by_dimension( meshSet, 2, faces );  // MB_CHK_ERR(rval);
00492 
00493         Range vertices;
00494         rval = mb->get_connectivity( faces, vertices );
00495 
00496         std::vector< double > coords;
00497         coords.resize( vertices.size() * 3 );
00498         rval = mb->get_coords( vertices, &coords[0] );
00499 
00500         int id;
00501         rval = mb->tag_get_data( idtag, &meshSet, 1, &id );
00502         printf( "surface %d  has %d faces and %d vertices \n", id, (int)faces.size(), (int)vertices.size() );
00503 
00504         std::vector< double > normals;
00505         normals.resize( 3 * vertices.size() );
00506 
00507         std::vector< double > on_coordinates;
00508         on_coordinates.resize( 3 * vertices.size() );
00509 
00510         for( int j = 0; j < (int)vertices.size(); j++ )
00511         {
00512             iGeom_getEntNrmlXYZ( geom, gentities[i], coords[3 * j], coords[3 * j + 1], coords[3 * j + 2],
00513                                  &normals[3 * j], &normals[3 * j + 1], &normals[3 * j + 2], &result );
00514             if( iBase_SUCCESS != result )
00515             {
00516                 printf( "Problem getting normals.\n" );
00517                 return 0;
00518             }
00519             EntityHandle vertex = vertices[j];
00520             rval                = mb->tag_get_data( idtag, &vertex, 1, &id );
00521             printf( "v: %4d  coor: %12.6f %12.6f %12.6f norm:  %12.6f %12.6f %12.6f \n", id, coords[3 * j],
00522                     coords[3 * j + 1], coords[3 * j + 2], normals[3 * j], normals[3 * j + 1], normals[3 * j + 2] );
00523             ;
00524         }
00525     }
00526 
00527     free( gentities );
00528     gentities = NULL;
00529     free( out_mentities );
00530     out_mentities = NULL;
00531     /* ok, we're done */
00532     return 1;
00533 }
00534 
00535 int main( int argc, char* argv[] )
00536 {
00537     /* Check command line arg */
00538     const char* geom_filename = DEFAULT_GEOM;
00539     const char* mesh_filename = DEFAULT_MESH;
00540 
00541     int result;
00542 
00543     iGeom_Instance geom;
00544     iMesh_Instance mesh;
00545     iRel_Instance assoc;
00546     iRel_PairHandle pair;
00547 
00548     printf( "Usage: %s %s\n", geom_filename, mesh_filename );
00549 
00550     if( argc == 2 && !strcmp( argv[1], "-h" ) )
00551     {
00552         printf( "Usage: %s <geom_filename> <mesh_filename>\n", argv[0] );
00553         return 1;
00554     }
00555     else if( argc == 2 )
00556     {
00557         geom_filename = argv[1];
00558         mesh_filename = argv[1];
00559     }
00560     else if( argc == 3 )
00561     {
00562         geom_filename = argv[1];
00563         mesh_filename = argv[2];
00564     }
00565 
00566     /* initialize the Geometry */
00567     iGeom_newGeom( 0, &geom, &result, 0 );
00568 
00569     /* initialize the Mesh */
00570     iMesh_newMesh( 0, &mesh, &result, 0 );
00571 
00572     /* initialize the Associate */
00573     iRel_create( 0, &assoc, &result, 0 );
00574 
00575     printf( "   load_geom_mesh: \n" );
00576     result = load_geom_mesh( geom_filename, mesh_filename, geom, mesh );
00577 
00578     printf( "   create_relation: \n" );
00579     result = create_relation( assoc, geom, mesh, &pair );
00580 
00581     printf( "   relate_geom_mesh: \n" );
00582     result = relate_geom_mesh( assoc, geom, mesh, pair );
00583 
00584     /* compute normals */
00585     printf( "   compute normals: \n" );
00586     result = compute_normals( assoc, geom, mesh, pair );
00587 
00588     /* summary */
00589 
00590     iRel_destroy( assoc, &result );
00591     iMesh_dtor( mesh, &result );
00592     iGeom_dtor( geom, &result );
00593 
00594     return 0;
00595 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines