moab-dev

Download

[email protected]

July 2014

13 participants
72 discussions

commit/MOAB: 4 new changesets
by commits-noreply＠bitbucket.org 30 Jun '14

30 Jun '14

4 new commits in MOAB: https://bitbucket.org/fathomteam/moab/commits/7d12eea5999e/ Changeset: 7d12eea5999e Branch: None User: iulian07 Date: 2014-06-29 04:26:05 Summary: add large mesh example Affected #: 2 files diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp new file mode 100644 index 0000000..6ef1dc9 --- /dev/null +++ b/examples/GenLargeMesh.cpp @@ -0,0 +1,105 @@ +/** @example GenLargeMesh.cpp \n + * \brief Create a large structured mesh, partitioned \n + * To run: mpiexec -np 2 GenLargeMesh \n + * + * It shows how to load create a mesh on the fly, on multiple + * processors, block wise + * Each processor will create its version of a block mesh, partitioned + * as AxBxC blocks. Each block will be with blockSize^3 hexahedrons, and will get a + * different PARALLEL_PARTITION tag + * + * The number of tasks will be MxNxK, and it must match the mpi size + * Each task will generate its mesh at location (m,n,k) + * + * By default M=2, N=1, K=1, so by default it should be launched on 2 procs + * By default, blockSize is 4, and A=2, B=2, C=2, so each task will generate locally + * blockSize^3 x A x B x C hexahedrons (value = 64x8 = 512 hexas, in 8 partitions) + * + * The total number f partitions will be A*B*C*M*N*K (default 16) + * + * Each part in partition will get a proper tag + * + * The vertices will get a proper global id, which will be used to resolve the + * shared entities + * The output will be written in parallel, and we will try sizes as big as we can + * (up to a billion vertices, because we use int for global ids) + * + * Within each partition, the hexas will be numbered contiguously, and also the + * vertices; The global id will be determined easily, for a vertex, but the entity + * handle space will be more interesting to handle, within a partition (we want + * contiguous handles within a partition) + * + * We may or may not use ScdInterface, because we want control over global id and + * entity handle within a partition + * + * + * mpiexec -np 2 GenLargeMesh + */ + +#include "moab/Core.hpp" +#include "moab/ProgOptions.hpp" +#include "moab/ParallelComm.hpp" +#include "moab/CN.hpp" + +#include <iostream> +#include <vector> + +using namespace moab; +using namespace std; +int main(int argc, char **argv) +{ + int A=2, B=2, C=2, M=2, N=1, K=1; + int blockSize = 4; + + MPI_Init(&argc, &argv); + + ProgOptions opts; + + opts.addOpt<int>(std::string("blockSize,b"), + std::string("Block size of mesh (default=4)"), &blockSize); + opts.addOpt<int>(std::string("xproc,M"), + std::string("Number of processors in x dir (default=2)"), &M); + opts.addOpt<int>(std::string("yproc,N"), + std::string("Number of processors in y dir (default=1)"), &N); + opts.addOpt<int>(std::string("zproc,K"), + std::string("Number of processors in z dir (default=1)"), &K); + + opts.addOpt<int>(std::string("xblocks,A"), + std::string("Number of blocks on a task in x dir (default=2)"), &A); + opts.addOpt<int>(std::string("yblocks,B"), + std::string("Number of blocks on a task in y dir (default=2)"), &B); + opts.addOpt<int>(std::string("zblocks,C"), + std::string("Number of blocks on a task in x dir (default=2)"), &C); + + + opts.parseCommandLine(argc, argv); + + Interface* mb = new Core; + if (NULL == mb) + { + MPI_Finalize(); + return 1; + } + + MPI_Comm comm ; + int rank, size; + MPI_Comm_rank( MPI_COMM_WORLD, &rank ); + MPI_Comm_rank( MPI_COMM_WORLD, &size ); + + if (M*N*K != size) + { + if (rank==0) + { + cout <<"M*N*K=" << M*N*K << " != size = "<< " size\n"; + } + MPI_Finalize(); + return 1; + } + + // the global id of each vertex will be, on task + + MPI_Finalize(); + return 0; +} + + diff --git a/examples/makefile b/examples/makefile index 787c8d3..f41def7 100644 --- a/examples/makefile +++ b/examples/makefile @@ -71,6 +71,9 @@ VisTags: VisTags.o ${MOAB_LIBDIR}/libMOAB.la ReadWriteTest: ReadWriteTest.o ${MOAB_LIBDIR}/libMOAB.la ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} +GenLargeMesh: GenLargeMesh.o ${MOAB_LIBDIR}/libMOAB.la + ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} + clean: rm -rf *.o *.mod *.h5m ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES} ${F90EXAMPLES} https://bitbucket.org/fathomteam/moab/commits/89f888c9c193/ Changeset: 89f888c9c193 Branch: None User: iulian07 Date: 2014-06-30 05:49:04 Summary: is_valid should be on the interface too not only on the Core.hpp How do we know if an entity handle is valid or not? is there another way? Affected #: 1 file diff --git a/src/moab/Interface.hpp b/src/moab/Interface.hpp index 09f4219..c8102b8 100644 --- a/src/moab/Interface.hpp +++ b/src/moab/Interface.hpp @@ -1981,6 +1981,8 @@ public: SetIterator *&set_iter) = 0; /**@}*/ + + virtual bool is_valid(EntityHandle entity) const = 0; }; //! predicate for STL algorithms. Returns true if the entity handle is https://bitbucket.org/fathomteam/moab/commits/e861c1df4f15/ Changeset: e861c1df4f15 Branch: None User: iulian07 Date: 2014-06-30 05:51:54 Summary: validate sharedEnts it is becoming important if some entities are deleted after they are resolved they are part of the sharedEnts vector, but they are not deleted from there. This is the advantage of keeping entities in lists in member data we should hunt them down :( Affected #: 2 files diff --git a/src/parallel/ParallelComm.cpp b/src/parallel/ParallelComm.cpp index 10ce709..a8571af 100644 --- a/src/parallel/ParallelComm.cpp +++ b/src/parallel/ParallelComm.cpp @@ -8789,7 +8789,22 @@ ErrorCode ParallelComm::settle_intersection_points(Range & edges, Range & shared return MB_SUCCESS; // end copy +} +ErrorCode ParallelComm::correct_shared_entities() +{ + + std::vector<EntityHandle> good_ents; + for (size_t i=0; i<sharedEnts.size(); i++) + { + if (mbImpl->is_valid(sharedEnts[i])) + { + good_ents.push_back(sharedEnts[i]); + } } + sharedEnts = good_ents; + + return MB_SUCCESS; +} void ParallelComm::print_pstatus(unsigned char pstat, std::string &ostr) { diff --git a/src/parallel/moab/ParallelComm.hpp b/src/parallel/moab/ParallelComm.hpp index 3424c51..2a16ac2 100644 --- a/src/parallel/moab/ParallelComm.hpp +++ b/src/parallel/moab/ParallelComm.hpp @@ -932,6 +932,12 @@ namespace moab { ErrorCode settle_intersection_points(Range & edges, Range & shared_edges_owned, std::vector<std::vector<EntityHandle> *> & extraNodesVec, double tolerance); + /* \brief check the shared entities if they exist anymore + * will check the shared ents array, and clean it if necessary + * + */ + ErrorCode correct_shared_entities(); + private: ErrorCode reduce_void(int tag_data_type, const MPI_Op mpi_op, int num_ents, void *old_vals, void *new_vals); https://bitbucket.org/fathomteam/moab/commits/11db00407a80/ Changeset: 11db00407a80 Branch: iulian07/largemesh User: iulian07 Date: 2014-06-30 05:54:01 Summary: improve the large mesh example this example pushes the limits of moab It will be used later to create meshes larger than 1 billion vertices Each processor creates a structure of AxBxC blocks, each block with blockSize ^ 3 hexas The example should be run on MxNxK array of processors. After resolving and sharing the vertices, mesh is written in parallel. each block will represent a partition It is different than test paralle_write_test because there the number of partitions created is equal to the number of processors. Here, each processor creates AxBxC partitions For a total of AxBxC x MxNxK partitions block size, A, B, .., K are options, with default values The total number of hexas will be AxBxC x MxNxK x block^3 edges and faces are removed after resolve sharing the file will have just vertices and hexas, and it will be already partitioned We will add some tags of interest, to make it better. The writer failed if I just deleted the edges and faces. I had to introduce a method to correct sharedEnts in PC Affected #: 1 file diff --git a/examples/GenLargeMesh.cpp b/examples/GenLargeMesh.cpp index 6ef1dc9..1e0fa73 100644 --- a/examples/GenLargeMesh.cpp +++ b/examples/GenLargeMesh.cpp @@ -40,12 +40,16 @@ #include "moab/ProgOptions.hpp" #include "moab/ParallelComm.hpp" #include "moab/CN.hpp" +#include "moab/ReadUtilIface.hpp" +#include "moab/MergeMesh.hpp" #include <iostream> #include <vector> using namespace moab; using namespace std; + +#define CHECKE(message) if(rval != MB_SUCCESS) { cout << (message) << "\n"; MPI_Finalize(); return 1; } int main(int argc, char **argv) { int A=2, B=2, C=2, M=2, N=1, K=1; @@ -58,7 +62,7 @@ int main(int argc, char **argv) opts.addOpt<int>(std::string("blockSize,b"), std::string("Block size of mesh (default=4)"), &blockSize); opts.addOpt<int>(std::string("xproc,M"), - std::string("Number of processors in x dir (default=2)"), &M); + std::string("Number of processors in x dir (default=1)"), &M); opts.addOpt<int>(std::string("yproc,N"), std::string("Number of processors in y dir (default=1)"), &N); opts.addOpt<int>(std::string("zproc,K"), @@ -80,23 +84,175 @@ int main(int argc, char **argv) MPI_Finalize(); return 1; } + ReadUtilIface *iface; + ErrorCode rval = mb->query_interface(iface); + CHECKE("Can't get reader interface"); - MPI_Comm comm ; int rank, size; MPI_Comm_rank( MPI_COMM_WORLD, &rank ); - MPI_Comm_rank( MPI_COMM_WORLD, &size ); + MPI_Comm_size( MPI_COMM_WORLD, &size ); if (M*N*K != size) { if (rank==0) { - cout <<"M*N*K=" << M*N*K << " != size = "<< " size\n"; + cout <<"M*N*K=" << M*N*K << " != size = "<< size << "\n"; } MPI_Finalize(); return 1; } - // the global id of each vertex will be, on task + // determine m, n, k for processor rank + int m,n,k; + k = rank/(M*N); + int leftover = rank%(M*N); + n = leftover/M; + m = leftover%M; + if (rank==size-1) + { + cout << "m, n, k for last rank:" << m << " " << n << " " << k << "\n"; + } + + // so there are a total of M * A * blockSize elements in x direction (so M * A * blockSize + 1 verts in x direction) + // so there are a total of N * B * blockSize elements in y direction (so N * B * blockSize + 1 verts in y direction) + // so there are a total of K * C * blockSize elements in z direction (so K * C * blockSize + 1 verts in z direction) + + // there are ( M * A blockSize ) * ( N * B * blockSize) * (K * C * blockSize ) hexas + // there are ( M * A * blockSize + 1) * ( N * B * blockSize + 1 ) * (K * C * blockSize + 1) vertices + // x is the first dimension that varies + + int NX = ( M * A * blockSize + 1); + int NY = ( N * B * blockSize + 1); + int NZ = ( K * C * blockSize + 1); + int blockSize1 = blockSize + 1;// used for vertices + int bsq = blockSize * blockSize; + int b1sq = blockSize1 * blockSize1; + // generate the block at (a, b, c); it will represent a partition , it will get a partition tag + + Tag global_id_tag; + mb->tag_get_handle("GLOBAL_ID", 1, MB_TYPE_INTEGER, + global_id_tag); + + Tag part_tag; + int dum_id=-1; + mb->tag_get_handle("PARALLEL_PARTITION", 1, MB_TYPE_INTEGER, + part_tag, MB_TAG_CREAT|MB_TAG_SPARSE, &dum_id); + + for (int a=0; a<A; a++) + { + for (int b=0; b<B; b++) + { + for (int c=0; c<C; c++) + { + // we will generate (block+1)^3 vertices, and block^3 hexas + // the global id of the vertices will come from m, n, k, a, b, c + // x will vary from m*A*block + a*block to m*A*block+(a+1)*block etc; + int num_nodes = (blockSize+1)*(blockSize+1)*(blockSize+1); + + vector<double*> arrays; + EntityHandle startv; + rval = iface->get_node_coords(3, num_nodes, 0, startv, arrays); + CHECKE("Can't get node coords."); + + int x = m*A*blockSize + a*blockSize; + int y = n*B*blockSize + b*blockSize; + int z = k*C*blockSize + c*blockSize; + int ix=0; + vector<int> gids(num_nodes); + Range verts(startv, startv+num_nodes-1); + for (int kk=0; kk<blockSize+1; kk++) + { + for (int jj=0; jj<blockSize+1; jj++) + { + for (int ii=0; ii<blockSize+1; ii++) + { + arrays[0][ix] = (double)x+ii; + arrays[1][ix] = (double)y+jj; + arrays[2][ix] = (double)z+kk; + gids[ix] = 1 + (z+kk) * (NX*NY) + (y+jj) * NX + (x+ii); + ix++; + } + } + } + mb->tag_set_data(global_id_tag, verts, &gids[0]); + int num_hexas = (blockSize)*(blockSize)*(blockSize); + EntityHandle starte; // connectivity + EntityHandle * conn; + rval = iface->get_element_connect(num_hexas, 8, MBHEX, 0, starte, conn); + CHECKE("Can't get hexa connectivity."); + + Range hexas(starte, starte+num_hexas-1); + // fill hexas + ix=0; + for (int kk=0; kk<blockSize; kk++) + { + for (int jj=0; jj<blockSize; jj++) + { + for (int ii=0; ii<blockSize; ii++) + { + EntityHandle corner=startv + ii + (jj*blockSize1) + kk * b1sq; + conn[ix]=corner; + conn[ix+1]=corner+1; + conn[ix+2]=corner+ 1 + blockSize1; + conn[ix+3]=corner+ blockSize1; + conn[ix+4]=corner + b1sq; + conn[ix+5]=corner+1 + b1sq; + conn[ix+6]=corner+ 1 + blockSize1 + b1sq; + conn[ix+7]=corner+ blockSize1 + b1sq; + ix+=8; + } + } + } + EntityHandle part_set; + rval = mb->create_meshset(MESHSET_SET, part_set); + CHECKE("Can't create mesh set."); + rval = mb->add_entities(part_set, hexas); + CHECKE("Can't add entities to set."); + int part_num= a + m*A + (b + n*B)*(M*A) + (c+k*C)*(M*A * N*B); + rval = mb->tag_set_data(part_tag, &part_set, 1, &part_num); + CHECKE("Can't set part tag on set"); + + + + } + } + } + // after the mesh is generated on each proc, merge the vertices + MergeMesh mm(mb); + Range allhexas; + rval = mb->get_entities_by_type(0, MBHEX, allhexas); + + CHECKE("Can't get all hexa elements."); + + rval = mm.merge_entities( allhexas, 0.0001); + CHECKE("Can't merge"); + ParallelComm* pcomm = ParallelComm::get_pcomm(mb, 0); + if (pcomm==NULL) + { + pcomm = new ParallelComm( mb, MPI_COMM_WORLD ); + } + rval = pcomm->resolve_shared_ents( 0, allhexas, 3, 0 ); + CHECKE("Can't resolve shared ents"); + + // delete all quads and edges + Range toDelete; + rval = mb->get_entities_by_dimension(0, 1, toDelete); + CHECKE("Can't get edges"); + rval = mb->delete_entities(toDelete); + CHECKE("Can't delete edges"); + + toDelete.clear(); + + rval = mb->get_entities_by_dimension(0, 2, toDelete); + CHECKE("Can't get faces"); + rval = mb->delete_entities(toDelete); + CHECKE("Can't delete faces"); + + rval = pcomm->correct_shared_entities() ; + CHECKE("Can't correct local shared") + + rval = mb->write_file("test1.h5m", 0, ";;PARALLEL=WRITE_PART"); + CHECKE("Can't write in parallel"); MPI_Finalize(); return 0; Repository URL: https://bitbucket.org/fathomteam/moab/ -- This is a commit notification from bitbucket.org. You are receiving this because you have the service enabled, addressing the recipient of this email.

4 5

commit/MOAB: 3 new changesets
by commits-noreply＠bitbucket.org 30 Jun '14

30 Jun '14

3 new commits in MOAB: https://bitbucket.org/fathomteam/moab/commits/dbf0a5098d31/ Changeset: dbf0a5098d31 Branch: None User: nray Date: 2014-07-01 05:28:15 Summary: Added the HalfFacetRep class and one test code. Affected #: 7 files diff --git a/MeshFiles/unittest/Makefile.am b/MeshFiles/unittest/Makefile.am index bc4f522..2bd2a8a 100644 --- a/MeshFiles/unittest/Makefile.am +++ b/MeshFiles/unittest/Makefile.am @@ -34,5 +34,6 @@ EXTRA_DIST = 125hex.g \ BedCrop2.h5m \ mpas_p8.h5m \ Homme_2pt.h5m \ - surfrandomtris-4part.h5m + surfrandomtris-4part.h5m \ + hexes_mixed.vtk diff --git a/MeshFiles/unittest/hexes_mixed.vtk b/MeshFiles/unittest/hexes_mixed.vtk new file mode 100755 index 0000000..8ee5c11 --- /dev/null +++ b/MeshFiles/unittest/hexes_mixed.vtk @@ -0,0 +1,45 @@ +# vtk DataFile Version 3.0 +This file was written using writevtk_unstr.m +ASCII +DATASET UNSTRUCTURED_GRID +POINTS 18 double +0 -1 0 +0 1 0 +1 -1 0 +1 1 0 +-1 1 0 +-1 -1 0 +0 -1 1 +0 1 1 +1 -1 1 +1 1 1 +-1 1 1 +-1 -1 1 +0 -1 -1 +0 1 -1 +1 -1 -1 +1 1 -1 +-1 1 -1 +-1 -1 -1 + +CELLS 9 59 +2 0 1 +4 0 1 7 6 +4 0 2 3 1 +4 0 12 13 1 +4 0 1 4 5 +8 0 2 3 1 6 8 9 7 +8 5 0 1 4 11 6 7 10 +8 12 14 15 13 0 2 3 1 +8 17 12 13 16 5 0 1 4 + +CELL_TYPES 9 +3 +9 +9 +9 +9 +12 +12 +12 +12 This diff is so big that we needed to truncate the remainder. https://bitbucket.org/fathomteam/moab/commits/8fa155080b49/ Changeset: 8fa155080b49 Branch: None User: nray Date: 2014-07-01 05:28:15 Summary: Added interface to AHF functionalities through the get_adjacencies functions. Affected #: 10 files diff --git a/configure.ac b/configure.ac index 5de30e3..87ad9ee 100644 --- a/configure.ac +++ b/configure.ac @@ -657,6 +657,13 @@ MB_OPTIONAL_TOOL([mbdepth], [yes]) MB_OPTIONAL_TOOL([refiner], [no]) MB_OPTIONAL_TOOL([h5mtools], [yes]) MB_OPTIONAL_TOOL([mbcslam], [no]) +MB_OPTIONAL_TOOL([ahf], [no]) + + +if test "xyes" = "x$ENABLE_ahf"; then + AM_CPPFLAGS="$AM_CPPFLAGS -DUSE_AHF" +fi + if test "xyes" = "x$ENABLE_refiner"; then if test "xyes" != "x$WITH_MPI"; then @@ -1333,3 +1340,27 @@ if test "x$WARN_PARALLEL_HDF5_NO_COMPLEX" = "xyes"; then *************************************************************************]) fi +################################################################################ +# AHF OPTIONS +################################################################################ +#AC_ARG_ENABLE([ahf], +#[AC_HELP_STRING([--enable-ahf],[use the half facet representation])], +#[ case "${enableval}" in +# yes) ahf1=true ;; +# no) ahf1=false ;; +# *) AC_MSG_ERROR([bad value]) ;; +#esac],[ahf1=false]) + #ENABLE_AHF=$enableval],[ENABLE_AHF=no] ) + +#if test "x$ENABLE_AHF" = "xyes"; then +# AM_CPPFLAGS="$AM_CPPFLAGS -DUSE_AHF" +#fi +#AM_CONDITIONAL([ENABLE_AHF],[test x$ahf1 = xtrue]) +#AC_SUBST(ENABLE_AHF) + + +#test "xyes" != "x$WITH_AHF" || AM_CPPFLAGS="$AM_CPPFLAGS -DUSE_AHF" +#AM_CONDITIONAL(USE_AHF, [test "xno" != "x$WITH_AHF"]) + + + diff --git a/src/Core.cpp b/src/Core.cpp index 0bd882f..3641fd5 100644 --- a/src/Core.cpp +++ b/src/Core.cpp @@ -56,6 +56,10 @@ #include <errno.h> #include <string.h> +#ifdef USE_AHF +#include "moab/HalfFacetRep.hpp" +#endif + #ifdef USE_MPI /* Leave ParallelComm.hpp before mpi.h or MPICH2 will fail * because its C++ headers do not like SEEK_* macros. @@ -272,6 +276,12 @@ ErrorCode Core::initialize() geom_dimension_tag(); globalId_tag(); +#ifdef USE_AHF + ahfRep = new HalfFacetRep(this); + if (!ahfRep) + return MB_MEMORY_ALLOCATION_FAILED; +#endif + return MB_SUCCESS; } @@ -322,6 +332,11 @@ void Core::deinitialize() if (mpiFinalize) MPI_Finalize(); #endif + +#ifdef USE_AHF + delete ahfRep; + ahfRep = 0; +#endif } ErrorCode Core::query_interface_type( const std::type_info& type, void*& ptr ) @@ -1338,6 +1353,80 @@ ErrorCode get_adjacencies_union( Core* gMB, return result; } + +/////////////////////////// +/// +#ifdef USE_AHF +template <typename ITER> static inline +ErrorCode get_adjacencies_intersection_ahf(Core *mb, + ITER begin, ITER end, + const int to_dimension, + std::vector<EntityHandle>& adj_entities ) +{ + const size_t SORT_THRESHOLD = 200; + std::vector<EntityHandle> temp_vec; + std::vector<EntityHandle>::iterator adj_it, w_it; + ErrorCode result = MB_SUCCESS; + + if (begin == end) { + adj_entities.clear(); // intersection + return MB_SUCCESS; + } + + // First iteration is a special case if input list is empty. + // Rather than returning nothing (intersecting with empty + // input list), we begin with the adjacencies for the first entity. + if (adj_entities.empty()) { + EntityType type = TYPE_FROM_HANDLE(*begin); + + if(to_dimension == 0 && type != MBPOLYHEDRON) + result = mb->get_connectivity(&(*begin), 1, adj_entities); + else + result = mb->a_half_facet_rep()->get_adjacencies(*begin, to_dimension, adj_entities); + if (MB_SUCCESS != result) + return result; + ++begin; + } + + for (ITER from_it = begin; from_it != end; from_it++) + { + // running results kept in adj_entities; clear temp_vec, which is working space + temp_vec.clear(); + + // get the next set of adjacencies + EntityType type = TYPE_FROM_HANDLE(*from_it); + if(to_dimension == 0 && type != MBPOLYHEDRON) + result = mb->get_connectivity(&(*from_it), 1, temp_vec); + else + result = mb->a_half_facet_rep()->get_adjacencies(*from_it, to_dimension, temp_vec); + if (MB_SUCCESS != result) + return result; + + // otherwise intersect with the current set of results + w_it = adj_it = adj_entities.begin(); + if (temp_vec.size()*adj_entities.size() < SORT_THRESHOLD) { + for (; adj_it != adj_entities.end(); ++adj_it) + if (std::find(temp_vec.begin(), temp_vec.end(), *adj_it) != temp_vec.end()) + { *w_it = *adj_it; ++w_it; } + } + else { + std::sort( temp_vec.begin(), temp_vec.end() ); + for (; adj_it != adj_entities.end(); ++adj_it) + if (std::binary_search(temp_vec.begin(), temp_vec.end(), *adj_it)) + { *w_it = *adj_it; ++w_it; } + } + adj_entities.erase( w_it, adj_entities.end() ); + + // we're intersecting, so if there are no more results, we're done + if (adj_entities.empty()) + break; + } + + return MB_SUCCESS; +} +#endif +/////////////////////////////////////////// + template <typename ITER> static inline ErrorCode get_adjacencies_intersection( Core* mb, ITER begin, ITER end, @@ -1445,37 +1534,116 @@ ErrorCode Core::get_adjacencies( const EntityHandle *from_entities, const int to_dimension, const bool create_if_missing, std::vector<EntityHandle> &adj_entities, - const int operation_type ) + const int operation_type +#ifdef USE_AHF + , const bool use_ahf) +#else + ) +#endif { - if (operation_type == Interface::INTERSECT) - return get_adjacencies_intersection( this, from_entities, from_entities+num_entities, - to_dimension, create_if_missing, adj_entities ); - else if (operation_type != Interface::UNION) - return MB_FAILURE; - // do union - ErrorCode result; - std::vector<EntityHandle> tmp_storage; - const EntityHandle* conn; - int len; - for (int i = 0; i < num_entities; ++i) { - if(to_dimension == 0 && TYPE_FROM_HANDLE(from_entities[0]) != MBPOLYHEDRON) { - result = get_connectivity(from_entities[i], conn, len, false, &tmp_storage); - adj_entities.insert( adj_entities.end(), conn, conn+len ); - if (MB_SUCCESS != result) - return result; +#ifdef USE_AHF + bool can_handle = true; + int source_dim = this->dimension_from_handle(from_entities[0]); + + if ((source_dim > to_dimension) && (to_dimension != 0)) + { + std::cout<<"Currently Not Supported by MOAB_AHF: Down Adjacencies"<<std::endl; + std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; + can_handle = false; } - else { - result = aEntityFactory->get_adjacencies(from_entities[i], to_dimension, - create_if_missing, adj_entities); - if (MB_SUCCESS != result) - return result; + else if (((source_dim == 0) && (to_dimension == 2))||((source_dim == 0) && (to_dimension == 3))) + { + std::cout<<"Currently Not Supported by MOAB_AHF: Vertex to face/cell"<<std::endl; + std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; + can_handle = false; + } + else if (TYPE_FROM_HANDLE(from_entities[0]) == MBPOLYHEDRON) + { + std::cout<<"Currently Not Supported by MOAB_AHF: Polyhedron Meshes"<<std::endl; + std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; + can_handle = false; + } + else if (create_if_missing) + { + std::cout<<"Currently Not Supporteded by MOAB_AHF: create_if_missing capability "<<std::endl; + std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; + can_handle = false; } - } - std::sort( adj_entities.begin(), adj_entities.end() ); - adj_entities.erase( std::unique( adj_entities.begin(), adj_entities.end() ), adj_entities.end() ); - return MB_SUCCESS; + if ((use_ahf) && (can_handle)) + { + + if (operation_type == Interface::INTERSECT) + return get_adjacencies_intersection_ahf(this, from_entities, from_entities+num_entities, + to_dimension, adj_entities ); + else if (operation_type != Interface::UNION) + return MB_FAILURE; + + // do union + ErrorCode result; + std::vector<EntityHandle> tmp_storage; + const EntityHandle* conn; + int len; + for (int i = 0; i < num_entities; ++i) { + if(to_dimension == 0 && TYPE_FROM_HANDLE(from_entities[0]) != MBPOLYHEDRON) { + result = get_connectivity(from_entities[i], conn, len, false, &tmp_storage); + adj_entities.insert( adj_entities.end(), conn, conn+len ); + if (MB_SUCCESS != result) + return result; + } + else { + result = ahfRep->get_adjacencies(from_entities[i], to_dimension, adj_entities); + if (MB_SUCCESS != result) + return result; + } + } + std::sort( adj_entities.begin(), adj_entities.end() ); + adj_entities.erase( std::unique( adj_entities.begin(), adj_entities.end() ), adj_entities.end() ); + + } + else + { + +#endif + + if (operation_type == Interface::INTERSECT) + return get_adjacencies_intersection( this, from_entities, from_entities+num_entities, + to_dimension, create_if_missing, adj_entities ); + else if (operation_type != Interface::UNION) + return MB_FAILURE; + + // do union + ErrorCode result; + std::vector<EntityHandle> tmp_storage; + const EntityHandle* conn; + int len; + for (int i = 0; i < num_entities; ++i) { + if(to_dimension == 0 && TYPE_FROM_HANDLE(from_entities[0]) != MBPOLYHEDRON) { + result = get_connectivity(from_entities[i], conn, len, false, &tmp_storage); + adj_entities.insert( adj_entities.end(), conn, conn+len ); + if (MB_SUCCESS != result) + return result; + } + else { + result = aEntityFactory->get_adjacencies(from_entities[i], to_dimension, + create_if_missing, adj_entities); + if (MB_SUCCESS != result) + return result; + } + } + std::sort( adj_entities.begin(), adj_entities.end() ); + adj_entities.erase( std::unique( adj_entities.begin(), adj_entities.end() ), adj_entities.end() ); + + //return MB_SUCCESS; + + +#ifdef USE_AHF + } +#endif + +return MB_SUCCESS; + } diff --git a/src/HalfFacetRep.cpp b/src/HalfFacetRep.cpp index 95efb7e..f3818a6 100755 --- a/src/HalfFacetRep.cpp +++ b/src/HalfFacetRep.cpp @@ -17,21 +17,32 @@ #pragma warning (disable : 4786) #endif +#include "moab/HalfFacetRep.hpp" #include <iostream> #include <assert.h> -#include <time.h> #include <vector> -#include <queue> -#include <stack> #include "moab/Core.hpp" #include "moab/Range.hpp" #include "moab/CN.hpp" -#include "moab/HalfFacetRep.hpp" namespace moab { const int MAXSIZE = 500; + HalfFacetRep::HalfFacetRep(Core *impl) + { + assert(NULL != impl); + mb = impl; + mInitAHFmaps = false; + } + + HalfFacetRep::~HalfFacetRep() + { + ErrorCode result; + result = deinitialize(); + //if (MB_SUCCESS != result) return result; + } + MESHTYPE HalfFacetRep::get_mesh_type(int nverts, int nedges, int nfaces, int ncells) { MESHTYPE mesh_type; @@ -58,7 +69,9 @@ namespace moab { * initialize * ******************************************************/ - ErrorCode HalfFacetRep::initialize(){ + ErrorCode HalfFacetRep::initialize() + { + mInitAHFmaps = true; ErrorCode error; @@ -436,36 +449,100 @@ namespace moab { * User interface for adjacency functions * ********************************************************/ - ErrorCode HalfFacetRep::get_upward_incidences(EntityHandle ent, int out_dim, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * lids) + ErrorCode HalfFacetRep::get_adjacencies(const EntityHandle source_entity, + const unsigned int target_dimension, + std::vector<EntityHandle> &target_entities) + { + + ErrorCode error; + + unsigned int source_dimension = mb->dimension_from_handle(source_entity); + + if (mInitAHFmaps == false) + { + error = initialize(); + if (MB_SUCCESS != error) return error; + } + + if ((source_dimension == 0) && (target_dimension == 1)) + { + error = get_up_adjacencies_1d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + + else if ((source_dimension == 1) && (target_dimension == 2)) + { + error = get_up_adjacencies_2d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + else if ((source_dimension == 1) && (target_dimension == 3)) + { + error = get_up_adjacencies_edg_3d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + else if ((source_dimension == 2) && (target_dimension ==3)) + { + error = get_up_adjacencies_face_3d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + else if (source_dimension == target_dimension) + { + if (target_dimension == 1) + { + error = get_neighbor_adjacencies_1d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + + else if (target_dimension == 2) + { + error = get_neighbor_adjacencies_2d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + else if (target_dimension == 3) + { + error = get_neighbor_adjacencies_3d(source_entity, target_entities); + if (MB_SUCCESS != error) return error; + } + } + return MB_SUCCESS; + } + + + ErrorCode HalfFacetRep::get_up_adjacencies(EntityHandle ent, + int out_dim, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * lids) { ErrorCode error; int in_dim = mb->dimension_from_handle(ent); if ((in_dim == 0) && (out_dim == 1)) { - error = get_upward_incidences_1d(ent, adjents, local_id, lids); + error = get_up_adjacencies_1d(ent, adjents, local_id, lids); if (MB_SUCCESS != error) return error; } else if ((in_dim == 1) && (out_dim == 2)) { - error = get_upward_incidences_2d(ent, adjents, local_id, lids); + error = get_up_adjacencies_2d(ent, adjents, local_id, lids); if (MB_SUCCESS != error) return error; } else if ((in_dim == 1) && (out_dim == 3)) { - error = get_upward_incidences_edg_3d(ent, adjents, local_id, lids); + error = get_up_adjacencies_edg_3d(ent, adjents, local_id, lids); if (MB_SUCCESS != error) return error; } else if ((in_dim == 2) && (out_dim ==3)) { - error = get_upward_incidences_face_3d(ent, adjents, local_id, lids); + error = get_up_adjacencies_face_3d(ent, adjents, local_id, lids); if (MB_SUCCESS != error) return error; } return MB_SUCCESS; } - ErrorCode HalfFacetRep::get_neighbor_adjacencies(EntityHandle ent, std::vector<EntityHandle> &adjents) + ErrorCode HalfFacetRep::get_neighbor_adjacencies(EntityHandle ent, + std::vector<EntityHandle> &adjents) { ErrorCode error; int in_dim = mb->dimension_from_handle(ent); @@ -619,8 +696,11 @@ namespace moab { return MB_SUCCESS; } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_1d( EntityHandle vid, std::vector< EntityHandle > &adjents, bool local_id, std::vector<int> * lvids){ - ErrorCode error; + ErrorCode HalfFacetRep::get_up_adjacencies_1d( EntityHandle vid, + std::vector< EntityHandle > &adjents, + bool local_id, + std::vector<int> * lvids){ + ErrorCode error; EntityHandle start_eid, eid, sibeid[2]; int start_lid, lid, siblid[2]; @@ -655,7 +735,8 @@ namespace moab { return MB_SUCCESS; } ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_neighbor_adjacencies_1d( EntityHandle eid, std::vector<EntityHandle> &adjents){ + ErrorCode HalfFacetRep::get_neighbor_adjacencies_1d( EntityHandle eid, + std::vector<EntityHandle> &adjents){ ErrorCode error; @@ -912,9 +993,65 @@ namespace moab { return MB_SUCCESS; } + /////////////////////////////////////////////////////////////////// + /* ErrorCode HalfFacetRep::get_up_adjacencies_2d(EntityHandle vid, std::vector<EntityHandle> &adjents) + { + ErrorCode error; + + EntityHandle fid; int lid; + error = mb->tag_get_data(v2he_fid, &vid, 1, &fid); + if (MB_SUCCESS != error) return error; + error = mb->tag_get_data(v2he_leid, &vid, 1, &lid); + if (MB_SUCCESS != error) return error; + + if (fid != 0) + { + adjents.push_back(fid); + + EntityHandle queue_fid[MAXSIZE], trackfaces[MAXSIZE]; + int queue_lid[MAXSIZE]; + for (int i = 0; i< MAXSIZE; i++) + { + queue_fid[i] = 0; + queue_lid[i] = 0; + trackfaces[i] = 0; + } + int qsize = 0, tcount = -1; + int num_qvals = 0; + error = gather_halfedges(vid, fid, lid, queue_fid, queue_lid, &qsize, trackfaces, &tcount); + if (MB_SUCCESS != error) return error; + + while (num_qvals < *qsize) + { + EntityHandle curfid = queue_fid[num_qvals]; + int curlid = queue_lid[num_qvals]; + num_qvals += 1; + + EntityHandle he2_fid; int he2_lid; + error = another_halfedge(vid, curfid, curlid, &he2_fid, &he2_lid); + if (MB_SUCCESS != error) return error; + bool found_ent = find_match_in_array(he2_fid, trackfaces, tcount[0]); + + if (found_ent) + continue; + tcount[0] += 1; + trackfaces[tcount[0]] = he2_fid; + + error = get_up_adjacencies_2d(he2_fid, he2_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); + if (MB_SUCCESS != error) return error; + + adjents.push_back(he2_fid); + + } + } + } +*/ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_2d( EntityHandle eid, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * leids) + ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle eid, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * leids) { // Given an explicit edge eid, find the incident faces. @@ -928,7 +1065,7 @@ namespace moab { // Step 2: If there is a corresponding half-edge, collect all sibling half-edges and store the incident faces. if (found) { - error = get_upward_incidences_2d(he_fid, he_lid, true, adjents, local_id, leids); + error = get_up_adjacencies_2d(he_fid, he_lid, true, adjents, local_id, leids); if (MB_SUCCESS != error) return error; } @@ -936,7 +1073,12 @@ namespace moab { } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_2d( EntityHandle fid, int leid, bool add_inent, std::vector<EntityHandle> &fids, bool local_id, std::vector<int> * leids) + ErrorCode HalfFacetRep::get_up_adjacencies_2d( EntityHandle fid, + int leid, + bool add_inent, + std::vector<EntityHandle> &fids, + bool local_id, + std::vector<int> * leids) { // Given an implicit half-edge <fid, leid>, find the incident half-edges. ErrorCode error; @@ -982,7 +1124,13 @@ namespace moab { } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_2d(EntityHandle fid, int lid, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount) + ErrorCode HalfFacetRep::get_up_adjacencies_2d(EntityHandle fid, + int lid, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount) { ErrorCode error; @@ -1035,7 +1183,9 @@ namespace moab { } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - bool HalfFacetRep::find_matching_halfedge( EntityHandle eid, EntityHandle *hefid, int *helid){ + bool HalfFacetRep::find_matching_halfedge( EntityHandle eid, + EntityHandle *hefid, + int *helid){ ErrorCode error; std::vector<EntityHandle> conn(2); error = mb->get_connectivity(&eid, 1, conn); @@ -1071,7 +1221,14 @@ namespace moab { return found; } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::gather_halfedges( EntityHandle vid, EntityHandle he_fid, int he_lid, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount) + ErrorCode HalfFacetRep::gather_halfedges( EntityHandle vid, + EntityHandle he_fid, + int he_lid, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount) { ErrorCode error; EntityHandle he2_fid = 0; int he2_lid = 0; @@ -1086,16 +1243,20 @@ namespace moab { tcount[0] += 1; trackfaces[tcount[0]] = he_fid; - error = get_upward_incidences_2d(he_fid, he_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); + error = get_up_adjacencies_2d(he_fid, he_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); if (MB_SUCCESS != error) return error; - error = get_upward_incidences_2d(he2_fid, he2_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); + error = get_up_adjacencies_2d(he2_fid, he2_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); if (MB_SUCCESS != error) return error; return MB_SUCCESS; } //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::another_halfedge( EntityHandle vid, EntityHandle he_fid, int he_lid, EntityHandle *he2_fid, int *he2_lid) + ErrorCode HalfFacetRep::another_halfedge( EntityHandle vid, + EntityHandle he_fid, + int he_lid, + EntityHandle *he2_fid, + int *he2_lid) { ErrorCode error; int nepf; @@ -1121,7 +1282,14 @@ namespace moab { } ///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - bool HalfFacetRep::collect_and_compare(std::vector<EntityHandle> &edg_vert, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount, EntityHandle *he_fid, int *he_lid) + bool HalfFacetRep::collect_and_compare(std::vector<EntityHandle> &edg_vert, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount, + EntityHandle *he_fid, + int *he_lid) { ErrorCode error; int nepf = local_maps_2d(*_faces.begin()); @@ -1162,7 +1330,7 @@ namespace moab { error = another_halfedge(edg_vert[0], curfid, curlid, &he2_fid, &he2_lid); if (MB_SUCCESS != error) return error; - error = get_upward_incidences_2d(he2_fid, he2_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); + error = get_up_adjacencies_2d(he2_fid, he2_lid, queue_fid, queue_lid, qsize, trackfaces, tcount); if (MB_SUCCESS != error) return error; counter += 1; @@ -1174,14 +1342,15 @@ namespace moab { } ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_neighbor_adjacencies_2d( EntityHandle fid, std::vector<EntityHandle> &adjents) + ErrorCode HalfFacetRep::get_neighbor_adjacencies_2d( EntityHandle fid, + std::vector<EntityHandle> &adjents) { ErrorCode error; if (fid != 0){ int nepf = local_maps_2d(fid); for (int lid = 0; lid < nepf; ++lid){ - error = get_upward_incidences_2d(fid, lid, false, adjents); + error = get_up_adjacencies_2d(fid, lid, false, adjents); if (MB_SUCCESS != error) return error; } } @@ -1212,7 +1381,7 @@ namespace moab { int id = nepf*(*f-firstF)+l; if (!trackF[id]) { - error = get_upward_incidences_2d(*f,l, false, adj_fids, true, &adj_lids); + error = get_up_adjacencies_2d(*f,l, false, adj_fids, true, &adj_lids); if (MB_SUCCESS != error) return error; total_edges -= adj_fids.size(); @@ -1530,7 +1699,10 @@ namespace moab { } ////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_edg_3d( EntityHandle eid, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * leids) + ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle eid, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * leids) { ErrorCode error; @@ -1542,7 +1714,7 @@ namespace moab { //Find all incident cells if (found) { - error =get_upward_incidences_edg_3d(cid, leid, adjents, local_id, leids); + error =get_up_adjacencies_edg_3d(cid, leid, adjents, local_id, leids); if (MB_SUCCESS != error) return error; } @@ -1550,7 +1722,11 @@ namespace moab { } ////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_edg_3d( EntityHandle cid, int leid, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * leids) + ErrorCode HalfFacetRep::get_up_adjacencies_edg_3d( EntityHandle cid, + int leid, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * leids) { ErrorCode error; @@ -1648,7 +1824,10 @@ namespace moab { } ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_face_3d( EntityHandle fid, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * lfids) + ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle fid, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * lfids) { ErrorCode error; @@ -1657,14 +1836,18 @@ namespace moab { bool found = find_matching_halfface(fid, &cid, &lid); if (found){ - error = get_upward_incidences_face_3d(cid, lid, adjents, local_id, lfids); + error = get_up_adjacencies_face_3d(cid, lid, adjents, local_id, lfids); if (MB_SUCCESS != error) return error; } return MB_SUCCESS; } /////////////////////////////////////////// - ErrorCode HalfFacetRep::get_upward_incidences_face_3d( EntityHandle cid, int lfid, std::vector<EntityHandle> &adjents, bool local_id, std::vector<int> * lfids) + ErrorCode HalfFacetRep::get_up_adjacencies_face_3d( EntityHandle cid, + int lfid, + std::vector<EntityHandle> &adjents, + bool local_id, + std::vector<int> * lfids) { ErrorCode error; @@ -1701,7 +1884,9 @@ namespace moab { return MB_SUCCESS; } ///////////////////////////////////////////////////////////////// - bool HalfFacetRep::find_matching_implicit_edge_in_cell( EntityHandle eid, EntityHandle *cid, int *leid) + bool HalfFacetRep::find_matching_implicit_edge_in_cell( EntityHandle eid, + EntityHandle *cid, + int *leid) { ErrorCode error; diff --git a/src/Makefile.am b/src/Makefile.am index e2bf713..6fac039 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -226,6 +226,12 @@ if HAVE_CGM AM_CPPFLAGS += @CGM_CPPFLAGS@ -DCGM @MOAB_CGM_DEFINES@ endif +# if ENABLE_ahf + # libMOAB_la_SOURCES += HalfFacetRep.cpp + # nobase_libMOAB_la_include_HEADERS += moab/HalfFacetRep.hpp +# AM_CPPFLAGS += -DUSE_AHF +# endif + MBCN_protos.h: MBCN.h $(top_srcdir)/itaps/mkprotos.sh $(AM_V_GEN)$(top_srcdir)/itaps/mkprotos.sh MBCN MOAB $< $@ MOAB_FCDefs.h diff --git a/src/moab/Core.hpp b/src/moab/Core.hpp index aa7d9da..941129a 100644 --- a/src/moab/Core.hpp +++ b/src/moab/Core.hpp @@ -35,6 +35,10 @@ class EntitySequence; class FileOptions; class SetIterator; +#ifdef USE_AHF +class HalfFacetRep; +#endif + #ifdef XPCOM_MB #define MBCORE_CID \ @@ -311,12 +315,18 @@ public: get_adjacencies( from_entities, MB_1D_ENTITY, adjacencies ); \endcode */ - virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, - const int num_entities, - const int to_dimension, - const bool create_if_missing, - std::vector<EntityHandle>& adj_entities, - const int operation_type = Interface::INTERSECT); + virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, + const int num_entities, + const int to_dimension, + const bool create_if_missing, + std::vector<EntityHandle>& adj_entities, + const int operation_type = Interface::INTERSECT + #ifdef USE_AHF + , const bool use_ahf = false); + #else + ); + #endif + virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, const int num_entities, @@ -1126,7 +1136,12 @@ public: //! return the a_entity_factory pointer AEntityFactory *a_entity_factory() { return aEntityFactory; } const AEntityFactory *a_entity_factory() const { return aEntityFactory; } - + +#ifdef USE_AHF + HalfFacetRep *a_half_facet_rep() { return ahfRep; } + const HalfFacetRep *a_half_facet_rep() const {return ahfRep; } +#endif + //! return set of registered IO tools ReaderWriterSet* reader_writer_set() { return readerWriterSet; } @@ -1347,6 +1362,10 @@ private: //! list of iterators std::vector<SetIterator*> setIterators; + +#ifdef USE_AHF + HalfFacetRep *ahfRep; +#endif }; diff --git a/src/moab/HalfFacetRep.hpp b/src/moab/HalfFacetRep.hpp index eca5e91..3d1069a 100755 --- a/src/moab/HalfFacetRep.hpp +++ b/src/moab/HalfFacetRep.hpp @@ -52,14 +52,15 @@ namespace moab { VOLUME_MIXED //Volume mesh with embedded curves and surfaces }; - + class Core; + class HalfFacetRep{ public: - HalfFacetRep(Interface *impl) : mb(impl) {} + HalfFacetRep(Core *impl); - ~HalfFacetRep() {} + ~HalfFacetRep(); // User interface functions @@ -74,6 +75,12 @@ namespace moab { //! Prints the tag values. ErrorCode print_tags(); + + ErrorCode get_adjacencies(const EntityHandle source_entity, + const unsigned int target_dimension, + std::vector<EntityHandle> &target_entities); + + //! Get the upward incidences associated with an entity. /** Given an entity of dimension d, gather all the incident D(>d) dimensional entities. * Parameters: @@ -84,7 +91,11 @@ namespace moab { * \param lids Vector in which the local id's are returned. */ - ErrorCode get_upward_incidences(EntityHandle ent, int out_dim, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * lids = NULL ); + ErrorCode get_up_adjacencies(EntityHandle ent, + int out_dim, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * lids = NULL ); //! Get the same-dimensional entities connected with an entity. /** Given an entity of dimension d, gather all the entities connected via d-1 dimensional entities. @@ -94,7 +105,8 @@ namespace moab { * \param adjents Vector in which the neighbor entities are returned. */ - ErrorCode get_neighbor_adjacencies(EntityHandle ent, std::vector<EntityHandle> &adjents); + ErrorCode get_neighbor_adjacencies(EntityHandle ent, + std::vector<EntityHandle> &adjents); // 1D Maps and queries @@ -127,7 +139,10 @@ namespace moab { * \param lvids Vector returning the local vertex id's */ - ErrorCode get_upward_incidences_1d(EntityHandle vid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * lvids = NULL); + ErrorCode get_up_adjacencies_1d(EntityHandle vid, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * lvids = NULL); //! Given an edge, finds vertex-connected neighbor edges /** Given an edge, it gathers all the incident edges of each vertex of the edge. @@ -136,7 +151,8 @@ namespace moab { * \param adjents Vector returning neighbor edges */ - ErrorCode get_neighbor_adjacencies_1d(EntityHandle eid, std::vector<EntityHandle> &adjents); + ErrorCode get_neighbor_adjacencies_1d(EntityHandle eid, + std::vector<EntityHandle> &adjents); // 2D Maps and queries @@ -161,6 +177,9 @@ namespace moab { ErrorCode determine_incident_halfedges(Range &faces); + /* ErrorCode get_up_adjacencies_2d(EntityHandle vid, + std::vector<EntityHandle> &adjents);*/ + //! Given an edge, finds the faces incident on it. /** Given an edge, it first finds a matching half-edge corresponding to eid, and then * collects all the incident half-edges/faces via the sibhes map. @@ -171,7 +190,10 @@ namespace moab { * \param leids Vector returning local edge ids */ - ErrorCode get_upward_incidences_2d(EntityHandle eid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * leids = NULL); + ErrorCode get_up_adjacencies_2d(EntityHandle eid, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * leids = NULL); //! Given a half-edge <fid, leid>, finds the faces incident on it. /** @@ -184,7 +206,12 @@ namespace moab { * \param leids Vector returning local edge ids */ - ErrorCode get_upward_incidences_2d(EntityHandle fid, int leid, bool add_inent, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * leids = NULL); + ErrorCode get_up_adjacencies_2d(EntityHandle fid, + int leid, + bool add_inent, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * leids = NULL); //! Given an edge, finds edge-connected neighbor face /** Given an face, it gathers all the neighbor faces of each local edge of the face. @@ -193,7 +220,8 @@ namespace moab { * \param adjents Vector returning neighbor faces */ - ErrorCode get_neighbor_adjacencies_2d(EntityHandle fid, std::vector<EntityHandle> &adjents); + ErrorCode get_neighbor_adjacencies_2d(EntityHandle fid, + std::vector<EntityHandle> &adjents); //! Given a range of faces, finds the total number of edges. @@ -229,7 +257,8 @@ namespace moab { * \param isborder: A dense tag over all vertices of size 1. Value is true for a border vertex, otherwise is false. */ - ErrorCode determine_border_vertices( Range &cells, Tag isborder); + ErrorCode determine_border_vertices( Range &cells, + Tag isborder); //! Given an edge, finds the cells incident on it. /** Given an edge, it first finds a matching local edge in a cell corresponding to eid, and then @@ -241,7 +270,10 @@ namespace moab { * \param leids Vector returning local edge ids */ - ErrorCode get_upward_incidences_edg_3d(EntityHandle eid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * leids = NULL); + ErrorCode get_up_adjacencies_edg_3d(EntityHandle eid, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * leids = NULL); //! Given a local edge <cid, leid>, finds the cells incident on it. /** Given a local edge, it gathers all the incident cells via the sibhfs map. @@ -253,7 +285,10 @@ namespace moab { * \param leids Vector returning local edge ids */ - ErrorCode get_upward_incidences_edg_3d(EntityHandle cid, int leid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * leids = NULL); + ErrorCode get_up_adjacencies_edg_3d(EntityHandle cid, + int leid, std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * leids = NULL); //! Given an face, finds the cells incident on it. /** Given an face, it first finds a matching half-face in a cell corresponding to face, and then @@ -265,7 +300,10 @@ namespace moab { * \param leids Vector returning local face ids */ - ErrorCode get_upward_incidences_face_3d(EntityHandle fid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * lfids = NULL); + ErrorCode get_up_adjacencies_face_3d(EntityHandle fid, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * lfids = NULL); //! Given a local face <cid, lfid>, finds the cells incident on it. /** Given a local face, it gathers all the incident cells via the sibhfs map. @@ -277,7 +315,11 @@ namespace moab { * \param lfids Vector returning local face ids */ - ErrorCode get_upward_incidences_face_3d(EntityHandle cid, int lfid, std::vector<EntityHandle> &adjents, bool local_id = false, std::vector<int> * lfids = NULL); + ErrorCode get_up_adjacencies_face_3d(EntityHandle cid, + int lfid, + std::vector<EntityHandle> &adjents, + bool local_id = false, + std::vector<int> * lfids = NULL); //! Given a cell, finds face-connected neighbor cells /** Given a cell, it gathers all the neighbor cells of each local face of the cell. @@ -286,12 +328,17 @@ namespace moab { * \param adjents Vector returning neighbor cells */ - ErrorCode get_neighbor_adjacencies_3d(EntityHandle cid, std::vector<EntityHandle> &adjents); + ErrorCode get_neighbor_adjacencies_3d(EntityHandle cid, + std::vector<EntityHandle> &adjents); protected: - Interface * mb; + Core * mb; + + HalfFacetRep(); + + bool mInitAHFmaps; enum { MAX_VERTICES = 8, @@ -351,14 +398,23 @@ namespace moab { * queue of half-edges, if they do not already exist in the queue. This function is used to increment the * search space for finding a matching half-edge. * Parameters: - * \param <EntityHandle he_fid, int he_lid> The query half-edge - * \param queue_fid, queue_lid, qsize Array of faces and local edge ids. qsize is the current size of the queue_fid. + * \param he_fid EntityHandle of query half-edge + * \param he_lid Local id of query half-edge + * \param queue_fid + * \param queue_lid + * \param qsize Array of faces and local edge ids. qsize is the current size of the queue_fid. * \param trackfaces Array containing faces. If fid of an incident half-edge doesn't belong to trackfaces, * the half-edge is added to the queue. * \param tcount Current size of trackfaces */ - ErrorCode get_upward_incidences_2d(EntityHandle he_fid, int he_lid, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount); + ErrorCode get_up_adjacencies_2d(EntityHandle he_fid, + int he_lid, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount); //! Given an edge, finds a matching half-edge in the surface. /** Given an edge eid, it first collects few half-edges belonging to one-ring neighborhood of @@ -369,7 +425,9 @@ namespace moab { * \param hefid, helid: Returns the matching half-edge corresponding to the query edge. */ - bool find_matching_halfedge( EntityHandle eid, EntityHandle *hefid, int *helid); + bool find_matching_halfedge( EntityHandle eid, + EntityHandle *hefid, + int *helid); //! Gather half-edges to a queue of half-edges. /** Given a vertex vid, and a half-edge <he_fid,he_lid>, add another half-edge in the same face sharing the vertex @@ -381,7 +439,14 @@ namespace moab { * \param trackfaces, tcount */ - ErrorCode gather_halfedges( EntityHandle vid, EntityHandle he_fid, int he_lid, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount); + ErrorCode gather_halfedges( EntityHandle vid, + EntityHandle he_fid, + int he_lid, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount); //! Obtains another half-edge belonging to the same face as the input half-edge /** It uses the local maps to find another half-edge that is either incident or outgoing depending @@ -392,7 +457,11 @@ namespace moab { * \param <EntityHandle he2_fid, int he2_lid>: Returns another half-edge in the same he_fid sharing vid. */ - ErrorCode another_halfedge( EntityHandle vid, EntityHandle he_fid, int he_lid, EntityHandle *he2_fid, int *he2_lid); + ErrorCode another_halfedge( EntityHandle vid, + EntityHandle he_fid, + int he_lid, + EntityHandle *he2_fid, + int *he2_lid); //! Collect and compare to find a matching half-edge with the given edge connectivity. /** Given edge connectivity, compare to an input list of half-edges to find a matching half-edge @@ -404,7 +473,14 @@ namespace moab { * \param <EntityHandle he_fid, int he_lid>: Returns matching half-edge */ - bool collect_and_compare(std::vector<EntityHandle> &edg_vert, EntityHandle *queue_fid, int *queue_lid, int *qsize, EntityHandle *trackfaces, int *tcount, EntityHandle *he_fid, int *he_lid); + bool collect_and_compare(std::vector<EntityHandle> &edg_vert, + EntityHandle *queue_fid, + int *queue_lid, + int *qsize, + EntityHandle *trackfaces, + int *tcount, + EntityHandle *he_fid, + int *he_lid); //! The local maps for 3D entities. @@ -469,7 +545,9 @@ namespace moab { * \param leid Returns the local id of the edge corresponding to the input edge w.r.t the incident cell. */ - bool find_matching_implicit_edge_in_cell( EntityHandle eid, EntityHandle *cid, int *leid); + bool find_matching_implicit_edge_in_cell( EntityHandle eid, + EntityHandle *cid, + int *leid); //! Given a face, finds a matching local face in an incident cell. /** Find a local face with the same connectivity as the input face, belonging to an incident cell. @@ -479,9 +557,15 @@ namespace moab { * \param lfid Returns the local id of the face corresponding to the input face w.r.t the incident cell. */ - bool find_matching_halfface(EntityHandle fid, EntityHandle *cid, int *leid); + bool find_matching_halfface(EntityHandle fid, + EntityHandle *cid, + int *leid); - bool find_match_in_array(EntityHandle ent, EntityHandle *ent_list, int count, bool get_index = false, int *index = NULL); + bool find_match_in_array(EntityHandle ent, + EntityHandle *ent_list, + int count, + bool get_index = false, + int *index = NULL); }; diff --git a/src/moab/Interface.hpp b/src/moab/Interface.hpp index 09f4219..bfd5fb0 100644 --- a/src/moab/Interface.hpp +++ b/src/moab/Interface.hpp @@ -624,7 +624,12 @@ public: const int to_dimension, const bool create_if_missing, std::vector<EntityHandle>& adj_entities, - const int operation_type = Interface::INTERSECT) = 0; + const int operation_type = Interface::INTERSECT + #ifdef USE_AHF + , const bool use_ahf = false) = 0; + #else + ) = 0; + #endif //! Get the adjacencies associated with a vector of entities to entities of a specfied dimension. /** Identical to vector-based get_adjacencies function, except results are returned in a diff --git a/test/Makefile.am b/test/Makefile.am index fe7edb6..34432d0 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -61,6 +61,10 @@ if HDF5_FILE lloyd_smoother_test endif +if ENABLE_ahf + TESTS += ahf_mbintf_test +endif + # merge_test \ # input files no longer exist? # test_tag_server \ # fails @@ -154,6 +158,11 @@ lloyd_smoother_test_SOURCES = lloyd_smoother_test.cpp ahf_test_SOURCES = test_ahf_moab.cpp +if ENABLE_ahf +ahf_mbintf_test_SOURCES = test_ahf_mb_interface.cpp +endif + + if PARALLEL moab_test_CPPFLAGS += -I$(top_srcdir)/src/parallel kd_tree_test_CPPFLAGS += -I$(top_srcdir)/src/parallel diff --git a/test/test_ahf_mb_interface.cpp b/test/test_ahf_mb_interface.cpp new file mode 100755 index 0000000..9bf76ef --- /dev/null +++ b/test/test_ahf_mb_interface.cpp @@ -0,0 +1,174 @@ +/*This function tests the AHF datastructures on CST meshes*/ +#include <iostream> +#include <vector> +#include <algorithm> +#include "moab/Core.hpp" +#include "moab/Range.hpp" +#include "moab/MeshTopoUtil.hpp" +#include "moab/HalfFacetRep.hpp" +#include "TestUtil.hpp" + +using namespace moab; + +#ifdef MESHDIR +static const char example[] = STRINGIFY(MESHDIR) "/hexes_mixed.vtk"; +#else +static const char example[] = "/hexes_mixed.vtk"; +#endif + +void ahf_mbintf_test() +{ + + Core moab; + Interface* mbImpl = &moab; + MeshTopoUtil mtu(mbImpl); + + ErrorCode error = mbImpl->load_file(example); + CHECK_ERR(error); + + /*Create ranges for handles of explicit elements of the mixed mesh*/ + Range verts, edges, faces, cells; + error = mbImpl->get_entities_by_dimension( 0, 0, verts); + error = mbImpl->get_entities_by_dimension( 0, 1, edges); + error = mbImpl->get_entities_by_dimension( 0, 2, faces); + error = mbImpl->get_entities_by_dimension( 0, 3, cells); + + //Perform queries + std::vector<EntityHandle> adjents; + Range mbents, ahfents; + + //1D Queries // + //IQ1: For every vertex, obtain incident edges + for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 1, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents ); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(),mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ1: For every edge, obtain neighbor edges + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 1, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + // 2D Queries + //IQ2: For every edge, obtain incident faces + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 2, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ2: For every face, obtain neighbor faces + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 2, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + // 3D Queries + // IQ 31: For every edge, obtain incident cells + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 3, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //IQ32: For every face, obtain incident cells + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 3, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ3: For every cell, obtain neighbor cells + for (Range::iterator i = cells.begin(); i != cells.end(); ++i) { + adjents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 3, false, adjents, Interface::INTERSECT, true); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + +} + +int main(int argc, char *argv[]) +{ + int result = 0; + + argv[0] = argv[argc - argc]; // Followed read_mpas_nc.cpp test for removing warnings in serial mode about unused variables. + + result += RUN_TEST(ahf_mbintf_test); + + return result; +} + diff --git a/test/test_ahf_moab.cpp b/test/test_ahf_moab.cpp index ff97fb1..93bd947 100755 --- a/test/test_ahf_moab.cpp +++ b/test/test_ahf_moab.cpp @@ -34,7 +34,7 @@ void ahf_test() error = mbImpl->get_entities_by_dimension( 0, 3, cells); // Create an ahf instance - HalfFacetRep ahf(mbImpl); + HalfFacetRep ahf(&moab); // Call the initialize function which creates the maps for each dimension ahf.initialize(); @@ -47,7 +47,7 @@ void ahf_test() //IQ1: For every vertex, obtain incident edges for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { adjents.clear(); - error = ahf.get_upward_incidences( *i, 1, adjents); + error = ahf.get_up_adjacencies( *i, 1, adjents); CHECK_ERR(error); mbents.clear(); error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents ); @@ -82,7 +82,7 @@ void ahf_test() //IQ2: For every edge, obtain incident faces for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { adjents.clear(); - error = ahf.get_upward_incidences( *i, 2, adjents); + error = ahf.get_up_adjacencies( *i, 2, adjents); CHECK_ERR(error); mbents.clear(); error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents); @@ -117,7 +117,7 @@ void ahf_test() // IQ 31: For every edge, obtain incident cells for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { adjents.clear(); - error = ahf.get_upward_incidences( *i, 3, adjents); + error = ahf.get_up_adjacencies( *i, 3, adjents); CHECK_ERR(error); mbents.clear(); error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); @@ -134,7 +134,7 @@ void ahf_test() //IQ32: For every face, obtain incident cells for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { adjents.clear(); - error = ahf.get_upward_incidences( *i, 3, adjents); + error = ahf.get_up_adjacencies( *i, 3, adjents); CHECK_ERR(error); mbents.clear(); error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); @@ -165,7 +165,7 @@ void ahf_test() CHECK(!mbents.size()); } - ahf.deinitialize(); + //ahf.deinitialize(); } https://bitbucket.org/fathomteam/moab/commits/58f615d2bcf2/ Changeset: 58f615d2bcf2 Branch: nray/feature_ahf User: nray Date: 2014-07-01 05:28:15 Summary: Removed tag_iterate from sibling half-facet maps. Changed input to test_ahf_moab to read in meshfile from command line. Added a test code for timing and memory requirements. Affected #: 9 files diff --git a/src/Core.cpp b/src/Core.cpp index 3641fd5..692495b 100644 --- a/src/Core.cpp +++ b/src/Core.cpp @@ -1540,6 +1540,12 @@ ErrorCode Core::get_adjacencies( const EntityHandle *from_entities, #else ) #endif +/*ErrorCode Core::get_adjacencies( const EntityHandle *from_entities, + const int num_entities, + const int to_dimension, + const bool create_if_missing, + std::vector<EntityHandle> &adj_entities, + const int operation_type)*/ { #ifdef USE_AHF @@ -1558,6 +1564,12 @@ ErrorCode Core::get_adjacencies( const EntityHandle *from_entities, std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; can_handle = false; } + else if (to_dimension == 4) + { + std::cout<<"Currently Not Supported by MOAB_AHF: meshsets"<<std::endl; + std::cout<<"Reverting to MOAB adjacency functionality"<<std::endl; + can_handle = false; + } else if (TYPE_FROM_HANDLE(from_entities[0]) == MBPOLYHEDRON) { std::cout<<"Currently Not Supported by MOAB_AHF: Polyhedron Meshes"<<std::endl; @@ -1571,7 +1583,9 @@ ErrorCode Core::get_adjacencies( const EntityHandle *from_entities, can_handle = false; } + if ((use_ahf) && (can_handle)) + //if (can_handle) { if (operation_type == Interface::INTERSECT) diff --git a/src/HalfFacetRep.cpp b/src/HalfFacetRep.cpp index f3818a6..449a48a 100755 --- a/src/HalfFacetRep.cpp +++ b/src/HalfFacetRep.cpp @@ -27,7 +27,7 @@ namespace moab { - const int MAXSIZE = 500; + const int MAXSIZE = 150; HalfFacetRep::HalfFacetRep(Core *impl) { @@ -573,15 +573,6 @@ namespace moab { { ErrorCode error; - EntityHandle start_edge = *edges.begin(); - int count, *sibhvs_lvid_ptr; EntityHandle *sibhvs_eid_ptr; - error = mb->tag_iterate(sibhvs_eid, edges.begin(), edges.end(), count, (void*&)sibhvs_eid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) edges.size()); - error = mb->tag_iterate(sibhvs_lvid, edges.begin(), edges.end(), count, (void*&)sibhvs_lvid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) edges.size()); - //Step 1: Create an index list storing the starting position for each vertex int nv = _verts.size(); int *is_index = new int[nv+1]; @@ -595,8 +586,10 @@ namespace moab { error = mb->get_connectivity(&*eid, 1, conn); if (MB_SUCCESS != error) return error; - is_index[conn[0]-*_verts.begin()+1] += 1; - is_index[conn[1]-*_verts.begin()+1] += 1; + int index = _verts.index(conn[0]); + is_index[index+1] += 1; + index = _verts.index(conn[1]); + is_index[index+1] += 1; } is_index[0] = 0; @@ -615,7 +608,7 @@ namespace moab { for (int j = 0; j< 2; j++) { - int v = conn[j] - *_verts.begin(); + int v = _verts.index(conn[j]); v2hv_map_eid[is_index[v]] = *eid; v2hv_map_lvid[is_index[v]] = j; is_index[v] += 1; @@ -633,26 +626,59 @@ namespace moab { error = mb->get_connectivity(&*eid, 1, conn); if (MB_SUCCESS != error) return error; + EntityHandle sibeid[2]; + int siblvid[2]; + + error = mb->tag_get_data(sibhvs_eid, &*eid, 1, sibeid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhvs_lvid, &*eid, 1, siblvid); + if (MB_SUCCESS != error) return error; + for (int k =0; k<2; k++) { - if (sibhvs_eid_ptr[2*(*eid-start_edge)+k] != 0) + + if (sibeid[k] != 0) continue; - int v = conn[k] - *_verts.begin(); + int v = _verts.index(conn[k]); int last = is_index[v+1] - 1; if (last > is_index[v]) { EntityHandle prev_eid = v2hv_map_eid[last]; int prev_lvid = v2hv_map_lvid[last]; + + EntityHandle seteid[2]; + int setlvid[2]; + error = mb->tag_get_data(sibhvs_eid, &prev_eid, 1, seteid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhvs_lvid, &prev_eid, 1, setlvid); + if (MB_SUCCESS != error) return error; + for (int i=is_index[v]; i<=last; i++) { EntityHandle cur_eid = v2hv_map_eid[i]; int cur_lvid = v2hv_map_lvid[i]; - sibhvs_eid_ptr[2*(prev_eid-start_edge) + prev_lvid] = cur_eid; - sibhvs_lvid_ptr[2*(prev_eid-start_edge) + prev_lvid] = cur_lvid; + + + seteid[prev_lvid] = cur_eid; + setlvid[prev_lvid] = cur_lvid; + + error = mb->tag_set_data(sibhvs_eid, &prev_eid, 1, seteid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_set_data(sibhvs_lvid, &prev_eid, 1, setlvid); + if (MB_SUCCESS != error) return error; prev_eid = cur_eid; prev_lvid = cur_lvid; + error = mb->tag_get_data(sibhvs_eid, &prev_eid, 1, seteid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhvs_lvid, &prev_eid, 1, setlvid); + if (MB_SUCCESS != error) return error; + } } } @@ -663,6 +689,7 @@ namespace moab { delete [] is_index; delete [] v2hv_map_eid; delete [] v2hv_map_lvid; + return MB_SUCCESS; } @@ -797,7 +824,7 @@ namespace moab { int nepf; if (type == MBTRI) nepf = 3; else if (type ==MBQUAD) nepf = 4; - + return nepf; } ///////////////////////////////////////////////////////////////////////////////// @@ -807,16 +834,13 @@ namespace moab { // next: Local ids of next edges // prev: Local ids of prev edges - if (nepf == 3) - { - next[0] = 1; next[1] = 2; next[2] = 0; - prev[0] = 2; prev[1] = 0; prev[2] = 1; - } - else if (nepf == 4) - { - next[0] = 1; next[1] = 2; next[2] = 3; next[3] = 0; - prev[0] = 3; prev[1] = 0; prev[2] = 1; prev[3] = 2; - } + for (int k=0; k<nepf-1; k++) + { + next[k]=k+1; + prev[k+1]=k; + } + next[nepf-1]=0; + prev[0] = nepf-1; return MB_SUCCESS; } @@ -826,15 +850,6 @@ namespace moab { ErrorCode error; EntityHandle start_face = *faces.begin(); - // Create pointers to the tags for direct access to memory - int count, *sibhes_leid_ptr; EntityHandle *sibhes_fid_ptr; - error = mb->tag_iterate(sibhes_fid, faces.begin(), faces.end(), count, (void*&)sibhes_fid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) faces.size()); - error = mb->tag_iterate(sibhes_leid, faces.begin(), faces.end(), count, (void*&)sibhes_leid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) faces.size()); - int nepf = local_maps_2d(start_face); int * next = new int[nepf]; int * prev = new int[nepf]; @@ -847,6 +862,7 @@ namespace moab { for (int i =0; i<nv+1; i++) is_index[i] = 0; + int index; std::vector<EntityHandle> conn(nepf); for (Range::iterator fid = faces.begin(); fid != faces.end(); ++fid) { @@ -856,7 +872,8 @@ namespace moab { for (int i = 0; i<nepf; i++) { - is_index[conn[i]-*_verts.begin()+1] += 1; + index = _verts.index(conn[i]); + is_index[index+1] += 1; } } is_index[0] = 0; @@ -877,7 +894,7 @@ namespace moab { for (int j = 0; j< nepf; j++) { - int v = conn[j] - *_verts.begin(); + int v = _verts.index(conn[j]); v2nv[is_index[v]] = conn[next[j]]; v2he_map_fid[is_index[v]] = *fid; v2he_map_leid[is_index[v]] = j; @@ -889,6 +906,7 @@ namespace moab { is_index[i+1] = is_index[i]; is_index[0] = 0; + //Step 3: Fill up sibling half-verts map for (Range::iterator fid = faces.begin(); fid != faces.end(); ++fid) { @@ -896,13 +914,22 @@ namespace moab { error = mb->get_connectivity(&*fid, 1, conn); if (MB_SUCCESS != error) return error; + EntityHandle *sibfid = new EntityHandle[nepf]; + int *sibleid = new int[nepf]; + + error = mb->tag_get_data(sibhes_fid, &*fid, 1, sibfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhes_leid, &*fid, 1, sibleid); + if (MB_SUCCESS != error) return error; + for (int k =0; k<nepf; k++) { - if (sibhes_fid_ptr[nepf*(*fid-start_face)+k] != 0) + if (sibfid[k] != 0) continue; - int v = conn[k] - *_verts.begin(); - int vn = conn[next[k]] - *_verts.begin(); + int v = _verts.index(conn[k]); + int vn = _verts.index(conn[next[k]]); EntityHandle first_fid = *fid; int first_leid = k; @@ -910,38 +937,87 @@ namespace moab { EntityHandle prev_fid = *fid; int prev_leid = k; - for (int index = is_index[vn]; index <= is_index[vn+1]-1; index++) + EntityHandle *setfid = new EntityHandle[nepf]; + int *setleid = new int[nepf]; + + for (index = is_index[vn]; index <= is_index[vn+1]-1; index++) { if (v2nv[index] == conn[k]) { EntityHandle cur_fid = v2he_map_fid[index]; int cur_leid = v2he_map_leid[index]; - sibhes_fid_ptr[nepf*(prev_fid-start_face)+prev_leid] = cur_fid; - sibhes_leid_ptr[nepf*(prev_fid-start_face)+prev_leid] = cur_leid; + + error = mb->tag_get_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + + setfid[prev_leid] = cur_fid; + setleid[prev_leid] = cur_leid; + + error = mb->tag_set_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_set_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + prev_fid = cur_fid; prev_leid = cur_leid; + } } - for (int index = is_index[v]; index <= is_index[v+1]-1; index++) + for (index = is_index[v]; index <= is_index[v+1]-1; index++) { if ((v2nv[index] == conn[next[k]])&&(v2he_map_fid[index] != *fid)) { + EntityHandle cur_fid = v2he_map_fid[index]; int cur_leid = v2he_map_leid[index]; - sibhes_fid_ptr[nepf*(prev_fid-start_face)+prev_leid] = cur_fid; - sibhes_leid_ptr[nepf*(prev_fid-start_face)+prev_leid] = cur_leid; + + error = mb->tag_get_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + + setfid[prev_leid] = cur_fid; + setleid[prev_leid] = cur_leid; + + error = mb->tag_set_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_set_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + prev_fid = cur_fid; prev_leid = cur_leid; } } if (prev_fid != first_fid){ - sibhes_fid_ptr[nepf*(prev_fid-start_face) + prev_leid] = first_fid; - sibhes_leid_ptr[nepf*(prev_fid-start_face) + prev_leid] = first_leid; - } + error = mb->tag_get_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + + setfid[prev_leid] = first_fid; + setleid[prev_leid] = first_leid; + + error = mb->tag_set_data(sibhes_fid, &prev_fid, 1, setfid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_set_data(sibhes_leid, &prev_fid, 1, setleid); + if (MB_SUCCESS != error) return error; + } + delete [] setfid; + delete [] setleid; } + delete [] sibfid; + delete [] sibleid; } delete [] next; @@ -1311,7 +1387,7 @@ namespace moab { error = mb->get_connectivity(&curfid, 1, conn); if (MB_SUCCESS != error) return error; - int id = next[curlid]; + int id = next[curlid]; if (((conn[curlid]==edg_vert[0])&&(conn[id]==edg_vert[1]))||((conn[curlid]==edg_vert[1])&&(conn[id]==edg_vert[0]))){ *he_fid = curfid; *he_lid = curlid; @@ -1336,7 +1412,7 @@ namespace moab { counter += 1; } - delete [] next; + delete [] next; delete [] prev; return found; } @@ -1417,7 +1493,7 @@ namespace moab { } ///////////////////////////////////////////// - HalfFacetRep::LocalMaps3D HalfFacetRep::lConnMap3D[4] = + const HalfFacetRep::LocalMaps3D HalfFacetRep::lConnMap3D[4] = { // Tet {4, 6, 4, {3,3,3,3}, {{0,1,3},{1,2,3},{2,0,3},{0,2,1}}, {3,3,3,3}, {{0,2,3},{0,1,3},{1,2,3},{0,1,2}}, {{0,1},{1,2},{2,0},{0,3},{1,3},{2,3}}, {{3,0},{3,1},{3,2},{0,2},{0,1},{1,2}}, {{0,4,3},{1,5,4},{2,3,5},{2,1,0}}, {{-1,0,2,3},{0,-1,1,4},{2,1,-1,5},{3,4,5,-1}}}, @@ -1439,16 +1515,6 @@ namespace moab { { ErrorCode error; EntityHandle start_cell = *cells.begin(); - - int count, *sibhfs_lfid_ptr; - EntityHandle *sibhfs_cid_ptr; - error = mb->tag_iterate(sibhfs_cid, cells.begin(), cells.end(), count, (void*&)sibhfs_cid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) cells.size()); - error = mb->tag_iterate(sibhfs_lfid, cells.begin(), cells.end(), count, (void*&)sibhfs_lfid_ptr); - if (MB_SUCCESS != error) return error; - assert(count == (int) cells.size()); - int index = get_index_from_type(start_cell); int nvpc = lConnMap3D[index].num_verts_in_cell; int nfpc = lConnMap3D[index].num_faces_in_cell; @@ -1459,6 +1525,7 @@ namespace moab { for (int i =0; i<nv+1; i++) is_index[i] = 0; + int vindex; std::vector<EntityHandle> conn(nvpc); for (Range::iterator cid = cells.begin(); cid != cells.end(); ++cid) { @@ -1476,7 +1543,8 @@ namespace moab { if (v <= conn[id]) v = conn[id]; } - is_index[v-*_verts.begin()+1] += 1; + vindex = _verts.index(v); + is_index[vindex+1] += 1; } } is_index[0] = 0; @@ -1518,7 +1586,7 @@ namespace moab { error = local_maps_2d(nvF, next, prev); if (MB_SUCCESS != error) return error; - int v = vmax - *_verts.begin(); + int v = _verts.index(vmax); v2oe_v1[is_index[v]] = vs[next[lv]]; v2oe_v2[is_index[v]] = vs[prev[lv]]; v2hf_map_cid[is_index[v]] = *cid; @@ -1542,9 +1610,18 @@ namespace moab { error = mb->get_connectivity(&*cid, 1, conn); if (MB_SUCCESS != error) return error; + EntityHandle *sibcid = new EntityHandle[nfpc]; + int *siblfid = new int[nfpc]; + + error = mb->tag_get_data(sibhfs_cid, &*cid, 1, sibcid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhfs_lfid, &*cid, 1, siblfid); + if (MB_SUCCESS != error) return error; + for (int i =0; i<nfpc; i++) { - if (sibhfs_cid_ptr[nfpc*(*cid-start_cell)+i] != 0) + if (sibcid[i] != 0) continue; @@ -1568,31 +1645,60 @@ namespace moab { error = local_maps_2d(nvF, next, prev); if (MB_SUCCESS != error) return error; - int v = vmax - *_verts.begin(); + int v = _verts.index(vmax); EntityHandle v1 = vs[prev[lv]]; EntityHandle v2 = vs[next[lv]]; + EntityHandle *setcid = new EntityHandle[nfpc]; + int *setlfid = new int[nfpc]; + for (int ind = is_index[v]; ind <= is_index[v+1]-1; ind++) { if ((v2oe_v1[ind] == v1)&&(v2oe_v2[ind] == v2)) { - + // Map to opposite hf EntityHandle cur_cid = v2hf_map_cid[ind]; int cur_lfid = v2hf_map_lfid[ind]; - sibhfs_cid_ptr[nfpc*(*cid-start_cell)+i] = cur_cid; - sibhfs_lfid_ptr[nfpc*(*cid-start_cell)+i] = cur_lfid; + error = mb->tag_get_data(sibhfs_cid, &*cid, 1, setcid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhfs_lfid, &*cid, 1, setlfid); + if (MB_SUCCESS != error) return error; + + setcid[i] = cur_cid; + setlfid[i] = cur_lfid; + + error = mb->tag_set_data(sibhfs_cid, &*cid, 1, setcid); + if (MB_SUCCESS != error) return error; + error = mb->tag_set_data(sibhfs_lfid, &*cid, 1, setlfid); + if (MB_SUCCESS != error) return error; + + //Map opposite hf to current cell + error = mb->tag_get_data(sibhfs_cid, &cur_cid, 1, setcid); + if (MB_SUCCESS != error) return error; + + error = mb->tag_get_data(sibhfs_lfid, &cur_cid, 1, setlfid); + if (MB_SUCCESS != error) return error; - sibhfs_cid_ptr[nfpc*(cur_cid-start_cell)+cur_lfid] = *cid; - sibhfs_lfid_ptr[nfpc*(cur_cid-start_cell)+cur_lfid] = i; + setcid[cur_lfid] = *cid; + setlfid[cur_lfid] = i; + error = mb->tag_set_data(sibhfs_cid, &cur_cid, 1, setcid); + if (MB_SUCCESS != error) return error; + error = mb->tag_set_data(sibhfs_lfid, &cur_cid, 1, setlfid); + if (MB_SUCCESS != error) return error; } } delete [] next; delete [] prev; delete [] vs; + delete [] setcid; + delete [] setlfid; } + delete sibcid; + delete siblfid; } @@ -2115,10 +2221,8 @@ namespace moab { bool HalfFacetRep::find_match_in_array(EntityHandle ent, EntityHandle *ent_list, int count, bool get_index, int *index) { bool found = false; - // std::cout<<"count = "<<count<<std::endl; for (int i = 0; i<= count; i++) { - // std::cout<<"For i = "<<i<<": ent-ent_list[i] = "<<(int)(ent-ent_list[i])<<std::endl; if (!((int)(ent - ent_list[i]))) { found = true; diff --git a/src/moab/Core.hpp b/src/moab/Core.hpp index 941129a..20aafdc 100644 --- a/src/moab/Core.hpp +++ b/src/moab/Core.hpp @@ -328,6 +328,16 @@ public: #endif + /*virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, + const int num_entities, + const int to_dimension, + const bool create_if_missing, + std::vector<EntityHandle>& adj_entities, + const int operation_type = Interface::INTERSECT);*/ + + + + virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, const int num_entities, const int to_dimension, diff --git a/src/moab/HalfFacetRep.hpp b/src/moab/HalfFacetRep.hpp index 3d1069a..344dbd1 100755 --- a/src/moab/HalfFacetRep.hpp +++ b/src/moab/HalfFacetRep.hpp @@ -533,7 +533,7 @@ namespace moab { int lookup_leids[MAX_VERTICES][MAX_VERTICES]; }; - static LocalMaps3D lConnMap3D[4]; + static const LocalMaps3D lConnMap3D[4]; int get_index_from_type(EntityHandle cid); diff --git a/src/moab/Interface.hpp b/src/moab/Interface.hpp index bfd5fb0..57252e9 100644 --- a/src/moab/Interface.hpp +++ b/src/moab/Interface.hpp @@ -630,6 +630,12 @@ public: #else ) = 0; #endif + /*virtual ErrorCode get_adjacencies(const EntityHandle *from_entities, + const int num_entities, + const int to_dimension, + const bool create_if_missing, + std::vector<EntityHandle>& adj_entities, + const int operation_type = Interface::INTERSECT) = 0; */ //! Get the adjacencies associated with a vector of entities to entities of a specfied dimension. /** Identical to vector-based get_adjacencies function, except results are returned in a diff --git a/test/Makefile.am b/test/Makefile.am index 34432d0..1354e3c 100644 --- a/test/Makefile.am +++ b/test/Makefile.am @@ -61,9 +61,9 @@ if HDF5_FILE lloyd_smoother_test endif -if ENABLE_ahf - TESTS += ahf_mbintf_test -endif +# if ENABLE_ahf +# TESTS += ahf_mbintf_test +# endif # merge_test \ # input files no longer exist? # test_tag_server \ # fails @@ -158,9 +158,9 @@ lloyd_smoother_test_SOURCES = lloyd_smoother_test.cpp ahf_test_SOURCES = test_ahf_moab.cpp -if ENABLE_ahf -ahf_mbintf_test_SOURCES = test_ahf_mb_interface.cpp -endif +# if ENABLE_ahf +# ahf_mbintf_test_SOURCES = test_ahf_mb_interface.cpp +# endif if PARALLEL diff --git a/test/perf/Makefile.am b/test/perf/Makefile.am index f946f94..32e2a57 100644 --- a/test/perf/Makefile.am +++ b/test/perf/Makefile.am @@ -7,13 +7,14 @@ AM_CPPFLAGS += -I$(top_srcdir)/src \ LDADD = $(top_builddir)/src/libMOAB.la -check_PROGRAMS = perf seqperf adj_time perftool +check_PROGRAMS = perf seqperf adj_time perftool ahf_mem_time noinst_PROGRAMS = perf_SOURCES = perf.cpp seqperf_SOURCES = seqperf.cpp adj_time_SOURCES = adj_time.cpp perftool_SOURCES = perftool.cpp +ahf_mem_time_SOURCES = test_ahf_mem_time.cpp if ENABLE_imesh LDADD += $(top_builddir)/itaps/imesh/libiMesh.la diff --git a/test/perf/test_ahf_mem_time.cpp b/test/perf/test_ahf_mem_time.cpp new file mode 100755 index 0000000..9c72ae7 --- /dev/null +++ b/test/perf/test_ahf_mem_time.cpp @@ -0,0 +1,281 @@ +/*This function tests the AHF datastructures on CST meshes*/ +#include <iostream> +#include <assert.h> +#include <time.h> +#include <vector> +#include "moab/Core.hpp" +#include "moab/Range.hpp" +#include "moab/MeshTopoUtil.hpp" +#include "moab/HalfFacetRep.hpp" +#include <sys/time.h> + +using namespace moab; + +#ifdef MESHDIR +std::string TestDir(STRINGIFY(MESHDIR)); +#else +std::string TestDir("."); +#endif + +std::string filename; + +double wtime() { + double y = -1; + struct timeval cur_time; + gettimeofday(&cur_time, NULL); + y = (double)(cur_time.tv_sec) + (double)(cur_time.tv_usec)*1.e-6; + return (y); +} + +int main(int argc, char **argv) +{ + // Read the input mesh + filename = TestDir + "/hexes_mixed.vtk"; + + if (argc==1) + std::cout<<"Using default input file:"<<filename<<std::endl; + else if (argc==2) + filename = argv[1]; + else { + std::cerr << "Usage: " << argv[0] << " [filename]" << std::endl; + return 1; + } + + ErrorCode error; + Core moab; + Interface* mbImpl = &moab; + MeshTopoUtil mtu(mbImpl); + + error = mbImpl->load_file( filename.c_str()); + if (MB_SUCCESS != error) { + std::cerr << filename <<": failed to load file." << std::endl; + return error; + } + + //Create ranges for handles of explicit elements of the mixed mesh + Range verts, edges, faces, cells; + error = mbImpl->get_entities_by_dimension( 0, 0, verts); + error = mbImpl->get_entities_by_dimension( 0, 1, edges); + error = mbImpl->get_entities_by_dimension( 0, 2, faces); + error = mbImpl->get_entities_by_dimension( 0, 3, cells); + + int nverts = verts.size(); + int nedges = edges.size(); + int nfaces = faces.size(); + int ncells = cells.size(); + + std::cout<<"nverts = "<<nverts<<", nedges = "<<nedges<<", nfaces = "<<nfaces<<", ncells = "<<ncells<<std::endl; + + + //Storage Costs before calling ahf functionalities + unsigned long sTotS, sTAS, sES, sAES, sAS, sAAS, sTS, sATS; + sTotS = sTAS = sES = sAES = sAS = sAAS = sTS = sATS = 0; + mbImpl->estimated_memory_use(NULL, 0, &sTotS, &sTAS, &sES, &sAES, &sAS, &sAAS, NULL, 0, &sTS, &sATS); + std::cout<<std::endl; + std::cout<<"Total storage = "<<sTotS<<std::endl; + std::cout<<"Total amortized storage = "<<sTAS<<std::endl; + std::cout<<"Entity storage = "<<sES<<std::endl; + std::cout<<"Amortized entity storage = "<<sAES<<std::endl; + std::cout<<"Adjacency storage = "<<sAS<<std::endl; + std::cout<<"Amortized adjacency storage = "<<sAAS<<std::endl; + std::cout<<"Tag storage = "<<sTS<<std::endl; + std::cout<<"Amortized tag storage = "<<sATS<<std::endl; + std::cout<<std::endl; + + + double time_start, time_elapsed, time_avg; + + // Create an ahf instance + HalfFacetRep ahf(&moab); + + // Call the initialize function which creates the maps for each dimension + time_start = wtime(); + + ahf.initialize(); + + time_elapsed = wtime() - time_start; + std::cout << "Time taken to construct the MDS = "<<time_elapsed<<" secs"<<std::endl; + + //Storage Costs after calling ahf initialize + unsigned long TotS, TAS, ES, AES, AS, AAS, TS, ATS; + TotS = TAS = ES = AES = AS = AAS = TS = ATS = 0; + mbImpl->estimated_memory_use(NULL, 0, &TotS, &TAS, &ES, &AES, &AS, &AAS, NULL, 0, &TS, &ATS); + std::cout<<std::endl; + std::cout<<"Total storage = "<<TotS<<std::endl; + std::cout<<"Total amortized storage = "<<TAS<<std::endl; + std::cout<<"Entity storage = "<<ES<<std::endl; + std::cout<<"Amortized entity storage = "<<AES<<std::endl; + std::cout<<"Adjacency storage = "<<AS<<std::endl; + std::cout<<"Amortized adjacency storage = "<<AAS<<std::endl; + std::cout<<"Tag storage = "<<TS<<std::endl; + std::cout<<"Amortized tag storage = "<<ATS<<std::endl; + std::cout<<std::endl; + + //Perform queries + std::vector<EntityHandle> adjents; + std::vector<int> lids; + Range mbents; + + //1D Queries // + //IQ1: For every vertex, obtain incident edges + std::cout<<"1D QUERIES"<<std::endl; + time_start = wtime(); + for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 1, adjents); + } + time_avg = (wtime()-time_start)/(double)verts.size(); + std::cout<<"MOAB_AHF: Average time taken to compute incident edges to a vertex = "<< time_avg<<" secs" <<std::endl; + + error = mbImpl->get_adjacencies( &*verts.begin(), 1, 1, false, mbents ); + time_start = wtime(); + for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents ); + } + time_avg = (wtime()-time_start)/(double)verts.size(); + std::cout<<"MOAB: Average time taken to compute incident edges to a vertex = "<< time_avg<<" secs" <<std::endl; + std::cout<<std::endl; + + + //NQ1: For every edge, obtain neighbor edges + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB_AHF: Average time taken to compute neighbor edges of an edge = "<<time_avg<<" secs" << std::endl; + + error = mtu.get_bridge_adjacencies( *edges.begin(), 0, 1, mbents); + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB: Average time taken to compute neighbor edges of an edge = "<<time_avg<<" secs" << std::endl; + std::cout<<std::endl; + + + // 2D Queries + //IQ21: For every edge, obtain incident faces + std::cout<<"2D QUERIES"<<std::endl; + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 2, adjents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB_AHF: Average time taken to compute incident faces on an edge = "<<time_avg<<" secs" <<std::endl; + + error = mbImpl->get_adjacencies( &*edges.begin(), 1, 2, false, mbents); + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB: Average time taken to compute incident faces on an edge = "<<time_avg<<" secs" <<std::endl; + std::cout<<std::endl; + + //NQ2: For every face, obtain neighbor faces + time_start = wtime(); + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + } + time_avg = (wtime()-time_start)/(double)faces.size(); + std::cout<<"MOAB_AHF: Average time taken to compute neighbor faces of a face = "<< time_avg<<" secs" <<std::endl; + + error = mtu.get_bridge_adjacencies( *faces.begin(), 1, 2, mbents); + time_start = wtime(); + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents); + } + time_avg = (wtime()-time_start)/(double)faces.size(); + std::cout<<"MOAB: Average time taken to compute neighbor faces of a face = "<< time_avg<<" secs" <<std::endl; + std::cout<<std::endl; + + + // 3D Queries + // IQ 31: For every edge, obtain incident cells + std::cout<<"3D QUERIES"<<std::endl; + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 3, adjents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB_AHF: Average time taken to compute incident cells on an edge = "<<time_avg <<" secs"<<std::endl; + + error = mbImpl->get_adjacencies(&*edges.begin(), 1, 3, false, mbents); + time_start = wtime(); + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + } + time_avg = (wtime()-time_start)/(double)edges.size(); + std::cout<<"MOAB: Average time taken to compute incident cells on an edge = "<<time_avg <<" secs"<<std::endl; + std::cout<<std::endl; + + //IQ32: For every face, obtain incident cells + time_start = wtime(); + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 3, adjents); + + } + time_avg = (wtime()-time_start)/(double)faces.size(); + std::cout<<"MOAB_AHF: Average time taken to compute incident cells on a face = "<<time_avg <<" secs"<<std::endl; + + error = mbImpl->get_adjacencies(&*faces.begin(), 1, 3, false, mbents); + time_start = wtime(); + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + } + time_avg = (wtime()-time_start)/(double)faces.size(); + std::cout<<"MOAB: Average time taken to compute incident cells on a face = "<<time_avg <<" secs"<<std::endl; + std::cout<<std::endl; + + //NQ3: For every cell, obtain neighbor cells + time_start = wtime(); + for (Range::iterator i = cells.begin(); i != cells.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + } + time_avg = (wtime()-time_start)/(double)cells.size(); + std::cout<<"MOAB_AHF: Average time taken to compute neighbor cells of a cell = "<< time_avg <<" secs" << std::endl; + + error = mtu.get_bridge_adjacencies( *cells.begin(), 2, 3, mbents); + time_start = wtime(); + for (Range::iterator i = cells.begin(); i != cells.end(); ++i) { + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents); + } + time_avg = (wtime()-time_start)/(double)cells.size(); + std::cout<<"MOAB: Average time taken to compute neighbor cells of a cell = "<< time_avg <<" secs" << std::endl; + std::cout<<std::endl; + + ahf.deinitialize(); + + //Storage Costs after calling ahf deinitialize + unsigned long eTotS, eTAS, eES, eAES, eAS, eAAS, eTS, eATS; + eTotS = eTAS = eES = eAES = eAS = eAAS = eTS = eATS = 0; + mbImpl->estimated_memory_use(NULL, 0, &eTotS, &eTAS, &eES, &eAES, &eAS, &eAAS, NULL, 0, &eTS, &eATS); + std::cout<<std::endl; + std::cout<<"Total storage = "<<eTotS<<std::endl; + std::cout<<"Total amortized storage = "<<eTAS<<std::endl; + std::cout<<"Entity storage = "<<eES<<std::endl; + std::cout<<"Amortized entity storage = "<<eAES<<std::endl; + std::cout<<"Adjacency storage = "<<eAS<<std::endl; + std::cout<<"Amortized adjacency storage = "<<eAAS<<std::endl; + std::cout<<"Tag storage = "<<eTS<<std::endl; + std::cout<<"Amortized tag storage = "<<eATS<<std::endl; + std::cout<<std::endl; + + return 0; +} + diff --git a/test/test_ahf_moab.cpp b/test/test_ahf_moab.cpp index 93bd947..e415186 100755 --- a/test/test_ahf_moab.cpp +++ b/test/test_ahf_moab.cpp @@ -11,172 +11,203 @@ using namespace moab; #ifdef MESHDIR -static const char example[] = STRINGIFY(MESHDIR) "/hexes_mixed.vtk"; +std::string TestDir(STRINGIFY(MESHDIR)); #else -static const char example[] = "/hexes_mixed.vtk"; +std::string TestDir("."); #endif -void ahf_test() -{ - - Core moab; - Interface* mbImpl = &moab; - MeshTopoUtil mtu(mbImpl); - - ErrorCode error = mbImpl->load_file(example); - CHECK_ERR(error); - - /*Create ranges for handles of explicit elements of the mixed mesh*/ - Range verts, edges, faces, cells; - error = mbImpl->get_entities_by_dimension( 0, 0, verts); - error = mbImpl->get_entities_by_dimension( 0, 1, edges); - error = mbImpl->get_entities_by_dimension( 0, 2, faces); - error = mbImpl->get_entities_by_dimension( 0, 3, cells); - - // Create an ahf instance - HalfFacetRep ahf(&moab); - - // Call the initialize function which creates the maps for each dimension - ahf.initialize(); - - //Perform queries - std::vector<EntityHandle> adjents; - Range mbents, ahfents; - - //1D Queries // - //IQ1: For every vertex, obtain incident edges - for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { - adjents.clear(); - error = ahf.get_up_adjacencies( *i, 1, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents ); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(),mbents.size()); - - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } - - //NQ1: For every edge, obtain neighbor edges - for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { - adjents.clear(); - error = ahf.get_neighbor_adjacencies( *i, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(), mbents.size()); - - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } - - // 2D Queries - //IQ2: For every edge, obtain incident faces - for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { - adjents.clear(); - error = ahf.get_up_adjacencies( *i, 2, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(), mbents.size()); +std::string filename; - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } +int number_tests_successful = 0; +int number_tests_failed = 0; - //NQ2: For every face, obtain neighbor faces - for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { - adjents.clear(); - error = ahf.get_neighbor_adjacencies( *i, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(), mbents.size()); - - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); +void handle_error_code(ErrorCode rv, int &number_failed, int &number_successful) +{ + if (rv == MB_SUCCESS) { + std::cout << "Success"; + number_successful++; + } else { + std::cout << "Failure"; + number_failed++; } +} - // 3D Queries - // IQ 31: For every edge, obtain incident cells - for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { - adjents.clear(); - error = ahf.get_up_adjacencies( *i, 3, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(), mbents.size()); - - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } - //IQ32: For every face, obtain incident cells - for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { - adjents.clear(); - error = ahf.get_up_adjacencies( *i, 3, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); - CHECK_ERR(error); - - CHECK_EQUAL(adjents.size(), mbents.size()); +ErrorCode ahf_test(Core *moab) +{ - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } + Interface* mbImpl = &*moab; + MeshTopoUtil mtu(mbImpl); - //NQ3: For every cell, obtain neighbor cells - for (Range::iterator i = cells.begin(); i != cells.end(); ++i) { - adjents.clear(); - error = ahf.get_neighbor_adjacencies( *i, adjents); - CHECK_ERR(error); - mbents.clear(); - error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents); + ErrorCode error = mbImpl->load_file(filename.c_str()); CHECK_ERR(error); - CHECK_EQUAL(adjents.size(), mbents.size()); - - std::sort(adjents.begin(), adjents.end()); - std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); - mbents = subtract(mbents, ahfents); - CHECK(!mbents.size()); - } - - //ahf.deinitialize(); + /*Create ranges for handles of explicit elements of the mixed mesh*/ + Range verts, edges, faces, cells; + error = mbImpl->get_entities_by_dimension( 0, 0, verts); + error = mbImpl->get_entities_by_dimension( 0, 1, edges); + error = mbImpl->get_entities_by_dimension( 0, 2, faces); + error = mbImpl->get_entities_by_dimension( 0, 3, cells); + + // Create an ahf instance + HalfFacetRep ahf(&*moab); + + // Call the initialize function which creates the maps for each dimension + ahf.initialize(); + + //ahf.print_tags(); + + //Perform queries + std::vector<EntityHandle> adjents; + Range mbents, ahfents; + + //1D Queries // + //IQ1: For every vertex, obtain incident edges + for (Range::iterator i = verts.begin(); i != verts.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 1, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 1, false, mbents ); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(),mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ1: For every edge, obtain neighbor edges + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 0, 1, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + // 2D Queries + //IQ2: For every edge, obtain incident faces + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 2, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies( &*i, 1, 2, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ2: For every face, obtain neighbor faces + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 1, 2, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + // 3D Queries + // IQ 31: For every edge, obtain incident cells + for (Range::iterator i = edges.begin(); i != edges.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 3, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //IQ32: For every face, obtain incident cells + for (Range::iterator i = faces.begin(); i != faces.end(); ++i) { + adjents.clear(); + error = ahf.get_up_adjacencies( *i, 3, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mbImpl->get_adjacencies(&*i, 1, 3, false, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + //NQ3: For every cell, obtain neighbor cells + for (Range::iterator i = cells.begin(); i != cells.end(); ++i) { + adjents.clear(); + error = ahf.get_neighbor_adjacencies( *i, adjents); + CHECK_ERR(error); + mbents.clear(); + error = mtu.get_bridge_adjacencies( *i, 2, 3, mbents); + CHECK_ERR(error); + + CHECK_EQUAL(adjents.size(), mbents.size()); + + std::sort(adjents.begin(), adjents.end()); + std::copy(adjents.begin(), adjents.end(), range_inserter(ahfents)); + mbents = subtract(mbents, ahfents); + CHECK(!mbents.size()); + } + + return MB_SUCCESS; } int main(int argc, char *argv[]) { - int result = 0; - - argv[0] = argv[argc - argc]; // Followed read_mpas_nc.cpp test for removing warnings in serial mode about unused variables. - - result += RUN_TEST(ahf_test); - - return result; + filename = TestDir + "/hexes_mixed.vtk"; + + if (argc==1) + std::cout<<"Using default input file:"<<filename<<std::endl; + else if (argc==2) + filename = argv[1]; + else { + std::cerr << "Usage: " << argv[0] << " [filename]" << std::endl; + return 1; + } + + Core moab; + ErrorCode result; + + std::cout<<" ahf_test: "; + result = ahf_test(&moab); + handle_error_code(result, number_tests_failed, number_tests_successful); + std::cout<<"\n"; + + return number_tests_failed; } Repository URL: https://bitbucket.org/fathomteam/moab/ -- This is a commit notification from bitbucket.org. You are receiving this because you have the service enabled, addressing the recipient of this email.

1 0