April 2010 - moab-dev - mailman3.cels.anl.gov

Re: [MOAB-dev] create ray_intersect_all_sets function
by hongjun＠mcs.anl.gov 06 Apr '10

06 Apr '10

There might be performance differences because multiple calling of "ray_intersect_sets" does tree searches and some triangle intersection checks multiple times. Moreover, for complex geometry, I think we don't know which surfaces will have multiple hits and where we have to call "ray_intersect_sets" again. Hong-Jun ----- Original Message ----- From: "Paul Wilson" <wilsonp(a)engr.wisc.edu> To: "Hongjun Kim" <hongjun(a)mcs.anl.gov> Cc: "Jason Kraftcheck" <kraftche(a)cae.wisc.edu>, moab-dev(a)mcs.anl.gov Sent: Monday, April 5, 2010 10:32:14 PM GMT -06:00 US/Canada Central Subject: Re: [MOAB-dev] create ray_intersect_all_sets function Hi all, I am not sure how this is different from calling ray_intersect_sets multiple times, one after the other? This would march along the ray and find all the intersections. Is it just a performance issue? Paul Hongjun Kim wrote: Your patch is working for me. Are you going to commit it or do you want me to do it? Thanks. Hong-Jun ----- Original Message ----- From: "Jason Kraftcheck" <kraftche(a)cae.wisc.edu> To: hongjun(a)mcs.anl.gov Cc: "Brandon Smith" <bsmith82(a)yahoo.com> , moab-dev(a)mcs.anl.gov Sent: Monday, April 5, 2010 3:00:14 PM GMT -06:00 US/Canada Central Subject: Re: [MOAB-dev] create ray_intersect_all_sets function hongjun(a)mcs.anl.gov wrote: To find all surface sets intersected (including multiple hits for one surface set), I have to give large tolerance. If it is not, it just find the first intersections for each surface set. Specifying a large tolerance to obtain all intersections is a suboptimal solution. You might as well just test all of the triangles against the ray. However, rather than adding 200 lines of duplicated functionality, why not just fix the existing code such that it can return multiple intersections. For example, see the attached patch. - jason Index: src/OrientedBoxTreeTool.cpp =================================================================== --- src/OrientedBoxTreeTool.cpp (revision 3740) +++ src/OrientedBoxTreeTool.cpp (working copy) @@ -950,6 +950,14 @@ void RayIntersectSets::add_intersection( double t, EntityHandle facet ) { + // If minTolInt is less than zero, return all intersections + if (minTolInt < 0 && t > -tol) { + intersections.push_back(t); + sets.push_back(lastSet); + facets.push_back(facet); + return; + } + // Check if the 'len' pointer is pointing into the intersection // list. If this is the case, then the list contains, at that // location, an intersection greater than the tolerance away from Index: src/moab/OrientedBoxTreeTool.hpp =================================================================== --- src/moab/OrientedBoxTreeTool.hpp (revision 3740) +++ src/moab/OrientedBoxTreeTool.hpp (working copy) @@ -233,6 +233,9 @@ * ray start point is less than this number, the next closest * intersection. If the desired result is only the closest * intersection, pass zero for this argument. + * This function will return all intersections, regardless + * of distance from the start of the ray, if this value + * is negative. *\param tolerance The tolerance to use in intersection checks. *\param ray_point The base point of the ray. *\param unit_ray_dir The ray direction vector (must be unit length) -- Paul Wilson -- ------------------------------------------------------------------ -- Paul P.H. Wilson 419 Engineering Research Building wilsonp(a)engr.wisc.edu 1500 Engineering Dr ph/fax: 608/263-0807 Madison, WI 53706 My calendar: http://bit.ly/pphw-calendar Computational Nuclear Engineering Research Group http://cnerg.engr.wisc.edu Associate Professor, Nuclear Engineering Engineering Physics Department http://www.engr.wisc.edu/ep Chair, Energy Analysis & Policy Certificate Nelson Institute for Environmental Studies http://nelson.wisc.edu/eap/ Contributing to the joy and improvement of all those around me

1 0

Re: [MOAB-dev] create ray_intersect_all_sets function
by hongjun＠mcs.anl.gov 05 Apr '10

05 Apr '10

> > I think my way of communication looks wrong. Sorry. I think I need > more > > explanations here as follows. > > > > To find all surface intersections to geomtry (not just only one > > intersection to geometry), > > What are "surface intersections to [geometry]"? I mean I want to find all sets corresponding to the geometric surfaces containing the intersected triangles (including multiple hits). > > > current "ray_intersect_sets" needs large > > distance tolerance to cover all geometries. The large tolerance > makes it > > take long time to find all surface intersections. (several minutes > for > > simple geometry) > > > > Why does it need a large tolerance? To find all surface sets intersected (including multiple hits for one surface set), I have to give large tolerance. If it is not, it just find the first intersections for each surface set. > > > > My version is similar to the function "ray_intersect_triangles" as > > follows. It does box intersection test at first to store > intersected > > boxes and surface sets. After that, it checks > > "GeomUtil::ray_tri_intersect" to do intersection check with > triangles in > > stored boxes. > > But that looses all tracking of geometry subtrees. If all you want is > the > intersect triangles, why not just use ray_intersect_triangles? What > "sets" > does your "ray_intersect_all_sets" function pass back? The tree > leaves? My version stores intersected boxes and surface sets not to lose geometry subtree tracking. Only difference with the current version is just stores the intersected boxes and surface sets, and then do "GeomUtil::ray_tri_intersect" check multiple boxes and surface sets. I need intersected distances, triangles and sets like in "ray_intersect_sets". As mentioned above, "sets" is meshset and it is the meshset corresponding to the geometric surfaces usually. > > > However, current "ray_intersect_sets" searches all leaf > > nodes within large tolerance and then do > "GeomUtil::ray_tri_intersect" > > check for almost all triangles. > > Yes, if you give it an extremely large tolerance, then by definition > the ray > will intersect almost everything. Why are you specifying a large > tolerance? To find all intersected geometry surface sets (including multiple hits), I have to give large tolerance. Current "ray_intersect_sets" looks it can not find multiple hits to one surface. I also attached my code. Thanks. Hong-Jun

4 3

r3744 - in MOAB/trunk/src: . moab
by kraftche＠cae.wisc.edu 05 Apr '10

05 Apr '10

Author: kraftche Date: 2010-04-05 17:23:20 -0500 (Mon, 05 Apr 2010) New Revision: 3744 Modified: MOAB/trunk/src/OrientedBoxTreeTool.cpp MOAB/trunk/src/moab/OrientedBoxTreeTool.hpp Log: get all ray/surface intersections Modified: MOAB/trunk/src/OrientedBoxTreeTool.cpp =================================================================== --- MOAB/trunk/src/OrientedBoxTreeTool.cpp 2010-04-05 21:29:15 UTC (rev 3743) +++ MOAB/trunk/src/OrientedBoxTreeTool.cpp 2010-04-05 22:23:20 UTC (rev 3744) @@ -950,6 +950,14 @@ void RayIntersectSets::add_intersection( double t, EntityHandle facet ) { + // If minTolInt is less than zero, return all intersections + if (minTolInt < 0 && t > -tol) { + intersections.push_back(t); + sets.push_back(lastSet); + facets.push_back(facet); + return; + } + // Check if the 'len' pointer is pointing into the intersection // list. If this is the case, then the list contains, at that // location, an intersection greater than the tolerance away from Modified: MOAB/trunk/src/moab/OrientedBoxTreeTool.hpp =================================================================== --- MOAB/trunk/src/moab/OrientedBoxTreeTool.hpp 2010-04-05 21:29:15 UTC (rev 3743) +++ MOAB/trunk/src/moab/OrientedBoxTreeTool.hpp 2010-04-05 22:23:20 UTC (rev 3744) @@ -233,6 +233,9 @@ * ray start point is less than this number, the next closest * intersection. If the desired result is only the closest * intersection, pass zero for this argument. + * This function will return all intersections, regardless + * of distance from the start of the ray, if this value + * is negative. *\param tolerance The tolerance to use in intersection checks. *\param ray_point The base point of the ray. *\param unit_ray_dir The ray direction vector (must be unit length)

1 0

r3743 - MOAB/trunk
by karpeev＠mcs.anl.gov 05 Apr '10

05 Apr '10

Author: karpeev Date: 2010-04-05 16:29:15 -0500 (Mon, 05 Apr 2010) New Revision: 3743 Modified: MOAB/trunk/configure.ac Log: Added mbparread support. Modified: MOAB/trunk/configure.ac =================================================================== --- MOAB/trunk/configure.ac 2010-04-05 20:37:28 UTC (rev 3742) +++ MOAB/trunk/configure.ac 2010-04-05 21:29:15 UTC (rev 3743) @@ -480,6 +480,7 @@ MB_OPTIONAL_TOOL([gsets], [yes]) MB_OPTIONAL_TOOL([mbdepth], [yes]) MB_OPTIONAL_TOOL([refiner], [no]) +MB_OPTIONAL_TOOL([mbparread], [yes]) if test "xyes" = "x$ENABLE_refiner"; then if test "xyes" != "x$WITH_MPI"; then

1 0

r3742 - MOAB/trunk/tools
by karpeev＠mcs.anl.gov 05 Apr '10

05 Apr '10

Author: karpeev Date: 2010-04-05 15:37:28 -0500 (Mon, 05 Apr 2010) New Revision: 3742 Added: MOAB/trunk/tools/parread.cpp Modified: MOAB/trunk/tools/Makefile.am Log: Added a simple parallel read tool that reports read times. Useful for benchmarking partitioned files. Modified: MOAB/trunk/tools/Makefile.am =================================================================== --- MOAB/trunk/tools/Makefile.am 2010-04-05 17:08:45 UTC (rev 3741) +++ MOAB/trunk/tools/Makefile.am 2010-04-05 20:37:28 UTC (rev 3742) @@ -113,6 +113,10 @@ dagmc_DIR = endif +if ENABLE_mbparread + bin_PROGRAMS += mbparread +endif + SUBDIRS = $(refiner_DIR) \ $(mbcoupler_DIR) \ $(mbzoltan_DIR) \ @@ -121,6 +125,9 @@ $(qvdual_DIR) \ $(dagmc_DIR) +mbparread_SOURCES = parread.cpp +mbparread_DEPENDENCIES = $(LDADD) + mbconvert_SOURCES = convert.cpp mbconvert_DEPENDENCIES = $(LDADD) Added: MOAB/trunk/tools/parread.cpp =================================================================== --- MOAB/trunk/tools/parread.cpp (rev 0) +++ MOAB/trunk/tools/parread.cpp 2010-04-05 20:37:28 UTC (rev 3742) @@ -0,0 +1,365 @@ +/** + * MOAB, a Mesh-Oriented datABase, is a software component for creating, + * storing and accessing finite element mesh data. + * + * Copyright 2004 Sandia Corporation. Under the terms of Contract + * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government + * retains certain rights in this software. + * + * This library is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public + * License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * + */ + +// If Microsoft compiler, then WIN32 +#if defined(_MSC_VER) || defined(__MINGW32__) +# ifndef WIN32 +# define WIN32 +# endif

1 0

r3741 - MOAB/trunk/src
by kraftche＠cae.wisc.edu 05 Apr '10

05 Apr '10

Author: kraftche Date: 2010-04-05 12:08:45 -0500 (Mon, 05 Apr 2010) New Revision: 3741 Modified: MOAB/trunk/src/AdaptiveKDTree.cpp Log: More tweaks to kD-tree code: o Fix possible inadvertent increase in tolerance value used for ray intersections as tree is descended. o Increase tolerance used in tree construction even more, because while false positives might result in a slightly less efficient tree, false negatives result in incorrect trees. Modified: MOAB/trunk/src/AdaptiveKDTree.cpp =================================================================== --- MOAB/trunk/src/AdaptiveKDTree.cpp 2010-04-02 19:15:01 UTC (rev 3740) +++ MOAB/trunk/src/AdaptiveKDTree.cpp 2010-04-05 17:08:45 UTC (rev 3741) @@ -797,7 +797,8 @@ const CartVect left_dim = 0.5*(left_max - box_min); const CartVect right_cen = 0.5*(box_max + right_min); const CartVect right_dim = 0.5*(box_max - right_min); - //const CartVect dim = box_max - box_min; + const CartVect dim = box_max - box_min; + const double max_tol = std::max(dim[0], std::max(dim[1], dim[2]))/10; // test each entity @@ -856,9 +857,9 @@ // sides of plane, do more precise test to ensure that it is really // in both. if (lo && ro) { - double tol = eps/1000; + double tol = eps; lo = ro = false; - while (!lo && !ro && tol <= eps) { + while (!lo && !ro && tol <= max_tol) { lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen, left_dim+CartVect(tol) ); ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i),right_cen,right_dim+CartVect(tol) ); tol *= 10.0; @@ -2091,7 +2092,18 @@ // Otherwise we want the right child (index == 1); list.push_back( NodeSeg( children[fwd_child], seg.beg, seg.end ) ); } - // Otherwise we must intersect the plane + // Otherwise we must intersect the plane. + // Note: be careful not to grow the segment if t is slightly + // outside the current segment, as doing so would effectively + // increase the tolerance as we descend the tree. + else if (t <= seg.beg) { + list.push_back( NodeSeg( children[1-fwd_child], seg.beg, seg.beg ) ); + list.push_back( NodeSeg( children[ fwd_child], seg.beg, seg.end ) ); + } + else if (t >= seg.end) { + list.push_back( NodeSeg( children[1-fwd_child], seg.beg, seg.end ) ); + list.push_back( NodeSeg( children[ fwd_child], seg.end, seg.end ) ); + } else { list.push_back( NodeSeg( children[1-fwd_child], seg.beg, t ) ); list.push_back( NodeSeg( children[ fwd_child], t, seg.end ) );

1 0

Re: [MOAB-dev] create ray_intersect_all_sets function
by hongjun＠mcs.anl.gov 05 Apr '10

05 Apr '10

I think my way of communication looks wrong. Sorry. I think I need more explanations here as follows. To find all surface intersections to geomtry (not just only one intersection to geometry), current "ray_intersect_sets" needs large distance tolerance to cover all geometries. The large tolerance makes it take long time to find all surface intersections. (several minutes for simple geometry) My version is similar to the function "ray_intersect_triangles" as follows. It does box intersection test at first to store intersected boxes and surface sets. After that, it checks "GeomUtil::ray_tri_intersect" to do intersection check with triangles in stored boxes. However, current "ray_intersect_sets" searches all leaf nodes within large tolerance and then do "GeomUtil::ray_tri_intersect" check for almost all triangles. I think current "ray_intersect_sets" function is best to find the first intersection set but not to find all intersection sets. That is one reason I made it as a new seperate function. ---------------------------------------------------------------------- ErrorCode OrientedBoxTreeTool::ray_intersect_all_sets( std::vector<double>& distances_out, std::vector<EntityHandle>& sets_out, std::vector<EntityHandle>& facets_out, EntityHandle root_set, double tolerance, unsigned min_tolerace_intersections, const double ray_point[3], const double unit_ray_dir[3], const double* ray_length, TrvStats* accum ) { Range boxes, box_sets; ErrorCode rval; rval = ray_intersect_boxes_sets( boxes, box_sets, root_set, tolerance, ray_point, unit_ray_dir, ray_length, accum ); if (MB_SUCCESS != rval) return rval; return ray_intersect_all_sets( distances_out, sets_out, facets_out, boxes, box_sets, tolerance, ray_point, unit_ray_dir, ray_length ); RayIntersectAllSets op( this, ray_point, unit_ray_dir, ray_length, tolerance, boxes, box_sets ); return preorder_traverse( root_set, op, accum ); } -------------------------------------------------------- > Why is your version faster? How does the returned data/functionality > differ > from the existing function? Did you look at the existing function in > a > profiler to see if it could easily be made faster? > > > To minimize possible errors with other codes, I made it as a new > function > > not modyfing existing "MBOrientedBoxTreeTool::ray_intersect_sets" > > function. > > > > Having two functions that do almost the same thing will result in a > confusing API. > > > Is it ok to all for me to commit it? Thanks. > > > > Does it do something significantly different than the existing > function? > Have you documented it in the header? Did you add any unit tests for > it? > > - jason

2 1

create ray_intersect_all_sets function
by hongjun＠mcs.anl.gov 05 Apr '10

05 Apr '10

Hi! All I can find all intersected surface sets with "MBOrientedBoxTreeTool::ray_intersect_sets" function and large tolerance, but it takes too long time. So, I made "MBOrientedBoxTreeTool::ray_intersect_all_sets" function to find intersected boxes at first and find all surface sets and triangles intersected with small tolerance like "MBOrientedBoxTreeTool::ray_intersect_triangles" function which is much faster. To minimize possible errors with other codes, I made it as a new function not modyfing existing "MBOrientedBoxTreeTool::ray_intersect_sets" function. Is it ok to all for me to commit it? Thanks. Hong-Jun

2 1

Problem with mbparallelcomm_test?
by Dmitry Karpeev 02 Apr '10

02 Apr '10

I ran mbconvert and mbparallelcomm_test on blocks-64.h5m and tjunc6RIB.h5m using 4 procs on cosmea. It appears that mbconvert has no problem reading these files with PARALLEL=BCAST, but mbparallelcomm_test isn't doing so well. I'm not sure what the difference between the two is (I suppose I should go and sift through the source code), but perhaps mbconvert doesn't have to resolve shared entities, while mbparallelcomm_test does? The mesh files are in ~karpeev/fathom/moab/data/blocks-64.h5m and ~karpeev/fathom/moab/data/tjunc6/tjunc6RIB_16.h5m on MCS machines or on cosmea. Any ideas about what's going on? Thanks. Dmitry. Here's the output: ============================================================================================================================ blocks-64/mbconvert ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --- Running LD_PRELOAD= mpiexec -machinefile mpi.nodes.17204 -n 4 /home/karpeev/fathom/moab/data/../unstable/build/tools/mbconvert '-O' 'PARALLEL=BCAST' '-O' 'PARTITION=MATERIAL_SET' '-f' 'VTK' 'blocks-64.h5m' '/dev/null' 2>&1 1> /home/karpeev/fathom/moab/data/./mbconvert_runs/bcast/blocks-64.4.out Read "blocks-64.h5m" Wrote "/dev/null" Read "blocks-64.h5m" Wrote "/dev/null" Read "Read "blocks-64.h5m" blocks-64.h5m" Wrote "/dev/null" Wrote "/dev/null" --- Done ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- blocks-64/mbparallelcomm_test ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --- Running LD_PRELOAD= mpiexec -machinefile mpi.nodes.17589 -n 4 /home/karpeev/fathom/moab/data/../unstable/build/test/parallel/mbparallelcomm_test '-3' '3' '1' 'tjunc6/blocks-64.h5m' 'MATERIAL_SET' 2>&1 1> /home/karpeev/fathom/moab/data/./mbparallelcomm_test_runs/bcast/blocks-64.4.out Couldn't read mesh; error message: Failed in step PARALLEL READ application called MPI_Abort(MPI_COMM_WORLD, 0) - process 0--- Done ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tjunc6RIB_16/mbconvert ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --- Running LD_PRELOAD= mpiexec -machinefile mpi.nodes.17340 -n 4 /home/karpeev/fathom/moab/data/../unstable/build/tools/mbconvert '-O' 'PARALLEL=BCAST' '-O' 'PARTITION=PARALLEL_PARTITION' '-f' 'VTK' 'tjunc6/tjunc6RIB_16.h5m' '/dev/null' 2>&1 1> /home/karpeev/fathom/moab/data/./mbconvert_runs/bcast/tjunc6RIB_16.4.out Read "tjunc6/tjunc6RIB_16.h5m" Read "tjunc6/tjunc6RIB_16.h5m" Read "tjunc6/tjunc6RIB_16.h5m" Read "tjunc6/tjunc6RIB_16.h5m" Wrote "/dev/null" Wrote "/dev/null" Wrote "/dev/null" Wrote "/dev/null" --- Done ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- tjunc6RIB_16/mbparallelcomm_test ------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- --- Running LD_PRELOAD= mpiexec -machinefile mpi.nodes.17453 -n 4 /home/karpeev/fathom/moab/data/../unstable/build/test/parallel/mbparallelcomm_test '-3' '3' '1' 'tjunc6/tjunc6RIB_16.h5m' 'PARALLEL_PARTITION' 2>&1 1> /home/karpeev/fathom/moab/data/./mbparallelcomm_test_runs/bcast/tjunc6RIB_16.4.out Couldn't read mesh; error message: Failed in step PARALLEL RESOLVE_SHARED_ENTS application called MPI_Abort(MPI_COMM_WORLD, 0) - process 0--- Done

2 2

r3740 - MOAB/trunk/src
by kraftche＠cae.wisc.edu 02 Apr '10

02 Apr '10

Author: kraftche Date: 2010-04-02 14:15:01 -0500 (Fri, 02 Apr 2010) New Revision: 3740 Modified: MOAB/trunk/src/AdaptiveKDTree.cpp Log: o Rewrite ray/node intersection in KDTree to be a) correct, b) faster, c) shorter and hopefully easier to understand, and d) better commented. o Change odd tolerance stuff in tree construction that results in invalid trees for cases where splits are parallel to triangles. Modified: MOAB/trunk/src/AdaptiveKDTree.cpp =================================================================== --- MOAB/trunk/src/AdaptiveKDTree.cpp 2010-04-01 22:14:56 UTC (rev 3739) +++ MOAB/trunk/src/AdaptiveKDTree.cpp 2010-04-02 19:15:01 UTC (rev 3740) @@ -23,6 +23,7 @@ #include "moab/GeomUtil.hpp" #include "moab/Range.hpp" #include "Internals.hpp" +#include <math.h> #include <assert.h> #include <algorithm> @@ -796,7 +797,7 @@ const CartVect left_dim = 0.5*(left_max - box_min); const CartVect right_cen = 0.5*(box_max + right_min); const CartVect right_dim = 0.5*(box_max - right_min); - const CartVect dim = box_max - box_min; + //const CartVect dim = box_max - box_min; // test each entity @@ -850,18 +851,16 @@ ro = true; } - // Triangle must be in at least one leaf. If test against plain + // Triangle must be in at least one leaf. If test against plane // identified that leaf, then we're done. If triangle is on both // sides of plane, do more precise test to ensure that it is really // in both. if (lo && ro) { - lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen, left_dim ); - ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i),right_cen,right_dim ); - double tol = std::numeric_limits<double>::epsilon(); - // didn't intersect either - tolerance issue - while (!lo && !ro && tol < eps) { - lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen,left_dim+tol*dim ); - ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i),right_cen,right_dim+tol*dim ); + double tol = eps/1000; + lo = ro = false; + while (!lo && !ro && tol <= eps) { + lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen, left_dim+CartVect(tol) ); + ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i),right_cen,right_dim+CartVect(tol) ); tol *= 10.0; } } @@ -1949,13 +1948,13 @@ }; ErrorCode AdaptiveKDTree::ray_intersect_triangles( EntityHandle root, - const double tol,

1 0