[MOAB-dev] commit/MOAB: 11 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Oct 2 11:38:42 CDT 2013
11 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/5012686b4467/
Changeset: 5012686b4467
Branch: None
User: tautges
Date: 2013-09-12 23:29:44
Summary: Merged fathomteam/moab into master
Affected #: 81 files
diff --git a/.gitignore b/.gitignore
index 967aaa7..d3721c3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,99 +1,132 @@
-MOABConfig.cmake
-moab.files
-config/test-driver
-moab.config
-moab.includes
-moab.creator*
-*.ccm
-*.cub
+*~
+*.a
aclocal.m4
autom4te.cache/
-config.h
-config.h.in
-config.log
-config.lt
-config.status
+bin
+bin/*
+*.ccm
config/config.guess
config/config.sub
config/depcomp
+config.h
+config.h.in
config/install-sh
config/libtool.m4
+config.log
+config.lt
config/ltmain.sh
+config/lt~obsolete.m4
config/ltoptions.m4
config/ltsugar.m4
config/ltversion.m4
-config/lt~obsolete.m4
config/missing
+config.status
+config/test-driver
configure
+.cproject
+*.cub
+.deps
doc/config.tex
doc/dev.dox
-doc/user.dox
doc/user/*
+doc/user.dox
examples/examples.make
+examples/FileRead
+examples/GeomSetHierarchy
+examples/GetEntities
+examples/*.h5m
+examples/HelloMoabPar
+examples/itaps/FindConnectF
+examples/itaps/ListSetsNTagsCXX
+examples/itaps/ListSetsNTagsF90
+examples/itaps/TagIterateC
+examples/itaps/TagIterateF
+examples/KDTree
+examples/ObbTree
+examples/ReduceExchangeTags
+examples/SetsNTags
+examples/SkinMesh
+examples/SurfArea
+examples/TestExodusII
+include
+include/*
itaps/iBase_f.h
itaps/igeom/FBiGeom-Defs.inc
+itaps/igeom/FBiGeom_protos.h
+itaps/igeom/testgeom
+itaps/igeom/testSmooth2
+itaps/igeom/testSmoothGeom
+itaps/imesh/FindAdjacencyF90
itaps/imesh/iMesh-Defs.inc
+itaps/imesh/iMesh_extensions_protos.h
itaps/imesh/iMeshP_extensions_protos.h
itaps/imesh/iMeshP_protos.h
-itaps/imesh/iMesh_extensions_protos.h
itaps/imesh/iMesh_protos.h
+itaps/imesh/MOAB_iMesh_extensions_tests
+itaps/imesh/MOAB_iMeshP_unit_tests
+itaps/imesh/MOAB_iMesh_unit_tests
+itaps/imesh/partest
+itaps/imesh/ScdMeshF77
+itaps/imesh/ScdMeshF90
+itaps/imesh/testc_cbind
+*.la
+*.la
+*.lai
+lib
+lib/*
+.libs
libtool
+*.lo
+*.log
+makefile
+Makefile
+*/Makefile
+*/**/Makefile
+*/**/Makefile.in
+*/Makefile.in
+Makefile.in
+moab.config
+MOABConfig.cmake
+moab.creator*
+moab.files
+moab.includes
moab.make
+*.o
+.project
+*.rej
+share/*
+share/doc/moab
+share/man/man1
+*.so
src/FCDefs.h
+src/io/mhdf/h5minfo
+src/io/mhdf/h5mvalidate
src/MBCN_protos.h
-src/MOAB_FCDefs.h
src/moab/EntityHandle.hpp
-src/moab/Version.h
+src/MOAB_FCDefs.h
src/moab/stamp-h2
src/moab/stamp-h3
+src/moab/Version.h
src/parallel/moab_mpi_config.h
src/parallel/stamp-h4
src/stamp-h5
stamp-h1
-tools/mbcoupler/tests/
-tools/mbzoltan/Config.moab
-tools/vtkMOABReader/CMakeLists.txt
-tools/vtkMOABReaderNew/CMakeLists.txt
-.deps
-Makefile.in
-Makefile
-*/Makefile.in
-*/Makefile
-*/**/Makefile.in
-*/**/Makefile
-.libs
-*.o
-*.log
-*.lo
-*.la
-*.a
-*.so
-include/*
-lib/*
-share/*
-bin/*
-*~
-examples/HelloMoabPar
-examples/TestExodusII
-itaps/igeom/FBiGeom_protos.h
-itaps/igeom/testSmooth2
-itaps/igeom/testSmoothGeom
-itaps/igeom/testgeom
-itaps/imesh/FindAdjacencyF90
-itaps/imesh/MOAB_iMeshP_unit_tests
-itaps/imesh/ScdMeshF77
-itaps/imesh/ScdMeshF90
-itaps/imesh/partest
test/adaptive_kd_tree_tests
test/bsp_tree_poly_test
test/bsp_tree_test
-test/*.gen
+test/CMakeLists.txt
test/coords_connect_iterate
test/cropvol_test
test/dual/dual_test
+test/elem_eval_test
test/file_options_test
+test/*.g
+test/*.g
+test/*.gen
+test/*.gen
test/geom_util_test
test/gttool_test
+test/gttool_test
test/h5file/dump_sets
test/h5file/h5legacy
test/h5file/h5partial
@@ -102,13 +135,14 @@ test/h5file/h5regression
test/h5file/h5sets_test
test/h5file/h5test
test/h5file/h5varlen
-test/*.g
test/homxform_test
+test/io/*.ccmg
+test/io/ccmio_test
test/io/cub_file_test
test/io/exodus_test
+test/io/*.g
test/io/gmsh_test
test/io/ideas_test
-test/io/*.g
test/io/nastran_test
test/io/read_cgm_test
test/io/read_nc
@@ -133,7 +167,6 @@ test/obb_test
test/oldinc/test_oldinc
test/read_mpas_nc
test/parallel/*.h5m
-test/parallel/*.vtk
test/parallel/mbparallelcomm_test
test/parallel/mhdf_parallel
test/parallel/mpastrvpart
@@ -142,6 +175,8 @@ test/parallel/par_intx_sph
test/parallel/parallel_hdf5_test
test/parallel/parallel_unit_tests
test/parallel/parallel_write_test
+test/parallel/par_coupler_test
+test/parallel/par_intx_sph
test/parallel/parmerge
test/parallel/partcheck
test/parallel/pcomm_serial
@@ -152,39 +187,68 @@ test/parallel/scdtest
test/parallel/structured3
test/parallel/uber_parallel_test
test/parallel/ucdtrvpart
+test/parallel/*.vtk
test/perf/adj_time
test/perf/perf
test/perf/perftool
+test/perf/point_in_elem
+test/perf/runtest
test/perf/seqperf
test/perf/tstt_perf_binding
+test/perf/point_location/elem_eval_time
+test/perf/point_location/point_location
test/range_test
test/reorder_test
test/scdseq_test
test/scd_test_partn
test/seq_man_test
+test/spatial_locator_test
+test/test_boundbox
test/tag_test
test/test_adj
test/test_prog_opt
+test/TestRunner.hpp
test/var_len_test
test/var_len_test_no_template
test/xform_test
+tools/dagmc/dagmc_preproc
tools/dagmc/pt_vol_test
+tools/dagmc/quads_to_tris
tools/dagmc/ray_fire_test
tools/dagmc/test_geom
tools/dagmc/update_coords
+tools/hexmodops
+tools/mbconvert
tools/mbcoupler/*.h5m
+tools/mbcoupler/*.g
+tools/mbcoupler/tests/
tools/mbcslam/case1_test
tools/mbcslam/intersect1.h5m
-tools/mbcslam/intx.vtk
tools/mbcslam/intx1.vtk
tools/mbcslam/intx_in_plane_test
tools/mbcslam/intx_on_sphere_test
+tools/mbcslam/intx.vtk
tools/mbcslam/lagr.h5m
-tools/mbcslam/spec_visu_test
tools/mbcslam/spectral.vtk
+tools/mbcslam/spec_visu_test
tools/mbcslam/spherical_area_test
-.project
-.cproject
-examples/*.h5m
-examples/ReduceExchangeTags
-
+tools/mbdepth
+tools/mbgsets
+tools/mbmem
+tools/mbsize
+tools/mbskin
+tools/mbsurfplot
+tools/mbtagprop
+tools/mbzoltan/Config.moab
+tools/spheredecomp
+tools/vtkMOABReader/CMakeLists.txt
+tools/vtkMOABReaderNew/CMakeLists.txt
+tools/vtkMOABReaderNew/CMakeCache.txt
+tools/vtkMOABReaderNew/CMakeFiles/*
+tools/vtkMOABReaderNew/cmake_install.cmake
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin.qrc
+tools/vtkMOABReaderNew/vtkMoabReaderPluginInit.cxx
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin_Plugin.cxx
+tools/vtkMOABReaderNew/vtkMoabReaderPlugin_Plugin.h
+tools/vtkMOABReaderNew/vtkSMvtkMoabReaderPluginInstantiator.cxx
+tools/vtkMOABReaderNew/vtkSMvtkMoabReaderPluginInstantiator.h
diff --git a/itaps/igeom/FBiGeom_MOAB.cpp b/itaps/igeom/FBiGeom_MOAB.cpp
index cb889bd..7fcd2b1 100644
--- a/itaps/igeom/FBiGeom_MOAB.cpp
+++ b/itaps/igeom/FBiGeom_MOAB.cpp
@@ -4,7 +4,7 @@
#include "moab/GeomTopoTool.hpp"
#include "moab/OrientedBoxTreeTool.hpp"
#include "moab/CartVect.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "MBTagConventions.hpp"
#include <stdlib.h>
#include <cstring>
diff --git a/itaps/imesh/iMeshP_MOAB.cpp b/itaps/imesh/iMeshP_MOAB.cpp
index e554fba..51254ca 100644
--- a/itaps/imesh/iMeshP_MOAB.cpp
+++ b/itaps/imesh/iMeshP_MOAB.cpp
@@ -4,7 +4,7 @@
#include "moab/Range.hpp"
#include "moab/CN.hpp"
#include "moab/MeshTopoUtil.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/ParallelComm.hpp"
#include "MBParallelConventions.h"
#include "MBIter.hpp"
diff --git a/itaps/imesh/iMesh_MOAB.cpp b/itaps/imesh/iMesh_MOAB.cpp
index d09fec1..68c0b72 100644
--- a/itaps/imesh/iMesh_MOAB.cpp
+++ b/itaps/imesh/iMesh_MOAB.cpp
@@ -4,7 +4,7 @@
#include "moab/CN.hpp"
#include "moab/MeshTopoUtil.hpp"
#include "moab/ScdInterface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "iMesh_MOAB.hpp"
#include "MBIter.hpp"
#include "MBTagConventions.hpp"
@@ -3303,7 +3303,7 @@ void iMesh_createStructuredMesh(iMesh_Instance instance,
ScdBox *scd_box;
rval = scdi->construct_box(HomCoord(local_dims[0], local_dims[1], (-1 != local_dims[2] ? local_dims[2] : 0), 1),
HomCoord(local_dims[3], local_dims[4], (-1 != local_dims[5] ? local_dims[5] : 0), 1),
- NULL, 0, scd_box);
+ NULL, 0, scd_box, NULL, NULL, (vert_gids ? true : false));
CHKERR(rval, "Trouble creating scd vertex sequence.");
// set the global box parameters
diff --git a/src/Core.cpp b/src/Core.cpp
index 57ff342..476fdca 100644
--- a/src/Core.cpp
+++ b/src/Core.cpp
@@ -83,7 +83,7 @@
#include "MBTagConventions.hpp"
#include "ExoIIUtil.hpp"
#include "EntitySequence.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#ifdef LINUX
# include <dlfcn.h>
# include <dirent.h>
@@ -1153,11 +1153,11 @@ ErrorCode Core::get_connectivity_by_type(const EntityType type,
ErrorCode Core::get_connectivity(const EntityHandle *entity_handles,
const int num_handles,
Range &connectivity,
- bool topological_connectivity) const
+ bool corners_only) const
{
std::vector<EntityHandle> tmp_connect;
ErrorCode result = get_connectivity(entity_handles, num_handles, tmp_connect,
- topological_connectivity);
+ corners_only);
if (MB_SUCCESS != result) return result;
std::sort( tmp_connect.begin(), tmp_connect.end() );
@@ -1169,7 +1169,7 @@ ErrorCode Core::get_connectivity(const EntityHandle *entity_handles,
ErrorCode Core::get_connectivity(const EntityHandle *entity_handles,
const int num_handles,
std::vector<EntityHandle> &connectivity,
- bool topological_connectivity,
+ bool corners_only,
std::vector<int> *offsets) const
{
connectivity.clear(); // this seems wrong as compared to other API functions,
@@ -1182,7 +1182,7 @@ ErrorCode Core::get_connectivity(const EntityHandle *entity_handles,
int len;
if (offsets) offsets->push_back(0);
for (int i = 0; i < num_handles; ++i) {
- rval = get_connectivity( entity_handles[i], conn, len, topological_connectivity, &tmp_storage );
+ rval = get_connectivity( entity_handles[i], conn, len, corners_only, &tmp_storage );
if (MB_SUCCESS != rval)
return rval;
connectivity.insert( connectivity.end(), conn, conn + len );
@@ -1195,7 +1195,7 @@ ErrorCode Core::get_connectivity(const EntityHandle *entity_handles,
ErrorCode Core::get_connectivity(const EntityHandle entity_handle,
const EntityHandle*& connectivity,
int& number_nodes,
- bool topological_connectivity,
+ bool corners_only,
std::vector<EntityHandle>* storage) const
{
ErrorCode status;
@@ -1222,7 +1222,7 @@ ErrorCode Core::get_connectivity(const EntityHandle entity_handle,
return static_cast<const ElementSequence*>(seq)->get_connectivity(entity_handle,
connectivity,
number_nodes,
- topological_connectivity,
+ corners_only,
storage);
}
@@ -2896,9 +2896,11 @@ ErrorCode Core::list_entity(const EntityHandle entity) const
if (multiple != 0)
std::cout << " (MULTIPLE = " << multiple << ")" << std::endl;
+ result = print_entity_tags(std::string(), entity, MB_TAG_DENSE);
+
std::cout << std::endl;
- return MB_SUCCESS;
+ return result;
}
ErrorCode Core::convert_entities( const EntityHandle meshset,
@@ -3598,16 +3600,21 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
}
// print all sparse tags
+ print_entity_tags(indent_prefix, ms_handle, MB_TAG_SPARSE);
+}
+
+ErrorCode Core::print_entity_tags(std::string indent_prefix, const EntityHandle handle, TagType tp) const
+{
std::vector<Tag> set_tags;
- ErrorCode result = this->tag_get_tags_on_entity(ms_handle, set_tags);
- std::cout << indent_prefix << "Sparse tags:" << std::endl;
+ ErrorCode result = this->tag_get_tags_on_entity(handle, set_tags);
+ std::cout << indent_prefix << (tp == MB_TAG_SPARSE ? "Sparse tags:" : "Dense tags:") << std::endl;
indent_prefix += " ";
for (std::vector<Tag>::iterator vit = set_tags.begin();
vit != set_tags.end(); vit++) {
TagType this_type;
result = this->tag_get_type(*vit, this_type);
- if (MB_SUCCESS != result || MB_TAG_SPARSE != this_type) continue;
+ if (MB_SUCCESS != result || tp != this_type) continue;
DataType this_data_type;
result = this->tag_get_data_type(*vit, this_data_type);
int this_size;
@@ -3622,7 +3629,7 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
if (MB_SUCCESS != result) continue;
switch (this_data_type) {
case MB_TYPE_INTEGER:
- result = this->tag_get_data(*vit, &ms_handle, 1, &int_vals[0]);
+ result = this->tag_get_data(*vit, &handle, 1, &int_vals[0]);
if (MB_SUCCESS != result) continue;
std::cout << indent_prefix << tag_name << " = ";
if (this_size < 10)
@@ -3631,7 +3638,7 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
std::cout << std::endl;
break;
case MB_TYPE_DOUBLE:
- result = this->tag_get_data(*vit, &ms_handle, 1, &dbl_vals[0]);
+ result = this->tag_get_data(*vit, &handle, 1, &dbl_vals[0]);
if (MB_SUCCESS != result) continue;
std::cout << indent_prefix << tag_name << " = ";
if (this_size < 10)
@@ -3640,7 +3647,7 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
std::cout << std::endl;
break;
case MB_TYPE_HANDLE:
- result = this->tag_get_data(*vit, &ms_handle, 1, &hdl_vals[0]);
+ result = this->tag_get_data(*vit, &handle, 1, &hdl_vals[0]);
if (MB_SUCCESS != result) continue;
std::cout << indent_prefix << tag_name << " = ";
if (this_size < 10)
@@ -3651,7 +3658,7 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
case MB_TYPE_OPAQUE:
if (NAME_TAG_SIZE == this_size) {
char dum_tag[NAME_TAG_SIZE];
- result = this->tag_get_data(*vit, &ms_handle, 1, &dum_tag);
+ result = this->tag_get_data(*vit, &handle, 1, &dum_tag);
if (MB_SUCCESS != result) continue;
// insert NULL just in case there isn't one
dum_tag[NAME_TAG_SIZE-1] = '\0';
@@ -3662,6 +3669,8 @@ void Core::print(const EntityHandle ms_handle, const char *prefix,
break;
}
}
+
+ return MB_SUCCESS;
}
ErrorCode Core::check_adjacencies()
diff --git a/src/FileOptions.cpp b/src/FileOptions.cpp
index f2a4477..380a491 100644
--- a/src/FileOptions.cpp
+++ b/src/FileOptions.cpp
@@ -18,7 +18,7 @@
*\date 2007-08-21
*/
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <ctype.h>
#include <stdlib.h>
diff --git a/src/FileOptions.hpp b/src/FileOptions.hpp
deleted file mode 100644
index 85b3eef..0000000
--- a/src/FileOptions.hpp
+++ /dev/null
@@ -1,248 +0,0 @@
-/*
- * 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.
- *
- */
-
-/**\file FileOptions.hpp
- *\author Jason Kraftcheck (kraftche at cae.wisc.edu)
- *\date 2007-08-21
- */
-
-#ifndef FILE_OPTIONS_HPP
-#define FILE_OPTIONS_HPP
-
-#include <string>
-#include <vector>
-#include "moab/Types.hpp"
-
-namespace moab {
-
-/**\brief Parse options string passed to file IO routines
- *
- * This is a utility class used by file-IO-related code to
- * parse the options string passed to Core::load_file and
- * Core::write_file
- */
-class FileOptions {
-public:
-
- /*\param options_string The concatenation of a list of
- * options, separated either by the default separator
- * (semicolon) with a custom separator specified at
- * the beginning of the string (semicolon followed by
- * destired separator character.)
- */
- FileOptions( const char* option_string );
-
- FileOptions( const FileOptions& copy );
- FileOptions& operator=( const FileOptions& copy );
-
- ~FileOptions();
-
- /**\brief Check for option with no value
- *
- * Check for an option w/out a value.
- *\param name The option name
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but has value
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_null_option( const char* name ) const;
-
-
- /**\brief Check for option with boolean (true/false, yes/no) value.
- *
- * Check for an option with a true/false value. Allowable values
- * are "true", "false", "yes", "no", "1", "0", "on", "off".
- *\param name The option name
- *\param default_value The value to return if the option is not specified.
- *\param value The resulting value. This argument is set to the passed
- * default value if the option is not found.
- *\return - MB_TYPE_OUT_OF_RANGE if options is found, but has an invalid value
- * - MB_SUCCESS otherwise
- */
- ErrorCode get_toggle_option( const char* name,
- bool default_value,
- bool& value ) const;
-
- /**\brief Check for option with an integer value.
- *
- * Check for an option with an integer value
- *\param name The option name
- *\param value Output. The value.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have an integer value
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_int_option( const char* name, int& value ) const;
-
- /**\brief Check for option with an integer value. Accept option with no value.
- *
- * Check for an option with an integer value.
- * If the option is found but has no value specified, then
- * pass back the user-specified default value.
- *
- *\NOTE: This function will not pass back the default_val, but will instead
- * return MB_ENTITY_NOT_FOUND if the option is not specified at all.
- * The default value is returned only when the option is specified,
- * but is specified w/out a value.
- *
- *\param name The option name
- *\param default_val The default value for the option.
- *\param value Output. The value.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found but has a value that cannot be parsed as an int
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_int_option( const char* name, int default_val, int& value ) const;
-
- /**\brief Check for option with a double value.
- *
- * Check for an option with a double value
- *\param name The option name
- *\param value Output. The value.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have a double value
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_real_option( const char* name, double& value ) const;
-
- /**\brief Check for option with any value.
- *
- * Check for an option with any value.
- *\param name The option name
- *\param value Output. The value.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have a value
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_str_option( const char* name, std::string& value ) const;
-
- /**\brief Check for option
- *
- * Check for an option
- *\param name The option name
- *\param value The option value, or an empty string if no value.
- *\return MB_SUCCESS or MB_ENTITY_NOT_FOUND
- */
- ErrorCode get_option( const char* name, std::string& value ) const;
-
- /**\brief Check the string value of an option
- *
- * Check which of a list of possible values a string option contains.
- *\param name The option name
- *\param values A NULL-terminated array of C-style strings enumerating
- * the possible option values.
- *\param index Output: The index into <code>values</code> for the
- * option value.
- *\return MB_SUCCESS if matched name and value.
- * MB_ENTITY_NOT_FOUND if the option was not specified
- * MB_FAILURE if the option value is not in the input <code>values</code> array.
- */
- ErrorCode match_option( const char* name, const char* const* values, int& index ) const;
-
- /**\brief Check if an option matches a string value
- *
- * Check if the value for an option is the passed string.
- *\param name The option name
- *\param value The expected value.
- *\return MB_SUCCESS if matched name and value.
- * MB_ENTITY_NOT_FOUND if the option was not specified
- * MB_FAILURE if the option value doesn't match the passed string/
- */
- ErrorCode match_option( const char* name, const char* value) const;
-
- /**\brief Check for option for which the value is a list of ints
- *
- * Check for an option which is an int list. The value is expected to
- * be a comma-separated list of int ranges, where an int range can be
- * either a single integer value or a range of integer values separated
- * by a dash ('-').
- *
- *\param name The option name
- *\param values Output. The list of integer values.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain an ID list
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_ints_option( const char* name, std::vector<int>& values) const;
-
- /**\brief Check for option for which the value is a list of doubles
- *
- * Check for an option which is a double list. The value is expected to
- * be a comma-separated list of double values
- *
- *\param name The option name
- *\param values Output. The list of double values.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain an ID list
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_reals_option( const char* name, std::vector<double>& values) const;
-
- /**\brief Check for option for which the value is a list of strings
- *
- * Check for an option which is a string list. The value is expected to
- * be a comma-separated list of string values, with no embedded spaces or commas.
- *
- *\param name The option name
- *\param values Output. The list of string values.
- *\return - MB_SUCCESS if option is found
- * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain a string list
- * - MB_ENTITY_NOT_FOUND if option is not found.
- */
- ErrorCode get_strs_option( const char* name, std::vector<std::string>& values) const;
-
- /** number of options */
- inline unsigned size() const
- { return mOptions.size(); }
-
- /** true if no options */
- inline bool empty() const
- { return mOptions.empty(); }
-
- /** Get list of options */
- void get_options( std::vector<std::string>& list ) const;
-
- /** Check if all options have been looked at */
- bool all_seen() const;
-
- /** Mark all options as seen. USed for non-root procs during bcast-delete read */
- void mark_all_seen() const;
-
- /** Get first unseen option */
- ErrorCode get_unseen_option( std::string& value ) const;
-
-private:
-
- /**\brief Check for option
- *
- * Check for an option
- *\param name The option name
- *\param value The option value, or an empty string if no value.
- *\return MB_SUCCESS or MB_ENTITY_NOT_FOUND
- */
- ErrorCode get_option( const char* name, const char*& value) const;
-
- char* mData;
- std::vector<const char*> mOptions;
- mutable std::vector<bool> mSeen;
-
- /** Case-insensitive compare of name with option value. */
- static bool compare( const char* name, const char* option );
-};
-
-} // namespace moab
-
-#endif
-
diff --git a/src/Makefile.am b/src/Makefile.am
index 7cebb7c..ab5077a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -41,6 +41,7 @@ libMOAB_la_SOURCES = \
BitPage.hpp \
BitTag.cpp \
BitTag.hpp \
+ BoundBox.cpp \
BSPTree.cpp \
BSPTreePoly.cpp \
CN.cpp \
@@ -57,7 +58,6 @@ libMOAB_la_SOURCES = \
Factory.cpp \
FBEngine.cpp \
FileOptions.cpp \
- FileOptions.hpp \
GeomUtil.cpp \
GeomTopoTool.cpp \
HigherOrderFactory.cpp \
@@ -75,16 +75,10 @@ libMOAB_la_SOURCES = \
OrientedBox.cpp \
OrientedBox.hpp \
OrientedBoxTreeTool.cpp \
- moab/point_locater/lotte/minmax.h \
- moab/point_locater/lotte/extrafindpt.h \
lotte/poly.c \
- moab/point_locater/lotte/poly.h \
lotte/findpt.c \
- moab/point_locater/lotte/findpt.h \
lotte/errmem.c \
- moab/point_locater/lotte/errmem.h \
lotte/tensor.c \
- moab/point_locater/lotte/tensor.h \
PolyElementSeq.cpp \
PolyElementSeq.hpp \
ProgOptions.cpp \
@@ -146,12 +140,14 @@ libMOAB_la_SOURCES = \
# The list of header files which are to be installed
nobase_libMOAB_la_include_HEADERS = \
moab/AdaptiveKDTree.hpp \
+ moab/BoundBox.hpp \
moab/BSPTree.hpp \
moab/BSPTreePoly.hpp \
moab/CN.hpp \
moab/CartVect.hpp \
moab/Compiler.hpp \
moab/Core.hpp \
+ moab/CpuTimer.hpp \
moab/DualTool.hpp \
moab/Error.hpp \
moab/GeomTopoTool.hpp \
@@ -160,6 +156,7 @@ nobase_libMOAB_la_include_HEADERS = \
moab/EntityType.hpp \
moab/EntityHandle.hpp \
moab/FBEngine.hpp \
+ moab/FileOptions.hpp \
moab/Forward.hpp \
moab/GeomUtil.hpp \
moab/Interface.hpp \
@@ -172,10 +169,6 @@ nobase_libMOAB_la_include_HEADERS = \
moab/point_locater/element_maps/linear_tet_map.hpp \
moab/point_locater/element_maps/spectral_hex_map.hpp \
moab/point_locater/element_maps/quadratic_hex_map.hpp \
- moab/point_locater/lotte/types.h \
- moab/point_locater/lotte/errmem.h \
- moab/point_locater/lotte/findpt.h \
- moab/point_locater/lotte/poly.h \
moab/point_locater/point_locater.hpp \
moab/point_locater/parametrizer.hpp \
moab/Matrix3.hpp \
diff --git a/src/ScdInterface.cpp b/src/ScdInterface.cpp
index 1441e51..4dd7c5b 100644
--- a/src/ScdInterface.cpp
+++ b/src/ScdInterface.cpp
@@ -145,11 +145,13 @@ ErrorCode ScdInterface::construct_box(HomCoord low, HomCoord high, const double
rval = create_scd_sequence(low, high, MBVERTEX, 0, new_box);
ERRORR(rval, "Trouble creating scd vertex sequence.");
- if (num_coords && coords) {
- // set the vertex coordinates
- double *xc, *yc, *zc;
- rval = new_box->get_coordinate_arrays(xc, yc, zc);
- ERRORR(rval, "Couldn't get vertex coordinate arrays.");
+ // set the vertex coordinates
+ double *xc, *yc, *zc;
+
+ rval = new_box->get_coordinate_arrays(xc, yc, zc);
+ ERRORR(rval, "Couldn't get vertex coordinate arrays.");
+
+ if (coords && num_coords) {
unsigned int i = 0;
for (int kl = low[2]; kl <= high[2]; kl++) {
@@ -165,6 +167,23 @@ ErrorCode ScdInterface::construct_box(HomCoord low, HomCoord high, const double
}
}
}
+ else {
+ unsigned int i = 0;
+ for (int kl = low[2]; kl <= high[2]; kl++) {
+ for (int jl = low[1]; jl <= high[1]; jl++) {
+ for (int il = low[0]; il <= high[0]; il++) {
+ xc[i] = (double) il;
+ if (new_box->box_size()[1])
+ yc[i] = (double) jl;
+ else yc[i] = 0.0;
+ if (new_box->box_size()[2])
+ zc[i] = (double) kl;
+ else zc[i] = 0.0;
+ i++;
+ }
+ }
+ }
+ }
// create element sequence
Core *mbcore = dynamic_cast<Core*>(mbImpl);
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
index e5b3f08..ca86127 100644
--- a/src/io/NCHelperEuler.cpp
+++ b/src/io/NCHelperEuler.cpp
@@ -1,6 +1,6 @@
#include "NCHelperEuler.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <cmath>
#include <sstream>
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
index dd555ff..a234450 100644
--- a/src/io/NCHelperFV.cpp
+++ b/src/io/NCHelperFV.cpp
@@ -1,6 +1,6 @@
#include "NCHelperFV.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <cmath>
#include <sstream>
diff --git a/src/io/NCHelperHOMME.cpp b/src/io/NCHelperHOMME.cpp
index f883bd7..4f2ac51 100644
--- a/src/io/NCHelperHOMME.cpp
+++ b/src/io/NCHelperHOMME.cpp
@@ -1,6 +1,6 @@
#include "NCHelperHOMME.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/SpectralMeshTool.hpp"
#include <cmath>
diff --git a/src/io/NCHelperMPAS.cpp b/src/io/NCHelperMPAS.cpp
index 45d1b00..c0605e8 100644
--- a/src/io/NCHelperMPAS.cpp
+++ b/src/io/NCHelperMPAS.cpp
@@ -1,6 +1,6 @@
#include "NCHelperMPAS.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/SpectralMeshTool.hpp"
#include "MBTagConventions.hpp"
diff --git a/src/io/ReadABAQUS.cpp b/src/io/ReadABAQUS.cpp
index 1c44057..b3a37d0 100644
--- a/src/io/ReadABAQUS.cpp
+++ b/src/io/ReadABAQUS.cpp
@@ -36,7 +36,7 @@
#include "moab/ReadUtilIface.hpp"
#include "AffineXform.hpp"
// #include "abaqus_order.h"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
namespace moab {
diff --git a/src/io/ReadCCMIO.cpp b/src/io/ReadCCMIO.cpp
index 3161425..ec53b6c 100644
--- a/src/io/ReadCCMIO.cpp
+++ b/src/io/ReadCCMIO.cpp
@@ -11,7 +11,7 @@
#include "MBTagConventions.hpp"
#include "Internals.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "ReadCCMIO.hpp"
#include "moab/MeshTopoUtil.hpp"
diff --git a/src/io/ReadCGM.cpp b/src/io/ReadCGM.cpp
index 02f748e..a610d4a 100644
--- a/src/io/ReadCGM.cpp
+++ b/src/io/ReadCGM.cpp
@@ -37,7 +37,7 @@
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "MBTagConventions.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/GeomTopoTool.hpp"
diff --git a/src/io/ReadDamsel.cpp b/src/io/ReadDamsel.cpp
index 8d1bcfe..f02f0ea 100644
--- a/src/io/ReadDamsel.cpp
+++ b/src/io/ReadDamsel.cpp
@@ -13,7 +13,7 @@
#include "moab/Range.hpp"
#include "moab/Error.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "MBTagConventions.hpp"
#include "EntitySequence.hpp"
#include "Internals.hpp"
diff --git a/src/io/ReadGCRM.cpp b/src/io/ReadGCRM.cpp
index 19d789e..0b86bd6 100644
--- a/src/io/ReadGCRM.cpp
+++ b/src/io/ReadGCRM.cpp
@@ -15,7 +15,7 @@
#include "moab/ReaderIface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "MBTagConventions.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#define ERRORR(rval, str) \
if (MB_SUCCESS != rval) {readMeshIface->report_error("%s", str); return rval;}
diff --git a/src/io/ReadHDF5.cpp b/src/io/ReadHDF5.cpp
index 916b69a..c707380 100644
--- a/src/io/ReadHDF5.cpp
+++ b/src/io/ReadHDF5.cpp
@@ -38,7 +38,7 @@
#include "MBTagConventions.hpp"
#include "ReadHDF5.hpp"
#include "moab/CN.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#ifdef HDF5_PARALLEL
# include <H5FDmpi.h>
# include <H5FDmpio.h>
diff --git a/src/io/ReadMCNP5.cpp b/src/io/ReadMCNP5.cpp
index c05ed41..b290f75 100644
--- a/src/io/ReadMCNP5.cpp
+++ b/src/io/ReadMCNP5.cpp
@@ -18,7 +18,7 @@
#include "moab/ReadUtilIface.hpp"
#include "Internals.hpp" // for MB_START_ID
#include "moab/Range.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <iostream>
#include <sstream>
diff --git a/src/io/ReadNASTRAN.cpp b/src/io/ReadNASTRAN.cpp
index 6534417..47a8780 100644
--- a/src/io/ReadNASTRAN.cpp
+++ b/src/io/ReadNASTRAN.cpp
@@ -27,7 +27,7 @@
#include "moab/ReadUtilIface.hpp"
#include "Internals.hpp" // for MB_START_ID
#include "moab/Range.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "FileTokenizer.hpp"
#include "MBTagConventions.hpp"
#include "moab/CN.hpp"
diff --git a/src/io/ReadNC.cpp b/src/io/ReadNC.cpp
index 50225df..a27a89d 100644
--- a/src/io/ReadNC.cpp
+++ b/src/io/ReadNC.cpp
@@ -3,7 +3,7 @@
#include "moab/ReadUtilIface.hpp"
#include "MBTagConventions.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#define ERRORR(rval, str) \
if (MB_SUCCESS != rval) { readMeshIface->report_error("%s", str); return rval; }
diff --git a/src/io/ReadNCDF.cpp b/src/io/ReadNCDF.cpp
index c69b3b5..01e98aa 100644
--- a/src/io/ReadNCDF.cpp
+++ b/src/io/ReadNCDF.cpp
@@ -38,7 +38,7 @@
#include "Internals.hpp"
#include "moab/ReadUtilIface.hpp"
#include "exodus_order.h"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/CartVect.hpp"
diff --git a/src/io/ReadSTL.cpp b/src/io/ReadSTL.cpp
index 4ba085d..59bd80b 100644
--- a/src/io/ReadSTL.cpp
+++ b/src/io/ReadSTL.cpp
@@ -25,7 +25,7 @@
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "SysUtil.hpp"
#include <errno.h>
diff --git a/src/io/ReadSmf.cpp b/src/io/ReadSmf.cpp
index 0970d71..eac232e 100644
--- a/src/io/ReadSmf.cpp
+++ b/src/io/ReadSmf.cpp
@@ -28,7 +28,7 @@
#include "Internals.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "AffineXform.hpp"
static inline int streq(const char *a,const char *b) { return strcmp(a,b)==0; }
diff --git a/src/io/ReadTemplate.cpp b/src/io/ReadTemplate.cpp
index 309a71e..97f70ce 100644
--- a/src/io/ReadTemplate.cpp
+++ b/src/io/ReadTemplate.cpp
@@ -8,7 +8,7 @@
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <cstdio>
#include <assert.h>
diff --git a/src/io/ReadTetGen.cpp b/src/io/ReadTetGen.cpp
index 72197fb..1d6609a 100644
--- a/src/io/ReadTetGen.cpp
+++ b/src/io/ReadTetGen.cpp
@@ -2,7 +2,7 @@
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "MBTagConventions.hpp"
#include <iostream>
#include <fstream>
diff --git a/src/io/ReadVtk.cpp b/src/io/ReadVtk.cpp
index ac3b9e4..2936a85 100644
--- a/src/io/ReadVtk.cpp
+++ b/src/io/ReadVtk.cpp
@@ -28,7 +28,7 @@
#include "Internals.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "FileTokenizer.hpp"
#include "VtkUtil.hpp"
diff --git a/src/io/Tqdcfr.cpp b/src/io/Tqdcfr.cpp
index b930ab6..717f43a 100644
--- a/src/io/Tqdcfr.cpp
+++ b/src/io/Tqdcfr.cpp
@@ -16,7 +16,7 @@
#include "Tqdcfr.hpp"
#include "moab/Core.hpp"
#include "moab/Range.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include <iostream>
#include <string>
diff --git a/src/io/WriteDamsel.hpp b/src/io/WriteDamsel.hpp
index 8902f08..fff9e8b 100644
--- a/src/io/WriteDamsel.hpp
+++ b/src/io/WriteDamsel.hpp
@@ -26,7 +26,7 @@
#include "moab/Range.hpp"
#include "moab/WriterIface.hpp"
#include "RangeSeqIntersectIter.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "DamselUtil.hpp"
#include "damsel.h"
diff --git a/src/io/WriteGmsh.cpp b/src/io/WriteGmsh.cpp
index aac3d38..7d4e2f7 100644
--- a/src/io/WriteGmsh.cpp
+++ b/src/io/WriteGmsh.cpp
@@ -5,7 +5,7 @@
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/WriteUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "GmshUtil.hpp"
#include <fstream>
diff --git a/src/io/WriteHDF5.cpp b/src/io/WriteHDF5.cpp
index 839c455..337cb96 100644
--- a/src/io/WriteHDF5.cpp
+++ b/src/io/WriteHDF5.cpp
@@ -48,17 +48,12 @@
#include "Internals.hpp"
#include "MBTagConventions.hpp"
#include "moab/CN.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/Version.h"
+#include "moab/CpuTimer.hpp"
#include "IODebugTrack.hpp"
#include "mhdf.h"
-#ifdef USE_MPI
-#define RUNTIME MPI_Wtime()
-#else
-#define RUNTIME (clock()/(double)CLOCKS_PER_SEC)
-#endif
-
/* Access HDF5 file handle for debugging
#include <H5Fpublic.h>
struct file { uint32_t magic; hid_t handle; };
@@ -100,15 +95,6 @@ struct file { uint32_t magic; hid_t handle; };
namespace moab {
-class CpuTimer {
-private:
- double atBirth, atLast;
-public:
- CpuTimer() : atBirth(RUNTIME), atLast(atBirth) {}
- double since_birth() { return (atLast = RUNTIME) - atBirth; };
- double elapsed() { double tmp = atLast; return (atLast = RUNTIME) - tmp; }
-};
-
template <typename T> inline
void VALGRIND_MAKE_VEC_UNDEFINED( std::vector<T>& v ) {
VALGRIND_MAKE_MEM_UNDEFINED( &v[0], v.size() * sizeof(T) );
@@ -601,7 +587,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
topState.end(result);
CHK_MB_ERR_0(result);
- times[GATHER_TIME] = timer.elapsed();
+ times[GATHER_TIME] = timer.time_elapsed();
//if (nodeSet.range.size() == 0)
// return error(MB_ENTITY_NOT_FOUND);
@@ -650,7 +636,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
if (MB_SUCCESS != result)
return error(result);
- times[CREATE_TIME] = timer.elapsed();
+ times[CREATE_TIME] = timer.time_elapsed();
dbgOut.tprint(1,"Writing Nodes.\n");
// Write node coordinates
@@ -662,7 +648,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
return error(result);
}
- times[COORD_TIME] = timer.elapsed();
+ times[COORD_TIME] = timer.time_elapsed();
dbgOut.tprint(1,"Writing connectivity.\n");
@@ -674,7 +660,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
if (MB_SUCCESS != result)
return error(result);
}
- times[CONN_TIME] = timer.elapsed();
+ times[CONN_TIME] = timer.time_elapsed();
dbgOut.tprint(1,"Writing sets.\n");
@@ -684,7 +670,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
return error(result);
debug_barrier();
- times[SET_TIME] = timer.elapsed();
+ times[SET_TIME] = timer.time_elapsed();
dbgOut.tprint(1,"Writing adjacencies.\n");
// Write adjacencies
@@ -701,7 +687,7 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
if (MB_SUCCESS != result)
return error(result);
}
- times[ADJ_TIME] = timer.elapsed();
+ times[ADJ_TIME] = timer.time_elapsed();
dbgOut.tprint(1,"Writing tags.\n");
@@ -716,9 +702,9 @@ ErrorCode WriteHDF5::write_file_impl( const char* filename,
if (MB_SUCCESS != result)
return error(result);
}
- times[TAG_TIME] = timer.elapsed();
+ times[TAG_TIME] = timer.time_elapsed();
- times[TOTAL_TIME] = timer.since_birth();
+ times[TOTAL_TIME] = timer.time_since_birth();
if (cputime) {
print_times( times );
@@ -1374,7 +1360,7 @@ ErrorCode WriteHDF5::write_sets( double* times )
mhdf_closeData( filePtr, table, &status );
CHK_MHDF_ERR_0(status);
- times[SET_PARENT] = timer.elapsed();
+ times[SET_PARENT] = timer.time_elapsed();
track.all_reduce();
}
@@ -1393,7 +1379,7 @@ ErrorCode WriteHDF5::write_sets( double* times )
mhdf_closeData( filePtr, table, &status );
CHK_MHDF_ERR_0(status);
- times[SET_CHILD] = timer.elapsed();
+ times[SET_CHILD] = timer.time_elapsed();
track.all_reduce();
}
@@ -1415,7 +1401,7 @@ ErrorCode WriteHDF5::write_sets( double* times )
mhdf_closeData( filePtr, table, &status );
CHK_MHDF_ERR_0(status);
- times[SET_CONTENT] = timer.elapsed();
+ times[SET_CONTENT] = timer.time_elapsed();
track.all_reduce();
}
assert( ranged_sets.size() + null_stripped_sets.size() == set_sizes.size() );
@@ -1544,7 +1530,7 @@ ErrorCode WriteHDF5::write_sets( double* times )
mhdf_closeData( filePtr, table, &status );
CHK_MHDF_ERR_0(status);
- times[SET_META] = timer.elapsed();
+ times[SET_META] = timer.time_elapsed();
track_meta.all_reduce();
return MB_SUCCESS;
@@ -1881,7 +1867,7 @@ ErrorCode WriteHDF5::write_tag( const TagDesc& tag_data,
if (array_len == MB_VARIABLE_LENGTH && tag_data.write_sparse) {
dbgOut.printf( 2, "Writing sparse data for var-len tag: \"%s\"\n", name.c_str() );
rval = write_var_len_tag( tag_data, name, moab_type, hdf5_type, elem_size );
- times[VARLEN_TAG_TIME] += timer.elapsed();
+ times[VARLEN_TAG_TIME] += timer.time_elapsed();
}
else {
int data_len = elem_size;
@@ -1890,7 +1876,7 @@ ErrorCode WriteHDF5::write_tag( const TagDesc& tag_data,
if (tag_data.write_sparse) {
dbgOut.printf( 2, "Writing sparse data for tag: \"%s\"\n", name.c_str() );
rval = write_sparse_tag( tag_data, name, moab_type, hdf5_type, data_len );
- times[SPARSE_TAG_TIME] += timer.elapsed();
+ times[SPARSE_TAG_TIME] += timer.time_elapsed();
}
for (size_t i = 0; MB_SUCCESS == rval && i < tag_data.dense_list.size(); ++i) {
const ExportSet* set = find( tag_data.dense_list[i] );
@@ -1901,7 +1887,7 @@ ErrorCode WriteHDF5::write_tag( const TagDesc& tag_data,
rval = write_dense_tag( tag_data, *set, name, moab_type, hdf5_type, data_len );
subState.end(rval);
}
- times[DENSE_TAG_TIME] += timer.elapsed();
+ times[DENSE_TAG_TIME] += timer.time_elapsed();
}
H5Tclose( hdf5_type );
diff --git a/src/io/WriteSTL.cpp b/src/io/WriteSTL.cpp
index 62644a6..45718d8 100644
--- a/src/io/WriteSTL.cpp
+++ b/src/io/WriteSTL.cpp
@@ -25,7 +25,7 @@
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/WriteUtilIface.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "SysUtil.hpp"
#include <stdio.h>
diff --git a/src/io/WriteSmf.cpp b/src/io/WriteSmf.cpp
index 0a09b8b..9d7f5f0 100644
--- a/src/io/WriteSmf.cpp
+++ b/src/io/WriteSmf.cpp
@@ -40,7 +40,7 @@
#include "MBTagConventions.hpp"
#include "moab/WriteUtilIface.hpp"
#include "Internals.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/Version.h"
namespace moab {
diff --git a/src/io/WriteVtk.cpp b/src/io/WriteVtk.cpp
index 03e95e4..bfda520 100644
--- a/src/io/WriteVtk.cpp
+++ b/src/io/WriteVtk.cpp
@@ -41,7 +41,7 @@
#include "MBTagConventions.hpp"
#include "moab/WriteUtilIface.hpp"
#include "Internals.hpp"
-#include "FileOptions.hpp"
+#include "moab/FileOptions.hpp"
#include "moab/Version.h"
#define INS_ID(stringvar, prefix, id) \
diff --git a/src/lotte/errmem.c b/src/lotte/errmem.c
index 0196d83..b890bce 100644
--- a/src/lotte/errmem.c
+++ b/src/lotte/errmem.c
@@ -1,6 +1,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <stdarg.h>
+#include "moab/FindPtFuncs.h"
void fail(const char *fmt, ...)
{
diff --git a/src/lotte/findpt.c b/src/lotte/findpt.c
index 00fe8b7..160c7ab 100644
--- a/src/lotte/findpt.c
+++ b/src/lotte/findpt.c
@@ -5,12 +5,10 @@
#include <float.h>
#include <string.h> /* for memcpy */
-#include "errmem.h"
-#include "types.h"
-#include "minmax.h"
-#include "poly.h"
-#include "tensor.h"
-#include "extrafindpt.h"
+#include "moab/FindPtFuncs.h"
+
+const unsigned opt_no_constraints_2 = 3+1;
+const unsigned opt_no_constraints_3 = 9+3+1;
/*--------------------------------------------------------------------------
Lobatto Polynomial Bounds
@@ -484,14 +482,6 @@ static void obbox_free_3(obbox_data_3 *p)
free(p);
}
-typedef struct {
- real x[2], A[4], axis_bnd[4];
-} obbox_2;
-
-typedef struct {
- real x[3], A[9], axis_bnd[6];
-} obbox_3;
-
int obbox_axis_test_2(const obbox_2 *p, const real x[2])
{
return (x[0]<p->axis_bnd[0] || x[0]>p->axis_bnd[1] ||
@@ -767,34 +757,6 @@ void obbox_calc_3(const obbox_data_3 *p, real tol,
--------------------------------------------------------------------------*/
-typedef struct {
- unsigned hash_n;
- real bnd[4]; /* bounds for all elements */
- real fac[2]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
- obbox_2 *obb; /* obb[nel] -- bounding box info for each element */
- uint *offset; /* hash table -- for cell i,j:
- uint index = j*hash_n+i,
- b = offset[index ],
- e = offset[index+1];
- elements in cell are
- offset[b], offset[b+1], ..., offset[e-1] */
- unsigned max; /* maximum # of elements in any cell */
-} hash_data_2;
-
-typedef struct {
- unsigned hash_n;
- real bnd[6]; /* bounds for all elements */
- real fac[3]; /* fac[i] = hash_n / (bnd[2*i+1]-bnd[2*i]) */
- obbox_3 *obb; /* obb[nel] -- bounding box info for each element */
- uint *offset; /* hash table -- for cell i,j,k:
- uint index = (k*hash_n+j)*hash_n+i,
- b = offset[index ],
- e = offset[index+1];
- elements in cell are
- offset[b], offset[b+1], ..., offset[e-1] */
- unsigned max; /* maximum # of elements in any cell */
-} hash_data_3;
-
static int ifloor(real x)
{
/*
@@ -1885,12 +1847,6 @@ double opt_findpt_2(opt_data_2 *p, const real *const elx[2],
--------------------------------------------------------------------------*/
-typedef struct {
- uint el;
- real r[3];
- real dist;
-} findpt_listel;
-
/* heap sort on A[0:n-1] with key A[i]->dist
precondition: n!=0 */
static void findpt_list_sort(findpt_listel **A, unsigned n)
@@ -1925,34 +1881,6 @@ static void findpt_list_sort(findpt_listel **A, unsigned n)
}
}
-typedef struct {
- const real *xw[2]; /* geometry data */
- real *z[2]; /* lobatto nodes */
- lagrange_data ld[2]; /* interpolation, derivative weights & data */
- unsigned nptel; /* nodes per element */
- hash_data_2 *hash; /* geometric hashing data */
- findpt_listel *list, **sorted, **end; /* pre-allocated list of elements to
- check (found by hashing), and
- pre-allocated list of pointers into
- the first list for sorting */
- opt_data_2 *od; /* data for the optimization algorithm */
- real *od_work;
-} findpt_data_2;
-
-typedef struct {
- const real *xw[3]; /* geometry data */
- real *z[3]; /* lobatto nodes */
- lagrange_data ld[3]; /* interpolation, derivative weights & data */
- unsigned nptel; /* nodes per element */
- hash_data_3 *hash; /* geometric hashing data */
- findpt_listel *list, **sorted, **end; /* pre-allocated list of elements to
- check (found by hashing), and
- pre-allocated list of pointers into
- the first list for sorting */
- opt_data_3 *od; /* data for the optimization algorithm */
- real *od_work;
-} findpt_data_3;
-
findpt_data_2 *findpt_setup_2(
const real *const xw[2], const unsigned n[2], uint nel,
uint max_hash_size, real bbox_tol)
diff --git a/src/lotte/poly.c b/src/lotte/poly.c
index 6b5cba6..b8a93b7 100644
--- a/src/lotte/poly.c
+++ b/src/lotte/poly.c
@@ -5,8 +5,7 @@
#include <string.h> /* for memcpy */
#include <float.h>
-#include "errmem.h"
-#include "types.h"
+#include "moab/FindPtFuncs.h"
/*
For brevity's sake, some names have been shortened
@@ -354,15 +353,6 @@ void lagrange_weights_deriv(const real *z, unsigned n,
lagrange_free(&p); // deallocate memory allocated by setup
--------------------------------------------------------------------------*/
-typedef struct {
- unsigned n; /* number of Lagrange nodes */
- const real *z; /* Lagrange nodes (user-supplied) */
- real *J, *D, *D2; /* weights for 0th,1st,2nd derivatives */
- real *J_z0, *D_z0, *D2_z0; /* ditto at z[0] (computed at setup) */
- real *J_zn, *D_zn, *D2_zn; /* ditto at z[n-1] (computed at setup) */
- real *w, *d, *u0, *v0, *u1, *v1, *u2, *v2; /* work data */
-} lagrange_data;
-
void lagrange_0(lagrange_data *p, real x)
{
unsigned i, n=p->n;
diff --git a/src/lotte/tensor.c b/src/lotte/tensor.c
index 19e54b7..510dc38 100644
--- a/src/lotte/tensor.c
+++ b/src/lotte/tensor.c
@@ -1,4 +1,4 @@
-#include "types.h"
+#include "moab/FindPtFuncs.h"
/*--------------------------------------------------------------------------
Matrix-Matrix Multiplication
diff --git a/src/moab/CartVect.hpp b/src/moab/CartVect.hpp
index 383f8e4..a249fef 100644
--- a/src/moab/CartVect.hpp
+++ b/src/moab/CartVect.hpp
@@ -1,18 +1,3 @@
-/**
- * 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.
- *
- */
-
#ifndef MB_CART_VECT_HPP
#define MB_CART_VECT_HPP
@@ -24,8 +9,6 @@ namespace moab {
/**
* \brief Cartesian Vector
- * \author Jason Kraftcheck
- * \date July, 2006
*/
class CartVect
{
@@ -67,6 +50,8 @@ class CartVect
{ d[0] *= s; d[1] *= s; d[2] *= s; return *this; }
inline CartVect& operator/=( double s )
{ d[0] /= s; d[1] /= s; d[2] /= s; return *this; }
+ inline bool operator==(const CartVect& v ) const
+ { return d[0] == v[0] && d[1] == v[1] && d[2] == v[2]; }
inline double length() const; //!< vector length
diff --git a/src/moab/Core.hpp b/src/moab/Core.hpp
index efa17e3..4bb86c8 100644
--- a/src/moab/Core.hpp
+++ b/src/moab/Core.hpp
@@ -221,7 +221,7 @@ public:
\param entity_handle EntityHandle to get connectivity of.
\param connectivity Vector in which connectivity of <em>entity_handle</em> is returned.
Should contain MeshVertices.
- \param topological_connectivity If true, higher order nodes are ignored.
+ \param corners_only If true, higher order nodes are ignored.
Example: \code
std::vector<EntityHandle> conn;
@@ -229,7 +229,7 @@ public:
virtual ErrorCode get_connectivity(const EntityHandle *entity_handles,
const int num_handles,
std::vector<EntityHandle> &connectivity,
- bool topological_connectivity = false,
+ bool corners_only = false,
std::vector<int> *offsets = NULL) const;
//! Gets the connectivity for a vector of elements
@@ -238,14 +238,14 @@ public:
virtual ErrorCode get_connectivity(const EntityHandle *entity_handles,
const int num_handles,
Range &connectivity,
- bool topological_connectivity = false) const;
+ bool corners_only = false) const;
//! Gets the connectivity for elements
/** Same as vector-based version except range is returned (unordered!)
*/
virtual ErrorCode get_connectivity( const Range& entity_handles,
Range &connectivity,
- bool topological_connectivity = false) const;
+ bool corners_only = false) const;
//! Gets a pointer to constant connectivity data of <em>entity_handle</em>
/** Sets <em>number_nodes</em> equal to the number of nodes of the <em>
@@ -273,7 +273,7 @@ public:
virtual ErrorCode get_connectivity( const EntityHandle entity_handle,
const EntityHandle *&connectivity,
int &num_nodes,
- bool topological_connectivity = false,
+ bool corners_only = false,
std::vector<EntityHandle>* storage = 0
) const;
@@ -1134,6 +1134,8 @@ public:
void print(const EntityHandle handle, const char *prefix,
bool first_call = true) const;
+ ErrorCode print_entity_tags(std::string indent_prefix, const EntityHandle handle, TagType tp) const;
+
virtual ErrorCode get_last_error(std::string& info) const;
virtual std::string get_error_string(const ErrorCode code) const;
diff --git a/src/moab/CpuTimer.hpp b/src/moab/CpuTimer.hpp
new file mode 100644
index 0000000..cdda903
--- /dev/null
+++ b/src/moab/CpuTimer.hpp
@@ -0,0 +1,63 @@
+#ifndef CPUTIMER_HPP
+#define CPUTIMER_HPP
+
+#ifdef USE_MPI
+# include "moab_mpi.h"
+#else
+# include <sys/resource.h>
+#endif
+
+namespace moab
+{
+
+class CpuTimer {
+private:
+ double tAtBirth, tAtLast;
+ double mAtBirth, mAtLast;
+ long rssAtBirth, rssAtLast;
+
+ double runtime();
+ long runmem();
+
+public:
+ CpuTimer() : tAtBirth(runtime()), tAtLast(tAtBirth) {}
+ double time_since_birth() { return (tAtLast = runtime()) - tAtBirth; };
+ double time_elapsed() { double tmp = tAtLast; return (tAtLast = runtime()) - tmp; }
+ long mem_since_birth() {return (mAtLast=runmem()) - mAtBirth;}
+ long mem_elapsed() {long tmp = mAtLast; return (mAtLast=runmem()) - tmp;}
+};
+
+ inline double CpuTimer::runtime()
+ {
+#if defined(_MSC_VER) || defined(__MINGW32__)
+ return (double)clock() / CLOCKS_PER_SEC;
+#elif defined(USE_MPI)
+ return MPI_Wtime();
+#else
+ struct rusage r_usage;
+ getrusage(RUSAGE_SELF, &r_usage);
+ double utime = (double)r_usage.ru_utime.tv_sec +
+ ((double)r_usage.ru_utime.tv_usec/1.e6);
+ double stime = (double)r_usage.ru_stime.tv_sec +
+ ((double)r_usage.ru_stime.tv_usec/1.e6);
+ return utime + stime;
+#endif
+ }
+
+ inline long CpuTimer::runmem()
+ {
+#if defined(_MSC_VER) || defined(__MINGW32__)
+ return 0;
+#elif defined(USE_MPI)
+ return 0;
+#else
+ struct rusage r_usage;
+ getrusage(RUSAGE_SELF, &r_usage);
+ mAtLast = r_usage.ru_maxrss;
+ return mAtLast;
+#endif
+ }
+
+}
+
+#endif
diff --git a/src/moab/FileOptions.hpp b/src/moab/FileOptions.hpp
new file mode 100644
index 0000000..85b3eef
--- /dev/null
+++ b/src/moab/FileOptions.hpp
@@ -0,0 +1,248 @@
+/*
+ * 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.
+ *
+ */
+
+/**\file FileOptions.hpp
+ *\author Jason Kraftcheck (kraftche at cae.wisc.edu)
+ *\date 2007-08-21
+ */
+
+#ifndef FILE_OPTIONS_HPP
+#define FILE_OPTIONS_HPP
+
+#include <string>
+#include <vector>
+#include "moab/Types.hpp"
+
+namespace moab {
+
+/**\brief Parse options string passed to file IO routines
+ *
+ * This is a utility class used by file-IO-related code to
+ * parse the options string passed to Core::load_file and
+ * Core::write_file
+ */
+class FileOptions {
+public:
+
+ /*\param options_string The concatenation of a list of
+ * options, separated either by the default separator
+ * (semicolon) with a custom separator specified at
+ * the beginning of the string (semicolon followed by
+ * destired separator character.)
+ */
+ FileOptions( const char* option_string );
+
+ FileOptions( const FileOptions& copy );
+ FileOptions& operator=( const FileOptions& copy );
+
+ ~FileOptions();
+
+ /**\brief Check for option with no value
+ *
+ * Check for an option w/out a value.
+ *\param name The option name
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but has value
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_null_option( const char* name ) const;
+
+
+ /**\brief Check for option with boolean (true/false, yes/no) value.
+ *
+ * Check for an option with a true/false value. Allowable values
+ * are "true", "false", "yes", "no", "1", "0", "on", "off".
+ *\param name The option name
+ *\param default_value The value to return if the option is not specified.
+ *\param value The resulting value. This argument is set to the passed
+ * default value if the option is not found.
+ *\return - MB_TYPE_OUT_OF_RANGE if options is found, but has an invalid value
+ * - MB_SUCCESS otherwise
+ */
+ ErrorCode get_toggle_option( const char* name,
+ bool default_value,
+ bool& value ) const;
+
+ /**\brief Check for option with an integer value.
+ *
+ * Check for an option with an integer value
+ *\param name The option name
+ *\param value Output. The value.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have an integer value
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_int_option( const char* name, int& value ) const;
+
+ /**\brief Check for option with an integer value. Accept option with no value.
+ *
+ * Check for an option with an integer value.
+ * If the option is found but has no value specified, then
+ * pass back the user-specified default value.
+ *
+ *\NOTE: This function will not pass back the default_val, but will instead
+ * return MB_ENTITY_NOT_FOUND if the option is not specified at all.
+ * The default value is returned only when the option is specified,
+ * but is specified w/out a value.
+ *
+ *\param name The option name
+ *\param default_val The default value for the option.
+ *\param value Output. The value.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found but has a value that cannot be parsed as an int
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_int_option( const char* name, int default_val, int& value ) const;
+
+ /**\brief Check for option with a double value.
+ *
+ * Check for an option with a double value
+ *\param name The option name
+ *\param value Output. The value.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have a double value
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_real_option( const char* name, double& value ) const;
+
+ /**\brief Check for option with any value.
+ *
+ * Check for an option with any value.
+ *\param name The option name
+ *\param value Output. The value.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not have a value
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_str_option( const char* name, std::string& value ) const;
+
+ /**\brief Check for option
+ *
+ * Check for an option
+ *\param name The option name
+ *\param value The option value, or an empty string if no value.
+ *\return MB_SUCCESS or MB_ENTITY_NOT_FOUND
+ */
+ ErrorCode get_option( const char* name, std::string& value ) const;
+
+ /**\brief Check the string value of an option
+ *
+ * Check which of a list of possible values a string option contains.
+ *\param name The option name
+ *\param values A NULL-terminated array of C-style strings enumerating
+ * the possible option values.
+ *\param index Output: The index into <code>values</code> for the
+ * option value.
+ *\return MB_SUCCESS if matched name and value.
+ * MB_ENTITY_NOT_FOUND if the option was not specified
+ * MB_FAILURE if the option value is not in the input <code>values</code> array.
+ */
+ ErrorCode match_option( const char* name, const char* const* values, int& index ) const;
+
+ /**\brief Check if an option matches a string value
+ *
+ * Check if the value for an option is the passed string.
+ *\param name The option name
+ *\param value The expected value.
+ *\return MB_SUCCESS if matched name and value.
+ * MB_ENTITY_NOT_FOUND if the option was not specified
+ * MB_FAILURE if the option value doesn't match the passed string/
+ */
+ ErrorCode match_option( const char* name, const char* value) const;
+
+ /**\brief Check for option for which the value is a list of ints
+ *
+ * Check for an option which is an int list. The value is expected to
+ * be a comma-separated list of int ranges, where an int range can be
+ * either a single integer value or a range of integer values separated
+ * by a dash ('-').
+ *
+ *\param name The option name
+ *\param values Output. The list of integer values.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain an ID list
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_ints_option( const char* name, std::vector<int>& values) const;
+
+ /**\brief Check for option for which the value is a list of doubles
+ *
+ * Check for an option which is a double list. The value is expected to
+ * be a comma-separated list of double values
+ *
+ *\param name The option name
+ *\param values Output. The list of double values.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain an ID list
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_reals_option( const char* name, std::vector<double>& values) const;
+
+ /**\brief Check for option for which the value is a list of strings
+ *
+ * Check for an option which is a string list. The value is expected to
+ * be a comma-separated list of string values, with no embedded spaces or commas.
+ *
+ *\param name The option name
+ *\param values Output. The list of string values.
+ *\return - MB_SUCCESS if option is found
+ * - MB_TYPE_OUT_OF_RANGE if options is found, but does not contain a string list
+ * - MB_ENTITY_NOT_FOUND if option is not found.
+ */
+ ErrorCode get_strs_option( const char* name, std::vector<std::string>& values) const;
+
+ /** number of options */
+ inline unsigned size() const
+ { return mOptions.size(); }
+
+ /** true if no options */
+ inline bool empty() const
+ { return mOptions.empty(); }
+
+ /** Get list of options */
+ void get_options( std::vector<std::string>& list ) const;
+
+ /** Check if all options have been looked at */
+ bool all_seen() const;
+
+ /** Mark all options as seen. USed for non-root procs during bcast-delete read */
+ void mark_all_seen() const;
+
+ /** Get first unseen option */
+ ErrorCode get_unseen_option( std::string& value ) const;
+
+private:
+
+ /**\brief Check for option
+ *
+ * Check for an option
+ *\param name The option name
+ *\param value The option value, or an empty string if no value.
+ *\return MB_SUCCESS or MB_ENTITY_NOT_FOUND
+ */
+ ErrorCode get_option( const char* name, const char*& value) const;
+
+ char* mData;
+ std::vector<const char*> mOptions;
+ mutable std::vector<bool> mSeen;
+
+ /** Case-insensitive compare of name with option value. */
+ static bool compare( const char* name, const char* option );
+};
+
+} // namespace moab
+
+#endif
+
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/a91108ec4055/
Changeset: a91108ec4055
Branch: None
User: tautges
Date: 2013-09-12 23:37:24
Summary: Merged fathomteam/moab into master
Affected #: 2 files
diff --git a/src/BoundBox.cpp b/src/BoundBox.cpp
new file mode 100644
index 0000000..cb34330
--- /dev/null
+++ b/src/BoundBox.cpp
@@ -0,0 +1,78 @@
+#include "moab/Interface.hpp"
+#include "moab/BoundBox.hpp"
+
+namespace moab
+{
+ ErrorCode BoundBox::update(Interface &iface, const Range& elems)
+ {
+ ErrorCode rval;
+ bMin = CartVect(HUGE_VAL);
+ bMax = CartVect(-HUGE_VAL);
+
+ CartVect coords;
+ EntityHandle const *conn, *conn2;
+ int len, len2;
+ Range::const_iterator i;
+
+ // vertices
+ const Range::const_iterator elem_begin = elems.lower_bound( MBEDGE );
+ for (i = elems.begin(); i != elem_begin; ++i) {
+ rval = iface.get_coords( &*i, 1, coords.array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+ update_min(coords.array());
+ update_max(coords.array());
+ }
+
+ // elements with vertex-handle connectivity list
+ const Range::const_iterator poly_begin = elems.lower_bound( MBPOLYHEDRON, elem_begin );
+ std::vector<EntityHandle> dum_vector;
+ for (i = elem_begin; i != poly_begin; ++i) {
+ rval = iface.get_connectivity( *i, conn, len, true, &dum_vector);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ for (int j = 0; j < len; ++j) {
+ rval = iface.get_coords( conn+j, 1, coords.array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+ update_min(coords.array());
+ update_max(coords.array());
+ }
+ }
+
+ // polyhedra
+ const Range::const_iterator set_begin = elems.lower_bound( MBENTITYSET, poly_begin );
+ for (i = poly_begin; i != set_begin; ++i) {
+ rval = iface.get_connectivity( *i, conn, len, true );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ for (int j = 0; j < len; ++j) {
+ rval = iface.get_connectivity( conn[j], conn2, len2 );
+ for (int k = 0; k < len2; ++k) {
+ rval = iface.get_coords( conn2+k, 1, coords.array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+ update_min(coords.array());
+ update_max(coords.array());
+ }
+ }
+ }
+
+ // sets
+ BoundBox box;
+ for (i = set_begin; i != elems.end(); ++i) {
+ Range tmp_elems;
+ rval = iface.get_entities_by_handle(*i, tmp_elems);
+ if (MB_SUCCESS != rval) return rval;
+ rval = box.update(iface, tmp_elems);
+ if (MB_SUCCESS != rval) return rval;
+
+ update(box);
+ }
+
+ return MB_SUCCESS;
+ }
+
+}
diff --git a/src/moab/BoundBox.hpp b/src/moab/BoundBox.hpp
new file mode 100644
index 0000000..9400944
--- /dev/null
+++ b/src/moab/BoundBox.hpp
@@ -0,0 +1,177 @@
+#ifndef BOUND_BOX_HPP
+#define BOUND_BOX_HPP
+
+#include "moab/Interface.hpp"
+#include "moab/CartVect.hpp"
+
+#include <cfloat>
+
+namespace moab {
+
+ class BoundBox {
+ public:
+ BoundBox() : bMin(DBL_MAX), bMax(-DBL_MAX) {}
+ BoundBox(const CartVect &min, const CartVect &max) :
+ bMin(min), bMax(max) {}
+ BoundBox(const double *corners);
+ ~BoundBox() {}
+
+ bool contains_point(const double *point, const double tol = 0.0) const;
+ bool contains_box(const BoundBox &b, const double tol = 0.0) const;
+ void compute_center(CartVect ¢er);
+ void update(const BoundBox &other_box);
+ void update(const double *coords);
+ ErrorCode update(Interface &iface, const Range& elems);
+ ErrorCode update(Interface &iface, const EntityHandle ent);
+ void update_min(const BoundBox &other_box);
+ void update_min(const double *coords);
+ void update_max(const BoundBox &other_box);
+ void update_max(const double *coords);
+
+ /** \brief Return the diagonal length of this box
+ */
+ double diagonal_length() const;
+
+ /** \brief Return the square of the diagonal length of this box
+ */
+ double diagonal_squared() const;
+
+ /** \brief Return square of distance from box, or zero if inside
+ * \param from_point Point from which you want distance_sq
+ */
+ double distance_squared(const double *from_point) const;
+
+ /** \brief Return distance from box, or zero if inside
+ * \param from_point Point from which you want distance
+ */
+ double distance(const double *from_point) const;
+
+ BoundBox &operator=(const BoundBox &from) {
+ bMin = from.bMin;
+ bMax = from.bMax;
+ return *this;
+ }
+ inline bool operator==(const BoundBox &box) const {
+ return (bMin == box.bMin && bMax == box.bMax);
+ }
+
+ CartVect bMin, bMax;
+ };
+
+ inline BoundBox::BoundBox(const double *corners)
+ {
+ // relies on CartVect being Plain Old Data, no virtual table
+ double *arr = bMin.array();
+ for (int i = 0; i < 6; i++)
+ arr[i] = corners[i];
+ }
+
+ inline bool BoundBox::contains_point(const double *point, const double tol) const {
+ if (point[0] < bMin[0]-tol || point[0] > bMax[0]+tol ||
+ point[1] < bMin[1]-tol || point[1] > bMax[1]+tol ||
+ point[2] < bMin[2]-tol || point[2] > bMax[2]+tol)
+ return false;
+ else return true;
+ }
+
+ inline bool BoundBox::contains_box(const BoundBox &b, const double tol) const {
+ if (b.bMin[0] < bMin[0]-tol || b.bMax[0] > bMax[0]+tol ||
+ b.bMin[1] < bMin[1]-tol || b.bMax[1] > bMax[1]+tol ||
+ b.bMin[2] < bMin[2]-tol || b.bMax[2] > bMax[2]+tol)
+ return false;
+
+ else return true;
+ }
+
+ inline void BoundBox::update(const BoundBox &other_box)
+ {
+ update_min(other_box);
+ update_max(other_box);
+ }
+
+ inline void BoundBox::update(const double *coords)
+ {
+ update_min(coords);
+ update_max(coords+3);
+ }
+
+ inline void BoundBox::update_min(const BoundBox &other_box)
+ {
+ bMin[0] = std::min(bMin[0], other_box.bMin[0]);
+ bMin[1] = std::min(bMin[1], other_box.bMin[1]);
+ bMin[2] = std::min(bMin[2], other_box.bMin[2]);
+ }
+
+ inline void BoundBox::update_min(const double *coords)
+ {
+ bMin[0] = std::min(bMin[0], coords[0]);
+ bMin[1] = std::min(bMin[1], coords[1]);
+ bMin[2] = std::min(bMin[2], coords[2]);
+ }
+
+ inline void BoundBox::update_max(const BoundBox &other_box)
+ {
+ bMax[0] = std::max(bMax[0], other_box.bMax[0]);
+ bMax[1] = std::max(bMax[1], other_box.bMax[1]);
+ bMax[2] = std::max(bMax[2], other_box.bMax[2]);
+ }
+
+ inline void BoundBox::update_max(const double *coords)
+ {
+ bMax[0] = std::max(bMax[0], coords[0]);
+ bMax[1] = std::max(bMax[1], coords[1]);
+ bMax[2] = std::max(bMax[2], coords[2]);
+ }
+
+ inline void BoundBox::compute_center(CartVect ¢er){
+ center = 0.5 * (bMin + bMax);
+ }
+
+ inline std::ostream &operator<<(std::ostream& out, const BoundBox &box) {
+ out << "Min: ";
+ out << box.bMin;
+ out << ", Max: ";
+ out << box.bMax;
+ return out;
+ }
+
+ inline ErrorCode BoundBox::update(Interface &iface, const EntityHandle ent)
+ {
+ Range tmp_range(ent, ent);
+ return update(iface, tmp_range);
+ }
+
+ inline double BoundBox::distance_squared(const double *from_point) const
+ {
+ double dist_sq = 0.0;
+ for (int i = 0; i < 3; ++i) {
+ if (from_point[i] < bMin[i])
+ dist_sq += (bMin[i] - from_point[i]) * (bMin[i] - from_point[i]);
+ else if (from_point[i] > bMax[i])
+ dist_sq += (bMax[i] - from_point[i]) * (bMax[i] - from_point[i]);
+ }
+ return dist_sq;
+ }
+
+ inline double BoundBox::distance(const double *from_point) const
+ {
+ double dist_sq = distance_squared(from_point);
+ return sqrt(dist_sq);
+ }
+
+ inline double BoundBox::diagonal_length() const
+ {
+ if (DBL_MAX == bMax[0] || DBL_MAX == bMax[1] || DBL_MAX == bMax[2] ||
+ DBL_MAX == bMin[0] || DBL_MAX == bMin[1] || DBL_MAX == bMin[2]) return DBL_MAX;
+ return (bMax - bMin).length();
+ }
+
+ inline double BoundBox::diagonal_squared() const
+ {
+ if (DBL_MAX == bMax[0] || DBL_MAX == bMax[1] || DBL_MAX == bMax[2] ||
+ DBL_MAX == bMin[0] || DBL_MAX == bMin[1] || DBL_MAX == bMin[2]) return DBL_MAX;
+ return (bMax - bMin).length_squared();
+ }
+}
+
+#endif
https://bitbucket.org/fathomteam/moab/commits/00c578e5a49f/
Changeset: 00c578e5a49f
Branch: None
User: tautges
Date: 2013-09-12 23:53:37
Summary: Adding a Lloyd smoother tool, and a test of it.
Affected #: 6 files
diff --git a/.gitignore b/.gitignore
index d3721c3..5aac80a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -155,6 +155,7 @@ test/io/vtk_test
test/kd_tree_test
test/kd_tree_time
test/kd_tree_tool
+test/lloyd_smoother_test
test/mbcn_test
test/mbfacet_test
test/mbground_test
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
new file mode 100644
index 0000000..49177d0
--- /dev/null
+++ b/src/LloydSmoother.cpp
@@ -0,0 +1,220 @@
+#include "moab/LloydSmoother.hpp"
+#include "moab/Error.hpp"
+#include "moab/Skinner.hpp"
+#include "moab/CN.hpp"
+#include "moab/CartVect.hpp"
+#include "moab/BoundBox.hpp"
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
+#include <iostream>
+
+#define RR(a) if (MB_SUCCESS != rval) { \
+ errorHandler->set_last_error(a); \
+ return rval;}
+
+namespace moab {
+
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
+ double at, double rt)
+ : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+ reportIts(0), numIts(0), iCreatedTag(false)
+{
+ ErrorCode rval = mbImpl->query_interface(errorHandler);
+ if (rval) throw rval;
+}
+
+LloydSmoother::~LloydSmoother()
+{
+ if (iCreatedTag && fixedTag) {
+ ErrorCode rval = mbImpl->tag_delete(fixedTag);
+ if (rval) throw rval;
+ }
+}
+
+ErrorCode LloydSmoother::perform_smooth()
+{
+ ErrorCode rval;
+
+ // first figure out tolerance to use
+ if (0 > absTol) {
+ // no tolerance set - get one relative to bounding box around elements
+ BoundBox bb;
+ rval = bb.update(*mbImpl, myElems);
+ RR("Failed to compute bounding box around elements.");
+ absTol = relTol * bb.diagonal_length();
+ }
+
+ // initialize if we need to
+ rval = initialize();
+ RR("Failed to initialize.");
+
+ // get all vertices
+ Range verts;
+ rval = mbImpl->get_adjacencies(myElems, 0, false, verts, Interface::UNION);
+ RR("Failed to get all vertices.");
+
+ // perform Lloyd relaxation:
+ // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
+
+ // get all verts coords into tag; don't need to worry about filtering out fixed verts,
+ // we'll just be setting to their fixed coords
+ std::vector<double> vcentroids(3*verts.size());
+ rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+
+ Tag centroid;
+ rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE);
+ RR("Couldn't get tag handle.");
+ rval = mbImpl->tag_set_data(centroid, verts, &vcentroids[0]); RR("Failed setting centroid tag.");
+
+ Range owned_verts, shared_owned_verts;
+#ifdef USE_MPI
+ // filter verts down to owned ones and get fixed tag for them
+ if (myPcomm->size() > 1) {
+ rval = myPcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
+ RR("Failed to filter on pstatus.");
+ // get shared owned verts, for exchanging tags
+ rval = myPcomm->filter_pstatus(owned_verts, PSTATUS_SHARED, PSTATUS_AND, -1, &shared_owned_verts);
+ RR("Failed to filter for shared owned.");
+ // workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
+ // for all shared entities
+ if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+ }
+ else
+ owned_verts = verts;
+#else
+ owned_verts = verts;
+#endif
+
+ std::vector<unsigned char> fix_tag(owned_verts.size());
+ rval = mbImpl->tag_get_data(fixedTag, owned_verts, &fix_tag[0]); RR("Failed to get fixed tag.");
+
+ // now fill vcentroids array with positions of just owned vertices, since those are the ones
+ // we're actually computing
+ vcentroids.resize(3*owned_verts.size());
+ rval = mbImpl->tag_get_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to get centroid tag.");
+
+ // some declarations for later iterations
+ std::vector<double> fcentroids(3*myElems.size()); // fcentroids for element centroids
+ std::vector<double> ctag(3*CN::MAX_NODES_PER_ELEMENT); // temporary coordinate storage for verts bounding an element
+ const EntityHandle *conn; // const ptr & size to elem connectivity
+ int nconn;
+ Range::iterator eit, vit; // for iterating over elems, verts
+ int e, v; // for indexing into centroid vectors
+ std::vector<EntityHandle> adj_elems; // used in vertex iteration
+
+ // 2. while !converged
+ double resid = DBL_MAX;
+ numIts = 0;
+ while (resid > absTol) {
+ numIts++;
+ resid = 0.0;
+
+ // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
+ for (eit = myElems.begin(), e = 0; eit != myElems.end(); eit++, e++) {
+ // get verts for this elem
+ rval = mbImpl->get_connectivity(*eit, conn, nconn); RR("Failed to get connectivity.");
+ // get centroid tags for those verts
+ rval = mbImpl->tag_get_data(centroid, conn, nconn, &ctag[0]); RR("Failed to get centroid.");
+ fcentroids[3*e+0] = fcentroids[3*e+1] = fcentroids[3*e+2] = 0.0;
+ for (v = 0; v < nconn; v++) {
+ fcentroids[3*e+0] += ctag[3*v+0];
+ fcentroids[3*e+1] += ctag[3*v+1];
+ fcentroids[3*e+2] += ctag[3*v+2];
+ }
+ for (v = 0; v < 3; v++) fcentroids[3*e+v] /= nconn;
+ }
+ rval = mbImpl->tag_set_data(centroid, myElems, &fcentroids[0]); RR("Failed to set elem centroid.");
+
+ // 2b. foreach owned vertex:
+ for (vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); vit++, v++) {
+ // if !fixed
+ if (fix_tag[v]) continue;
+ // vertex centroid = sum(cell centroids)/ncells
+ adj_elems.clear();
+ rval = mbImpl->get_adjacencies(&(*vit), 1, 2, false, adj_elems); RR("Failed getting adjs.");
+ rval = mbImpl->tag_get_data(centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0]);
+ RR("Failed to get elem centroid.");
+ double vnew[] = {0.0, 0.0, 0.0};
+ for (e = 0; e < (int)adj_elems.size(); e++) {
+ vnew[0] += fcentroids[3*e+0];
+ vnew[1] += fcentroids[3*e+1];
+ vnew[2] += fcentroids[3*e+2];
+ }
+ for (e = 0; e < 3; e++) vnew[e] /= adj_elems.size();
+ double delta = (CartVect(vnew)-CartVect(&vcentroids[3*v])).length();
+ resid = (v ? std::max(resid, delta) : delta);
+ for (e = 0; e < 3; e++) vcentroids[3*v+e] = vnew[e];
+ }
+
+ // set the centroid tag; having them only in vcentroids array isn't enough, as vertex centroids are
+ // accessed randomly in loop over faces
+ rval = mbImpl->tag_set_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to set vertex centroid.");
+
+#ifdef USE_MPI
+ // 2c. exchange tags on owned verts
+ if (myPcomm->size() > 1) {
+ rval = pcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
+ }
+#endif
+
+ if (reportIts && !(numIts%reportIts)) {
+ double global_max = resid;
+#ifdef USE_MPI
+ // global reduce for maximum delta, then report it
+ if (myPcomm->size() > 1)
+ MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
+ if (!pcomm->rank())
+#endif
+ std::cout << "Max residual = " << global_max << std::endl;
+ }
+
+ } // end while
+
+ // write the tag back onto vertex coordinates
+ rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode LloydSmoother::initialize()
+{
+ ErrorCode rval = MB_SUCCESS;
+
+ // first, check for tag; if we don't have it, make one and mark skin as fixed
+ if (!fixedTag) {
+ unsigned char fixed = 0x0;
+ rval = mbImpl->tag_get_handle("", 1, MB_TYPE_OPAQUE, fixedTag, MB_TAG_DENSE|MB_TAG_CREAT, &fixed);
+ RR("Trouble making fixed tag.");
+ iCreatedTag = true;
+
+ // get the skin; get facets, because we might need to filter on shared entities
+ Skinner skinner(mbImpl);
+ Range skin, skin_verts;
+ rval = skinner.find_skin(0, myElems, false, skin);
+ RR("Unable to find skin.");
+
+#ifdef USE_MPI
+ // need to do a little extra if we're working in parallel
+ if (myPcomm) {
+ // filter out ghost and interface facets
+ rval = myPcomm->filter_pstatus(skin, PSTATUS_GHOST|PSTATUS_INTERFACE, PSTATUS_NOT);
+ RR("Failed to filter on shared status.");
+ }
+#endif
+ // get the vertices from those entities
+ rval = mbImpl->get_adjacencies(skin, 0, false, skin_verts, Interface::UNION);
+ RR("Trouble getting vertices.");
+
+ // mark them fixed
+ std::vector<unsigned char> marks(skin.size(), 0x1);
+ rval = mbImpl->tag_set_data(fixedTag, skin_verts, &marks[0]);
+ RR("Unable to set tag on skin.");
+ }
+
+ return MB_SUCCESS;
+}
+
+}
diff --git a/src/Makefile.am b/src/Makefile.am
index ab5077a..3fb67b1 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -63,6 +63,7 @@ libMOAB_la_SOURCES = \
HigherOrderFactory.cpp \
HomXform.cpp \
Internals.hpp \
+ LloydSmoother.cpp \
MBCNArrays.hpp \
MergeMesh.cpp \
MeshSet.cpp \
@@ -161,6 +162,7 @@ nobase_libMOAB_la_include_HEADERS = \
moab/GeomUtil.hpp \
moab/Interface.hpp \
moab/EigenDecomp.hpp \
+ moab/LloydSmoother.hpp \
moab/point_locater/tree/common_tree.hpp \
moab/point_locater/tree/element_tree.hpp \
moab/point_locater/tree/bvh_tree.hpp \
diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
new file mode 100644
index 0000000..5ab9771
--- /dev/null
+++ b/src/moab/LloydSmoother.hpp
@@ -0,0 +1,145 @@
+/** @class LloydSmoother.cpp \n
+ * \brief Perform Lloyd relaxation on a mesh and its dual \n
+ *
+ * Briefly, Lloyd relaxation is a technique to smooth out a mesh. The centroid of each cell is computed from its
+ * vertex positions, then vertices are placed at the average of their connected cells' centroids, and the process
+ * iterates until convergence.
+ *
+ * In the parallel algorithm, an extra ghost layer of cells is exchanged. This allows us to compute the centroids
+ * for boundary cells on each processor where they appear; this eliminates the need for one round of data exchange
+ * (for those centroids) between processors. New vertex positions must be sent from owning processors to processors
+ * sharing those vertices. Convergence is measured as the maximum distance moved by any vertex.
+ *
+ */
+
+#ifndef LLOYDSMOOTHER_HPP
+#define LLOYDSMOOTHER_HPP
+
+#include "moab/Interface.hpp"
+
+namespace moab {
+
+class ParallelComm;
+class Error;
+
+class LloydSmoother
+{
+public:
+
+ /* \brief Constructor
+ * Bare constructor, data input to this class through methods.
+ * \param impl The MOAB instance for this smoother
+ */
+ LloydSmoother(Interface *impl);
+
+ /* \brief Constructor
+ * Convenience constructor, data input directly
+ * \param impl The MOAB instance for this smoother
+ * \param pcomm The ParallelComm instance by which this mesh is parallel
+ * \param elems The mesh to be smoothed
+ * \param fixed_tag The tag marking which vertices are fixed
+ * \param abs_tol Absolute tolerance measuring convergence
+ * \param rel_tol Relative tolerance measuring convergence
+ */
+ LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+ double abs_tol = -1.0, double rel_tol = 1.0e-6);
+
+ /* \brief Destructor
+ */
+ ~LloydSmoother();
+
+ /* \brief perform smoothing operation
+ */
+ ErrorCode perform_smooth();
+
+ /* \brief get instance
+ */
+ Interface *mb_impl() {return mbImpl;}
+
+ /* \brief get/set ParallelComm
+ */
+ ParallelComm *pcomm() {return myPcomm;}
+
+ /* \brief get/set ParallelComm
+ */
+ void pcomm(ParallelComm *pc) {myPcomm = pc;}
+
+ /* \brief get/set elements
+ */
+ Range &elems() {return myElems;}
+
+ /* \brief get/set elements
+ */
+ const Range &elems() const {return myElems;}
+
+ /* \brief get/set fixed tag
+ */
+ Tag fixed_tag() {return fixedTag;}
+
+ /* \brief get/set fixed tag
+ */
+ void fixed_tag(Tag fixed) {fixedTag = fixed;}
+
+ /* \brief get/set tolerance
+ */
+ double abs_tol() {return absTol;}
+
+ /* \brief get/set tolerance
+ */
+ void abs_tol(double tol) {absTol = tol;}
+
+ /* \brief get/set tolerance
+ */
+ double rel_tol() {return relTol;}
+
+ /* \brief get/set tolerance
+ */
+ void rel_tol(double tol) {relTol = tol;}
+
+ /* \brief get/set numIts
+ */
+ int num_its() {return numIts;}
+ void num_its(int num) {numIts = num;}
+
+ /* \brief get/set reportIts
+ */
+ int report_its() {return reportIts;}
+ void report_its(int num) {reportIts = num;}
+
+
+private:
+
+ //- initialize some things in certain cases
+ ErrorCode initialize();
+
+ //- MOAB instance
+ Interface *mbImpl;
+
+ //- ParallelComm
+ ParallelComm *myPcomm;
+
+ //- elements to smooth
+ Range myElems;
+
+ //- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
+ Tag fixedTag;
+
+ //- tolerances
+ double absTol, relTol;
+
+ //- error handler for this class
+ Error *errorHandler;
+
+ //- number of iterations between reporting
+ int reportIts;
+
+ //- number of iterations taken during smooth
+ int numIts;
+
+ //- keep track of whether I created the fixed tag
+ bool iCreatedTag;
+};
+
+} // namespace moab
+
+#endif
diff --git a/test/Makefile.am b/test/Makefile.am
index 9566338..dffb8f9 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -45,7 +45,8 @@ TESTS = range_test \
reorder_test \
test_prog_opt \
coords_connect_iterate \
- test_boundbox
+ test_boundbox \
+ lloyd_smoother_test
if HDF5_FILE
TESTS += mbfacet_test \
@@ -139,6 +140,7 @@ coords_connect_iterate_SOURCES = coords_connect_iterate.cpp
coords_connect_iterate_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS)
test_boundbox_SOURCES = test_boundbox.cpp
+lloyd_smoother_test_SOURCES = lloyd_smoother_test.cpp
if PARALLEL
moab_test_CPPFLAGS += -I$(top_srcdir)/src/parallel
diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
new file mode 100644
index 0000000..8d88112
--- /dev/null
+++ b/test/lloyd_smoother_test.cpp
@@ -0,0 +1,41 @@
+#include "moab/Core.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "TestUtil.hpp"
+
+#ifdef MESHDIR
+std::string TestDir( STRINGIFY(MESHDIR) );
+#else
+std::string TestDir(".");
+#endif
+
+std::string filename = TestDir + "/surfrandomtris-4part.h5m";
+
+using namespace moab;
+
+int main(int argc, char**argv)
+{
+ if (argc > 1) filename = std::string(argv[1]);
+ Core mb;
+ ErrorCode rval = mb.load_file(filename.c_str());
+ CHECK_ERR(rval);
+
+ Range elems;
+ rval = mb.get_entities_by_dimension(0, 3, elems);
+ CHECK_ERR(rval);
+ if (elems.empty()) {
+ rval = mb.get_entities_by_dimension(0, 2, elems);
+ CHECK_ERR(rval);
+ }
+ if (elems.empty()) {
+ std::cout << "Mesh must have faces or regions for this test." << std::endl;
+ CHECK_ERR(MB_FAILURE);
+ }
+
+ LloydSmoother ll(&mb, NULL, elems);
+ ll.report_its(10);
+ rval = ll.perform_smooth();
+ CHECK_ERR(rval);
+
+ std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+}
+
https://bitbucket.org/fathomteam/moab/commits/da78fb058113/
Changeset: da78fb058113
Branch: None
User: tautges
Date: 2013-09-13 00:36:32
Summary: Adding ability to smooth based on a tag instead of coords.
Affected #: 3 files
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index 49177d0..09546b5 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -17,9 +17,9 @@
namespace moab {
-LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
- double at, double rt)
- : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ctag,
+ Tag ftag, double at, double rt)
+ : mbImpl(impl), myPcomm(pc), myElems(elms), coordsTag(ctag), fixedTag(ftag), absTol(at), relTol(rt),
reportIts(0), numIts(0), iCreatedTag(false)
{
ErrorCode rval = mbImpl->query_interface(errorHandler);
@@ -62,7 +62,12 @@ ErrorCode LloydSmoother::perform_smooth()
// get all verts coords into tag; don't need to worry about filtering out fixed verts,
// we'll just be setting to their fixed coords
std::vector<double> vcentroids(3*verts.size());
- rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+ if (!coordsTag) {
+ rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+ }
+ else {
+ rval = mbImpl->tag_get_data(coordsTag, verts, &vcentroids[0]); RR("Failed to get vert coords tag values.");
+ }
Tag centroid;
rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE);
@@ -174,7 +179,13 @@ ErrorCode LloydSmoother::perform_smooth()
} // end while
// write the tag back onto vertex coordinates
- rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+ if (!coordsTag) {
+ rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+ }
+ else {
+ rval = mbImpl->tag_set_data(coordsTag, owned_verts, &vcentroids[0]);
+ RR("Failed to set vert coords tag values.");
+ }
return MB_SUCCESS;
}
diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
index 5ab9771..4f8570a 100644
--- a/src/moab/LloydSmoother.hpp
+++ b/src/moab/LloydSmoother.hpp
@@ -35,13 +35,15 @@ public:
/* \brief Constructor
* Convenience constructor, data input directly
* \param impl The MOAB instance for this smoother
- * \param pcomm The ParallelComm instance by which this mesh is parallel
+ * \param pc The ParallelComm instance by which this mesh is parallel
* \param elems The mesh to be smoothed
+ * \param cds_tag If specified, this tag is used to get/set coordinates, rather than
+ * true vertex coordinates
* \param fixed_tag The tag marking which vertices are fixed
* \param abs_tol Absolute tolerance measuring convergence
* \param rel_tol Relative tolerance measuring convergence
*/
- LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+ LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag cds_tag = 0, Tag fixed_tag = 0,
double abs_tol = -1.0, double rel_tol = 1.0e-6);
/* \brief Destructor
@@ -80,6 +82,14 @@ public:
*/
void fixed_tag(Tag fixed) {fixedTag = fixed;}
+ /* \brief get/set coords tag
+ */
+ Tag coords_tag() {return coordsTag;}
+
+ /* \brief get/set coords tag
+ */
+ void coords_tag(Tag coords) {coordsTag = coords;}
+
/* \brief get/set tolerance
*/
double abs_tol() {return absTol;}
@@ -121,6 +131,9 @@ private:
//- elements to smooth
Range myElems;
+ //- tag for coordinates; if zero, true vertex coords are used
+ Tag coordsTag;
+
//- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
Tag fixedTag;
diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 8d88112..662732b 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,5 +1,6 @@
#include "moab/Core.hpp"
#include "moab/LloydSmoother.hpp"
+#include "moab/CartVect.hpp"
#include "TestUtil.hpp"
#ifdef MESHDIR
@@ -31,11 +32,47 @@ int main(int argc, char**argv)
CHECK_ERR(MB_FAILURE);
}
- LloydSmoother ll(&mb, NULL, elems);
+ // get the vertex positions and set on an intermediate tag, to test smoothing for
+ // a tag instead of for coords
+ Tag ctag;
+ rval = mb.tag_get_handle("vcentroid", 3, MB_TYPE_DOUBLE, ctag, MB_TAG_CREAT|MB_TAG_DENSE);
+ CHECK_ERR(rval);
+ Range verts;
+ rval = mb.get_entities_by_dimension(0, 0, verts);
+ CHECK_ERR(rval);
+ std::vector<double> coords(3*verts.size());
+ rval = mb.get_coords(verts, &coords[0]);
+ CHECK_ERR(rval);
+ rval = mb.tag_set_data(ctag, verts, &coords[0]);
+ CHECK_ERR(rval);
+
+ LloydSmoother ll(&mb, NULL, elems, ctag);
ll.report_its(10);
rval = ll.perform_smooth();
CHECK_ERR(rval);
-
std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+
+ // now, set vertex coords to almost their converged positions, then re-smooth; should take fewer
+ // iterations
+ std::vector<double> new_coords(3*verts.size());
+ rval = mb.tag_get_data(ctag, verts, &new_coords[0]);
+ CHECK_ERR(rval);
+ unsigned int i;
+ Range::iterator vit;
+ for (vit = verts.begin(), i = 0; vit != verts.end(); vit++, i+=3) {
+ CartVect old_pos(&coords[i]), new_pos(&new_coords[i]);
+ CartVect almost_pos = old_pos + .99 * (new_pos - old_pos);
+ almost_pos.get(&new_coords[i]);
+ }
+ rval = mb.set_coords(verts, &new_coords[0]);
+ CHECK_ERR(rval);
+
+ LloydSmoother ll2(&mb, NULL, elems);
+ ll2.report_its(10);
+ rval = ll2.perform_smooth();
+ CHECK_ERR(rval);
+ std::cout << "Mesh smoothed in " << ll2.num_its() << " iterations." << std::endl;
+
+
}
https://bitbucket.org/fathomteam/moab/commits/c1d28fa90048/
Changeset: c1d28fa90048
Branch: None
User: tautges
Date: 2013-09-17 00:29:43
Summary: Adding deformed mesh example.
Affected #: 2 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
new file mode 100644
index 0000000..ddb3f3a
--- /dev/null
+++ b/examples/DeformMeshRemap.cpp
@@ -0,0 +1,180 @@
+/** @example DeformMeshRemap.cpp
+ * Description: Account for mesh deformation of a solid due to structural mechanics\n
+ * In this example there are two meshes, a "master" and "slave" mesh. In the master mesh,
+ * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
+ * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
+ * remapped to those new positions. Then mesh positions and state variables are transferred
+ * to the slave mesh, mimicing another mesh used by some other physics.
+ *
+ * To run: ./DeformMeshRemap [<master_meshfile><slave_meshfile>]\n
+ * (default values can run if users don't specify the mesh files)
+ */
+
+#include "moab/Core.hpp"
+#include "moab/Range.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "MBTagConventions.hpp"
+
+#include <iostream>
+#include <assert.h>
+
+using namespace moab;
+using namespace std;
+
+#ifndef MESH_DIR
+#define MESH_DIR "."
+#endif
+
+ErrorCode read_file(string &fname, EntityHandle &seth,
+ Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems);
+void deform_func(double *xold, double *xnew);
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
+ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
+ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+string trimeshf = string(MESH_DIR) + string("/rodtri.g");
+string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
+const int SOLID_SETNO = 100, FLUID_SETNO = 200;
+
+Interface *mb;
+#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
+
+
+int main(int /*argc*/, char **/*argv*/) {
+
+ EntityHandle master, slave;
+ ErrorCode rval;
+
+ mb = new Core();
+
+ // read master/slave files and get fluid/solid material sets
+ Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
+ rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+ rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+
+ // deform the master's solid mesh, put results in a new tag
+ Tag xnew;
+ rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
+ rval = write_to_coords(solid_elems[0], xnew); RR("");
+ rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+
+ // smooth the master mesh
+ LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
+ rval = ll->perform_smooth();
+ RR("Failed in lloyd smoothing.");
+ cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+
+ rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+
+ delete ll;
+ delete mb;
+
+ return MB_SUCCESS;
+}
+
+ErrorCode write_to_coords(Range &elems, Tag tagh)
+{
+ // write the tag to coordinates
+ Range verts;
+ ErrorCode rval = mb->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get adj vertices.");
+ std::vector<double> coords(3*verts.size());
+ rval = mb->tag_get_data(tagh, verts, &coords[0]);
+ RR("Failed to get tag data.");
+ rval = mb->set_coords(verts, &coords[0]);
+ RR("Failed to set coordinates.");
+ return MB_SUCCESS;
+}
+
+void deform_func(double *xold, double *xnew)
+{
+ const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
+ // function: origin is at middle base of rod, and is .5 high
+ // top of rod is (0,.55) on left and (.2,.6) on right
+ double delx = 0.5*RODWIDTH;
+
+ double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
+ xnew[0] = xold[0] + yfrac * delx;
+ xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
+}
+
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew)
+{
+ // deform elements with an analytic function
+
+ // create the tag
+ ErrorCode rval = mb->tag_get_handle("", 3, MB_TYPE_DOUBLE, xnew, MB_TAG_CREAT|MB_TAG_DENSE);
+ RR("Failed to create xnew tag.");
+
+ // get all the vertices and coords in the fluid, set xnew to them
+ Range verts;
+ rval = mb->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ std::vector<double> coords(3*verts.size(), 0.0);
+ rval = mb->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ RR("Failed to set xnew tag on fluid verts.");
+
+ // get all the vertices and coords in the solid
+ verts.clear();
+ rval = mb->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ coords.resize(3*verts.size(), 0.0);
+ rval = mb->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ unsigned int num_verts = verts.size();
+ for (unsigned int i = 0; i < num_verts; i++)
+ deform_func(&coords[3*i], &coords[3*i]);
+
+ // set the new tag to those coords
+ rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ RR("Failed to set tag data.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode read_file(string &fname, EntityHandle &seth,
+ Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems)
+{
+ // create meshset
+ ErrorCode rval = mb->create_meshset(0, seth);
+ RR("Couldn't create master/slave set.");
+ rval = mb->load_file(fname.c_str(), &seth);
+ RR("Couldn't load master/slave mesh.");
+
+ // get material sets for solid/fluid
+ Tag tagh;
+ rval = mb->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+ const void *setno_ptr = &SOLID_SETNO;
+ rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, solids);
+ if (solids.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any solid sets.");
+
+ // get solid entities, and dimension
+ Range tmp_range;
+ for (Range::iterator rit = solids.begin(); rit != solids.end(); rit++) {
+ rval = mb->get_entities_by_handle(*rit, tmp_range, true);
+ RR("Failed to get entities in solid.");
+ }
+ int dim = mb->dimension_from_handle(*tmp_range.rbegin());
+ assert(dim > 0 && dim < 4);
+
+ solid_elems = tmp_range.subset_by_dimension(dim);
+
+ setno_ptr = &FLUID_SETNO;
+ rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, fluids);
+ if (fluids.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any fluid sets.");
+
+ for (Range::iterator rit = fluids.begin(); rit != fluids.end(); rit++) {
+ rval = mb->get_entities_by_dimension(*rit, dim, fluid_elems, true);
+ RR("Failed to get entities in fluid.");
+ }
+ if (mb->dimension_from_handle(*fluid_elems.begin()) != dim) {
+ rval = MB_FAILURE;
+ RR("Fluid and solid elements must be same dimension.");
+ }
+
+ return MB_SUCCESS;
+}
diff --git a/examples/makefile b/examples/makefile
index 3b848d2..fbd4412 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -4,7 +4,7 @@ include ${MOAB_DIR}/lib/iMesh-Defs.inc
.SUFFIXES: .o .cpp .F90
-# MESH_DIR is the top-level MOAB source directory, used to locate mesh files that come with MOAB source
+# MESH_DIR is the directory containing mesh files that come with MOAB source
MESH_DIR="../MeshFiles/unittest"
EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles
@@ -47,6 +47,9 @@ HelloParMOAB: HelloParMOAB.o
TestExodusII: TestExodusII.o
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+DeformMeshRemap: DeformMeshRemap.o
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+
clean:
rm -rf *.o ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES}
https://bitbucket.org/fathomteam/moab/commits/a8f966f4e3ff/
Changeset: a8f966f4e3ff
Branch: None
User: tautges
Date: 2013-09-18 00:44:23
Summary: Adding option processing to DeformMeshRemap, and making
parallel optional in LloydRelaxation example.
Affected #: 2 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index ddb3f3a..55ac3cc 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -13,6 +13,7 @@
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/LloydSmoother.hpp"
+#include "moab/ProgOptions.hpp"
#include "MBTagConventions.hpp"
#include <iostream>
@@ -32,39 +33,48 @@ ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
ErrorCode write_to_coords(Range &elems, Tag tagh);
-string trimeshf = string(MESH_DIR) + string("/rodtri.g");
-string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
const int SOLID_SETNO = 100, FLUID_SETNO = 200;
Interface *mb;
#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
-int main(int /*argc*/, char **/*argv*/) {
+int main(int argc, char **argv) {
EntityHandle master, slave;
ErrorCode rval;
+ ProgOptions po("Deformed mesh options");
+ po.addOpt<std::string> ("master,m", "Specify the master meshfile name" );
+ po.addOpt<std::string> ("slave,s", "Specify the slave meshfile name" );
+ po.parseCommandLine(argc, argv);
+ std::string foo;
+ string masterf, slavef;
+ if(!po.getOpt("master", &masterf))
+ masterf = string(MESH_DIR) + string("/rodquad.g");
+ if(!po.getOpt("slave", &slavef))
+ slavef = string(MESH_DIR) + string("/rodtri.g");
+
mb = new Core();
// read master/slave files and get fluid/solid material sets
Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
- rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
- rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+ rval = read_file(masterf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+ rval = read_file(slavef, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
// deform the master's solid mesh, put results in a new tag
Tag xnew;
rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
- rval = write_to_coords(solid_elems[0], xnew); RR("");
- rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+ if (debug) write_and_save(solid_elems[0], master, xnew, "deformed.vtk");
// smooth the master mesh
LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
rval = ll->perform_smooth();
RR("Failed in lloyd smoothing.");
cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+ if (debug) write_and_save(fluid_elems[0], master, xnew, "smoothed.vtk");
- rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+ // map new locations to slave
delete ll;
delete mb;
@@ -72,6 +82,13 @@ int main(int /*argc*/, char **/*argv*/) {
return MB_SUCCESS;
}
+ErrorCode write_and_save(Range &ents, EntithHanlde seth, Tag tagh, const char *filename)
+{
+ rval = write_to_coords(ents, tagh); RR("");
+ rval = mb->write_file("deformed.vtk", NULL, NULL, &seth, 1); RR("");
+ return rval;
+}
+
ErrorCode write_to_coords(Range &elems, Tag tagh)
{
// write the tag to coordinates
diff --git a/examples/LloydRelaxation.cpp b/examples/LloydRelaxation.cpp
index 0cf9d1e..34002cd 100644
--- a/examples/LloydRelaxation.cpp
+++ b/examples/LloydRelaxation.cpp
@@ -14,8 +14,10 @@
* in the current directory (H5M format must be used since the file is written in parallel).
*/
-#include "moab/ParallelComm.hpp"
-#include "MBParallelConventions.h"
+#ifdef USE_MPI
+# include "moab/ParallelComm.hpp"
+# include "MBParallelConventions.h"
+#endif
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
@@ -30,7 +32,7 @@ string test_file_name = string(MESH_DIR) + string("/surfrandomtris-4part.h5m");
#define RC if (MB_SUCCESS != rval) return rval
-ErrorCode perform_lloyd_relaxation(ParallelComm *pc, Range &verts, Range &cells, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &cells, Tag fixed,
int num_its, int report_its);
int main(int argc, char **argv)
@@ -38,7 +40,9 @@ int main(int argc, char **argv)
int num_its = 10;
int report_its = 1;
+#ifdef USE_MPI
MPI_Init(&argc, &argv);
+#endif
// need option handling here for input filename
if (argc > 1){
@@ -49,10 +53,13 @@ int main(int argc, char **argv)
// get MOAB and ParallelComm instances
Interface *mb = new Core;
if (NULL == mb) return 1;
+ int nprocs = 1;
+
+#ifdef USE_MPI
// get the ParallelComm instance
- ParallelComm* pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
- int nprocs = pcomm->size();
-
+ ParallelComm *pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
+ nprocs = pcomm->size();
+#endif
string options;
if (nprocs > 1) // if reading in parallel, need to tell it how
options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
@@ -76,34 +83,45 @@ int main(int argc, char **argv)
// ok to mark non-owned skin vertices too, I won't move those anyway
// use MOAB's skinner class to find the skin
Skinner skinner(mb);
- rval = skinner.find_skin(faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
+ rval = skinner.find_skin(0, faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
std::vector<int> fix_tag(skin_verts.size(), 1); // initialized to 1 to indicate fixed
rval = mb->tag_set_data(fixed, skin_verts, &fix_tag[0]); RC;
// now perform the Lloyd relaxation
- rval = perform_lloyd_relaxation(pcomm, verts, faces, fixed, num_its, report_its); RC;
+ rval = perform_lloyd_relaxation(mb, verts, faces, fixed, num_its, report_its); RC;
// delete fixed tag, since we created it here
rval = mb->tag_delete(fixed); RC;
// output file, using parallel write
- rval = mb->write_file("lloydfinal.h5m", NULL, "PARALLEL=WRITE_PART"); RC;
+
+#ifdef USE_MPI
+ options = "PARALLEL=WRITE_PART";
+#endif
+
+ rval = mb->write_file("lloydfinal.h5m", NULL, options.c_str()); RC;
// delete MOAB instance
delete mb;
+#ifdef USE_MPI
MPI_Finalize();
+#endif
return 0;
}
-ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &faces, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &faces, Tag fixed,
int num_its, int report_its)
{
ErrorCode rval;
- Interface *mb = pcomm->get_moab();
- int nprocs = pcomm->size();
+ int nprocs = 1;
+
+#ifdef USE_MPI
+ ParallelComm *pcomm = ParallelComm::get_pcomm(mb, 0);
+ nprocs = pcomm->size();
+#endif
// perform Lloyd relaxation:
// 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
@@ -120,8 +138,10 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// filter verts down to owned ones and get fixed tag for them
Range owned_verts, shared_owned_verts;
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
if (rval != MB_SUCCESS) return rval;
+#endif
}
else
owned_verts = verts;
@@ -133,11 +153,13 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
vcentroids.resize(3*owned_verts.size());
rval = mb->tag_get_data(centroid, owned_verts, &vcentroids[0]); RC;
+#ifdef USE_MPI
// get shared owned verts, for exchanging tags
rval = pcomm->get_shared_entities(-1, shared_owned_verts, 0, false, true); RC;
// workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
// for all shared entities
if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+#endif
// some declarations for later iterations
std::vector<double> fcentroids(3*faces.size()); // fcentroids for face centroids
@@ -195,16 +217,22 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// 2c. exchange tags on owned verts
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->exchange_tags(centroid, shared_owned_verts); RC;
+#endif
}
if (!(nit%report_its)) {
// global reduce for maximum delta, then report it
double global_max = mxdelta;
+ int myrank = 0;
+#ifdef USE_MPI
if (nprocs > 1)
MPI_Reduce(&mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm());
- if (1 == nprocs || !pcomm->rank())
+ myrank = pcomm->rank();
+#endif
+ if (1 == nprocs || !myrank)
cout << "Max delta = " << global_max << endl;
}
}
https://bitbucket.org/fathomteam/moab/commits/13f6ad9ea2c7/
Changeset: 13f6ad9ea2c7
Branch: None
User: tautges
Date: 2013-09-18 01:09:49
Summary: Merge branch '00c578e' into deformed-mesh-remap
Conflicts:
src/Makefile.am
test/Makefile.am
Affected #: 6 files
diff --git a/.gitignore b/.gitignore
index d3721c3..5aac80a 100644
--- a/.gitignore
+++ b/.gitignore
@@ -155,6 +155,7 @@ test/io/vtk_test
test/kd_tree_test
test/kd_tree_time
test/kd_tree_tool
+test/lloyd_smoother_test
test/mbcn_test
test/mbfacet_test
test/mbground_test
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
new file mode 100644
index 0000000..49177d0
--- /dev/null
+++ b/src/LloydSmoother.cpp
@@ -0,0 +1,220 @@
+#include "moab/LloydSmoother.hpp"
+#include "moab/Error.hpp"
+#include "moab/Skinner.hpp"
+#include "moab/CN.hpp"
+#include "moab/CartVect.hpp"
+#include "moab/BoundBox.hpp"
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
+#include <iostream>
+
+#define RR(a) if (MB_SUCCESS != rval) { \
+ errorHandler->set_last_error(a); \
+ return rval;}
+
+namespace moab {
+
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
+ double at, double rt)
+ : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+ reportIts(0), numIts(0), iCreatedTag(false)
+{
+ ErrorCode rval = mbImpl->query_interface(errorHandler);
+ if (rval) throw rval;
+}
+
+LloydSmoother::~LloydSmoother()
+{
+ if (iCreatedTag && fixedTag) {
+ ErrorCode rval = mbImpl->tag_delete(fixedTag);
+ if (rval) throw rval;
+ }
+}
+
+ErrorCode LloydSmoother::perform_smooth()
+{
+ ErrorCode rval;
+
+ // first figure out tolerance to use
+ if (0 > absTol) {
+ // no tolerance set - get one relative to bounding box around elements
+ BoundBox bb;
+ rval = bb.update(*mbImpl, myElems);
+ RR("Failed to compute bounding box around elements.");
+ absTol = relTol * bb.diagonal_length();
+ }
+
+ // initialize if we need to
+ rval = initialize();
+ RR("Failed to initialize.");
+
+ // get all vertices
+ Range verts;
+ rval = mbImpl->get_adjacencies(myElems, 0, false, verts, Interface::UNION);
+ RR("Failed to get all vertices.");
+
+ // perform Lloyd relaxation:
+ // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
+
+ // get all verts coords into tag; don't need to worry about filtering out fixed verts,
+ // we'll just be setting to their fixed coords
+ std::vector<double> vcentroids(3*verts.size());
+ rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+
+ Tag centroid;
+ rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE);
+ RR("Couldn't get tag handle.");
+ rval = mbImpl->tag_set_data(centroid, verts, &vcentroids[0]); RR("Failed setting centroid tag.");
+
+ Range owned_verts, shared_owned_verts;
+#ifdef USE_MPI
+ // filter verts down to owned ones and get fixed tag for them
+ if (myPcomm->size() > 1) {
+ rval = myPcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
+ RR("Failed to filter on pstatus.");
+ // get shared owned verts, for exchanging tags
+ rval = myPcomm->filter_pstatus(owned_verts, PSTATUS_SHARED, PSTATUS_AND, -1, &shared_owned_verts);
+ RR("Failed to filter for shared owned.");
+ // workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
+ // for all shared entities
+ if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+ }
+ else
+ owned_verts = verts;
+#else
+ owned_verts = verts;
+#endif
+
+ std::vector<unsigned char> fix_tag(owned_verts.size());
+ rval = mbImpl->tag_get_data(fixedTag, owned_verts, &fix_tag[0]); RR("Failed to get fixed tag.");
+
+ // now fill vcentroids array with positions of just owned vertices, since those are the ones
+ // we're actually computing
+ vcentroids.resize(3*owned_verts.size());
+ rval = mbImpl->tag_get_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to get centroid tag.");
+
+ // some declarations for later iterations
+ std::vector<double> fcentroids(3*myElems.size()); // fcentroids for element centroids
+ std::vector<double> ctag(3*CN::MAX_NODES_PER_ELEMENT); // temporary coordinate storage for verts bounding an element
+ const EntityHandle *conn; // const ptr & size to elem connectivity
+ int nconn;
+ Range::iterator eit, vit; // for iterating over elems, verts
+ int e, v; // for indexing into centroid vectors
+ std::vector<EntityHandle> adj_elems; // used in vertex iteration
+
+ // 2. while !converged
+ double resid = DBL_MAX;
+ numIts = 0;
+ while (resid > absTol) {
+ numIts++;
+ resid = 0.0;
+
+ // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
+ for (eit = myElems.begin(), e = 0; eit != myElems.end(); eit++, e++) {
+ // get verts for this elem
+ rval = mbImpl->get_connectivity(*eit, conn, nconn); RR("Failed to get connectivity.");
+ // get centroid tags for those verts
+ rval = mbImpl->tag_get_data(centroid, conn, nconn, &ctag[0]); RR("Failed to get centroid.");
+ fcentroids[3*e+0] = fcentroids[3*e+1] = fcentroids[3*e+2] = 0.0;
+ for (v = 0; v < nconn; v++) {
+ fcentroids[3*e+0] += ctag[3*v+0];
+ fcentroids[3*e+1] += ctag[3*v+1];
+ fcentroids[3*e+2] += ctag[3*v+2];
+ }
+ for (v = 0; v < 3; v++) fcentroids[3*e+v] /= nconn;
+ }
+ rval = mbImpl->tag_set_data(centroid, myElems, &fcentroids[0]); RR("Failed to set elem centroid.");
+
+ // 2b. foreach owned vertex:
+ for (vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); vit++, v++) {
+ // if !fixed
+ if (fix_tag[v]) continue;
+ // vertex centroid = sum(cell centroids)/ncells
+ adj_elems.clear();
+ rval = mbImpl->get_adjacencies(&(*vit), 1, 2, false, adj_elems); RR("Failed getting adjs.");
+ rval = mbImpl->tag_get_data(centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0]);
+ RR("Failed to get elem centroid.");
+ double vnew[] = {0.0, 0.0, 0.0};
+ for (e = 0; e < (int)adj_elems.size(); e++) {
+ vnew[0] += fcentroids[3*e+0];
+ vnew[1] += fcentroids[3*e+1];
+ vnew[2] += fcentroids[3*e+2];
+ }
+ for (e = 0; e < 3; e++) vnew[e] /= adj_elems.size();
+ double delta = (CartVect(vnew)-CartVect(&vcentroids[3*v])).length();
+ resid = (v ? std::max(resid, delta) : delta);
+ for (e = 0; e < 3; e++) vcentroids[3*v+e] = vnew[e];
+ }
+
+ // set the centroid tag; having them only in vcentroids array isn't enough, as vertex centroids are
+ // accessed randomly in loop over faces
+ rval = mbImpl->tag_set_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to set vertex centroid.");
+
+#ifdef USE_MPI
+ // 2c. exchange tags on owned verts
+ if (myPcomm->size() > 1) {
+ rval = pcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
+ }
+#endif
+
+ if (reportIts && !(numIts%reportIts)) {
+ double global_max = resid;
+#ifdef USE_MPI
+ // global reduce for maximum delta, then report it
+ if (myPcomm->size() > 1)
+ MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
+ if (!pcomm->rank())
+#endif
+ std::cout << "Max residual = " << global_max << std::endl;
+ }
+
+ } // end while
+
+ // write the tag back onto vertex coordinates
+ rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode LloydSmoother::initialize()
+{
+ ErrorCode rval = MB_SUCCESS;
+
+ // first, check for tag; if we don't have it, make one and mark skin as fixed
+ if (!fixedTag) {
+ unsigned char fixed = 0x0;
+ rval = mbImpl->tag_get_handle("", 1, MB_TYPE_OPAQUE, fixedTag, MB_TAG_DENSE|MB_TAG_CREAT, &fixed);
+ RR("Trouble making fixed tag.");
+ iCreatedTag = true;
+
+ // get the skin; get facets, because we might need to filter on shared entities
+ Skinner skinner(mbImpl);
+ Range skin, skin_verts;
+ rval = skinner.find_skin(0, myElems, false, skin);
+ RR("Unable to find skin.");
+
+#ifdef USE_MPI
+ // need to do a little extra if we're working in parallel
+ if (myPcomm) {
+ // filter out ghost and interface facets
+ rval = myPcomm->filter_pstatus(skin, PSTATUS_GHOST|PSTATUS_INTERFACE, PSTATUS_NOT);
+ RR("Failed to filter on shared status.");
+ }
+#endif
+ // get the vertices from those entities
+ rval = mbImpl->get_adjacencies(skin, 0, false, skin_verts, Interface::UNION);
+ RR("Trouble getting vertices.");
+
+ // mark them fixed
+ std::vector<unsigned char> marks(skin.size(), 0x1);
+ rval = mbImpl->tag_set_data(fixedTag, skin_verts, &marks[0]);
+ RR("Unable to set tag on skin.");
+ }
+
+ return MB_SUCCESS;
+}
+
+}
diff --git a/src/Makefile.am b/src/Makefile.am
index dc0c3f2..4249153 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -65,6 +65,7 @@ libMOAB_la_SOURCES = \
HigherOrderFactory.cpp \
HomXform.cpp \
Internals.hpp \
+ LloydSmoother.cpp \
MBCNArrays.hpp \
MergeMesh.cpp \
MeshSet.cpp \
@@ -167,6 +168,7 @@ nobase_libMOAB_la_include_HEADERS = \
moab/Forward.hpp \
moab/GeomUtil.hpp \
moab/Interface.hpp \
+ moab/LloydSmoother.hpp \
moab/point_locater/tree/common_tree.hpp \
moab/point_locater/tree/element_tree.hpp \
moab/point_locater/tree/bvh_tree.hpp \
diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
new file mode 100644
index 0000000..5ab9771
--- /dev/null
+++ b/src/moab/LloydSmoother.hpp
@@ -0,0 +1,145 @@
+/** @class LloydSmoother.cpp \n
+ * \brief Perform Lloyd relaxation on a mesh and its dual \n
+ *
+ * Briefly, Lloyd relaxation is a technique to smooth out a mesh. The centroid of each cell is computed from its
+ * vertex positions, then vertices are placed at the average of their connected cells' centroids, and the process
+ * iterates until convergence.
+ *
+ * In the parallel algorithm, an extra ghost layer of cells is exchanged. This allows us to compute the centroids
+ * for boundary cells on each processor where they appear; this eliminates the need for one round of data exchange
+ * (for those centroids) between processors. New vertex positions must be sent from owning processors to processors
+ * sharing those vertices. Convergence is measured as the maximum distance moved by any vertex.
+ *
+ */
+
+#ifndef LLOYDSMOOTHER_HPP
+#define LLOYDSMOOTHER_HPP
+
+#include "moab/Interface.hpp"
+
+namespace moab {
+
+class ParallelComm;
+class Error;
+
+class LloydSmoother
+{
+public:
+
+ /* \brief Constructor
+ * Bare constructor, data input to this class through methods.
+ * \param impl The MOAB instance for this smoother
+ */
+ LloydSmoother(Interface *impl);
+
+ /* \brief Constructor
+ * Convenience constructor, data input directly
+ * \param impl The MOAB instance for this smoother
+ * \param pcomm The ParallelComm instance by which this mesh is parallel
+ * \param elems The mesh to be smoothed
+ * \param fixed_tag The tag marking which vertices are fixed
+ * \param abs_tol Absolute tolerance measuring convergence
+ * \param rel_tol Relative tolerance measuring convergence
+ */
+ LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+ double abs_tol = -1.0, double rel_tol = 1.0e-6);
+
+ /* \brief Destructor
+ */
+ ~LloydSmoother();
+
+ /* \brief perform smoothing operation
+ */
+ ErrorCode perform_smooth();
+
+ /* \brief get instance
+ */
+ Interface *mb_impl() {return mbImpl;}
+
+ /* \brief get/set ParallelComm
+ */
+ ParallelComm *pcomm() {return myPcomm;}
+
+ /* \brief get/set ParallelComm
+ */
+ void pcomm(ParallelComm *pc) {myPcomm = pc;}
+
+ /* \brief get/set elements
+ */
+ Range &elems() {return myElems;}
+
+ /* \brief get/set elements
+ */
+ const Range &elems() const {return myElems;}
+
+ /* \brief get/set fixed tag
+ */
+ Tag fixed_tag() {return fixedTag;}
+
+ /* \brief get/set fixed tag
+ */
+ void fixed_tag(Tag fixed) {fixedTag = fixed;}
+
+ /* \brief get/set tolerance
+ */
+ double abs_tol() {return absTol;}
+
+ /* \brief get/set tolerance
+ */
+ void abs_tol(double tol) {absTol = tol;}
+
+ /* \brief get/set tolerance
+ */
+ double rel_tol() {return relTol;}
+
+ /* \brief get/set tolerance
+ */
+ void rel_tol(double tol) {relTol = tol;}
+
+ /* \brief get/set numIts
+ */
+ int num_its() {return numIts;}
+ void num_its(int num) {numIts = num;}
+
+ /* \brief get/set reportIts
+ */
+ int report_its() {return reportIts;}
+ void report_its(int num) {reportIts = num;}
+
+
+private:
+
+ //- initialize some things in certain cases
+ ErrorCode initialize();
+
+ //- MOAB instance
+ Interface *mbImpl;
+
+ //- ParallelComm
+ ParallelComm *myPcomm;
+
+ //- elements to smooth
+ Range myElems;
+
+ //- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
+ Tag fixedTag;
+
+ //- tolerances
+ double absTol, relTol;
+
+ //- error handler for this class
+ Error *errorHandler;
+
+ //- number of iterations between reporting
+ int reportIts;
+
+ //- number of iterations taken during smooth
+ int numIts;
+
+ //- keep track of whether I created the fixed tag
+ bool iCreatedTag;
+};
+
+} // namespace moab
+
+#endif
diff --git a/test/Makefile.am b/test/Makefile.am
index e4d591f..622c86c 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -48,7 +48,8 @@ TESTS = range_test \
coords_connect_iterate \
elem_eval_test \
spatial_locator_test \
- test_boundbox
+ test_boundbox \
+ lloyd_smoother_test
if HDF5_FILE
TESTS += mbfacet_test \
@@ -144,6 +145,7 @@ coords_connect_iterate_SOURCES = coords_connect_iterate.cpp
coords_connect_iterate_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS)
test_boundbox_SOURCES = test_boundbox.cpp
+lloyd_smoother_test_SOURCES = lloyd_smoother_test.cpp
if PARALLEL
moab_test_CPPFLAGS += -I$(top_srcdir)/src/parallel
diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
new file mode 100644
index 0000000..8d88112
--- /dev/null
+++ b/test/lloyd_smoother_test.cpp
@@ -0,0 +1,41 @@
+#include "moab/Core.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "TestUtil.hpp"
+
+#ifdef MESHDIR
+std::string TestDir( STRINGIFY(MESHDIR) );
+#else
+std::string TestDir(".");
+#endif
+
+std::string filename = TestDir + "/surfrandomtris-4part.h5m";
+
+using namespace moab;
+
+int main(int argc, char**argv)
+{
+ if (argc > 1) filename = std::string(argv[1]);
+ Core mb;
+ ErrorCode rval = mb.load_file(filename.c_str());
+ CHECK_ERR(rval);
+
+ Range elems;
+ rval = mb.get_entities_by_dimension(0, 3, elems);
+ CHECK_ERR(rval);
+ if (elems.empty()) {
+ rval = mb.get_entities_by_dimension(0, 2, elems);
+ CHECK_ERR(rval);
+ }
+ if (elems.empty()) {
+ std::cout << "Mesh must have faces or regions for this test." << std::endl;
+ CHECK_ERR(MB_FAILURE);
+ }
+
+ LloydSmoother ll(&mb, NULL, elems);
+ ll.report_its(10);
+ rval = ll.perform_smooth();
+ CHECK_ERR(rval);
+
+ std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+}
+
https://bitbucket.org/fathomteam/moab/commits/5048df95db0f/
Changeset: 5048df95db0f
Branch: None
User: tautges
Date: 2013-09-18 01:11:03
Summary: Merge branch 'da78fb058' into deformed-mesh-remap
Affected #: 3 files
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index 49177d0..09546b5 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -17,9 +17,9 @@
namespace moab {
-LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
- double at, double rt)
- : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ctag,
+ Tag ftag, double at, double rt)
+ : mbImpl(impl), myPcomm(pc), myElems(elms), coordsTag(ctag), fixedTag(ftag), absTol(at), relTol(rt),
reportIts(0), numIts(0), iCreatedTag(false)
{
ErrorCode rval = mbImpl->query_interface(errorHandler);
@@ -62,7 +62,12 @@ ErrorCode LloydSmoother::perform_smooth()
// get all verts coords into tag; don't need to worry about filtering out fixed verts,
// we'll just be setting to their fixed coords
std::vector<double> vcentroids(3*verts.size());
- rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+ if (!coordsTag) {
+ rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+ }
+ else {
+ rval = mbImpl->tag_get_data(coordsTag, verts, &vcentroids[0]); RR("Failed to get vert coords tag values.");
+ }
Tag centroid;
rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE);
@@ -174,7 +179,13 @@ ErrorCode LloydSmoother::perform_smooth()
} // end while
// write the tag back onto vertex coordinates
- rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+ if (!coordsTag) {
+ rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+ }
+ else {
+ rval = mbImpl->tag_set_data(coordsTag, owned_verts, &vcentroids[0]);
+ RR("Failed to set vert coords tag values.");
+ }
return MB_SUCCESS;
}
diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
index 5ab9771..4f8570a 100644
--- a/src/moab/LloydSmoother.hpp
+++ b/src/moab/LloydSmoother.hpp
@@ -35,13 +35,15 @@ public:
/* \brief Constructor
* Convenience constructor, data input directly
* \param impl The MOAB instance for this smoother
- * \param pcomm The ParallelComm instance by which this mesh is parallel
+ * \param pc The ParallelComm instance by which this mesh is parallel
* \param elems The mesh to be smoothed
+ * \param cds_tag If specified, this tag is used to get/set coordinates, rather than
+ * true vertex coordinates
* \param fixed_tag The tag marking which vertices are fixed
* \param abs_tol Absolute tolerance measuring convergence
* \param rel_tol Relative tolerance measuring convergence
*/
- LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+ LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag cds_tag = 0, Tag fixed_tag = 0,
double abs_tol = -1.0, double rel_tol = 1.0e-6);
/* \brief Destructor
@@ -80,6 +82,14 @@ public:
*/
void fixed_tag(Tag fixed) {fixedTag = fixed;}
+ /* \brief get/set coords tag
+ */
+ Tag coords_tag() {return coordsTag;}
+
+ /* \brief get/set coords tag
+ */
+ void coords_tag(Tag coords) {coordsTag = coords;}
+
/* \brief get/set tolerance
*/
double abs_tol() {return absTol;}
@@ -121,6 +131,9 @@ private:
//- elements to smooth
Range myElems;
+ //- tag for coordinates; if zero, true vertex coords are used
+ Tag coordsTag;
+
//- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
Tag fixedTag;
diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 8d88112..662732b 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,5 +1,6 @@
#include "moab/Core.hpp"
#include "moab/LloydSmoother.hpp"
+#include "moab/CartVect.hpp"
#include "TestUtil.hpp"
#ifdef MESHDIR
@@ -31,11 +32,47 @@ int main(int argc, char**argv)
CHECK_ERR(MB_FAILURE);
}
- LloydSmoother ll(&mb, NULL, elems);
+ // get the vertex positions and set on an intermediate tag, to test smoothing for
+ // a tag instead of for coords
+ Tag ctag;
+ rval = mb.tag_get_handle("vcentroid", 3, MB_TYPE_DOUBLE, ctag, MB_TAG_CREAT|MB_TAG_DENSE);
+ CHECK_ERR(rval);
+ Range verts;
+ rval = mb.get_entities_by_dimension(0, 0, verts);
+ CHECK_ERR(rval);
+ std::vector<double> coords(3*verts.size());
+ rval = mb.get_coords(verts, &coords[0]);
+ CHECK_ERR(rval);
+ rval = mb.tag_set_data(ctag, verts, &coords[0]);
+ CHECK_ERR(rval);
+
+ LloydSmoother ll(&mb, NULL, elems, ctag);
ll.report_its(10);
rval = ll.perform_smooth();
CHECK_ERR(rval);
-
std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+
+ // now, set vertex coords to almost their converged positions, then re-smooth; should take fewer
+ // iterations
+ std::vector<double> new_coords(3*verts.size());
+ rval = mb.tag_get_data(ctag, verts, &new_coords[0]);
+ CHECK_ERR(rval);
+ unsigned int i;
+ Range::iterator vit;
+ for (vit = verts.begin(), i = 0; vit != verts.end(); vit++, i+=3) {
+ CartVect old_pos(&coords[i]), new_pos(&new_coords[i]);
+ CartVect almost_pos = old_pos + .99 * (new_pos - old_pos);
+ almost_pos.get(&new_coords[i]);
+ }
+ rval = mb.set_coords(verts, &new_coords[0]);
+ CHECK_ERR(rval);
+
+ LloydSmoother ll2(&mb, NULL, elems);
+ ll2.report_its(10);
+ rval = ll2.perform_smooth();
+ CHECK_ERR(rval);
+ std::cout << "Mesh smoothed in " << ll2.num_its() << " iterations." << std::endl;
+
+
}
https://bitbucket.org/fathomteam/moab/commits/a1b3900d27b0/
Changeset: a1b3900d27b0
Branch: None
User: tautges
Date: 2013-09-18 01:12:54
Summary: Merge branch 'c1d28fa' into deformed-mesh-remap
Conflicts:
examples/makefile
Affected #: 2 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
new file mode 100644
index 0000000..ddb3f3a
--- /dev/null
+++ b/examples/DeformMeshRemap.cpp
@@ -0,0 +1,180 @@
+/** @example DeformMeshRemap.cpp
+ * Description: Account for mesh deformation of a solid due to structural mechanics\n
+ * In this example there are two meshes, a "master" and "slave" mesh. In the master mesh,
+ * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
+ * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
+ * remapped to those new positions. Then mesh positions and state variables are transferred
+ * to the slave mesh, mimicing another mesh used by some other physics.
+ *
+ * To run: ./DeformMeshRemap [<master_meshfile><slave_meshfile>]\n
+ * (default values can run if users don't specify the mesh files)
+ */
+
+#include "moab/Core.hpp"
+#include "moab/Range.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "MBTagConventions.hpp"
+
+#include <iostream>
+#include <assert.h>
+
+using namespace moab;
+using namespace std;
+
+#ifndef MESH_DIR
+#define MESH_DIR "."
+#endif
+
+ErrorCode read_file(string &fname, EntityHandle &seth,
+ Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems);
+void deform_func(double *xold, double *xnew);
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
+ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
+ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+string trimeshf = string(MESH_DIR) + string("/rodtri.g");
+string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
+const int SOLID_SETNO = 100, FLUID_SETNO = 200;
+
+Interface *mb;
+#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
+
+
+int main(int /*argc*/, char **/*argv*/) {
+
+ EntityHandle master, slave;
+ ErrorCode rval;
+
+ mb = new Core();
+
+ // read master/slave files and get fluid/solid material sets
+ Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
+ rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+ rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+
+ // deform the master's solid mesh, put results in a new tag
+ Tag xnew;
+ rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
+ rval = write_to_coords(solid_elems[0], xnew); RR("");
+ rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+
+ // smooth the master mesh
+ LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
+ rval = ll->perform_smooth();
+ RR("Failed in lloyd smoothing.");
+ cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+
+ rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+
+ delete ll;
+ delete mb;
+
+ return MB_SUCCESS;
+}
+
+ErrorCode write_to_coords(Range &elems, Tag tagh)
+{
+ // write the tag to coordinates
+ Range verts;
+ ErrorCode rval = mb->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get adj vertices.");
+ std::vector<double> coords(3*verts.size());
+ rval = mb->tag_get_data(tagh, verts, &coords[0]);
+ RR("Failed to get tag data.");
+ rval = mb->set_coords(verts, &coords[0]);
+ RR("Failed to set coordinates.");
+ return MB_SUCCESS;
+}
+
+void deform_func(double *xold, double *xnew)
+{
+ const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
+ // function: origin is at middle base of rod, and is .5 high
+ // top of rod is (0,.55) on left and (.2,.6) on right
+ double delx = 0.5*RODWIDTH;
+
+ double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
+ xnew[0] = xold[0] + yfrac * delx;
+ xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
+}
+
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew)
+{
+ // deform elements with an analytic function
+
+ // create the tag
+ ErrorCode rval = mb->tag_get_handle("", 3, MB_TYPE_DOUBLE, xnew, MB_TAG_CREAT|MB_TAG_DENSE);
+ RR("Failed to create xnew tag.");
+
+ // get all the vertices and coords in the fluid, set xnew to them
+ Range verts;
+ rval = mb->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ std::vector<double> coords(3*verts.size(), 0.0);
+ rval = mb->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ RR("Failed to set xnew tag on fluid verts.");
+
+ // get all the vertices and coords in the solid
+ verts.clear();
+ rval = mb->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ coords.resize(3*verts.size(), 0.0);
+ rval = mb->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ unsigned int num_verts = verts.size();
+ for (unsigned int i = 0; i < num_verts; i++)
+ deform_func(&coords[3*i], &coords[3*i]);
+
+ // set the new tag to those coords
+ rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ RR("Failed to set tag data.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode read_file(string &fname, EntityHandle &seth,
+ Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems)
+{
+ // create meshset
+ ErrorCode rval = mb->create_meshset(0, seth);
+ RR("Couldn't create master/slave set.");
+ rval = mb->load_file(fname.c_str(), &seth);
+ RR("Couldn't load master/slave mesh.");
+
+ // get material sets for solid/fluid
+ Tag tagh;
+ rval = mb->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+ const void *setno_ptr = &SOLID_SETNO;
+ rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, solids);
+ if (solids.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any solid sets.");
+
+ // get solid entities, and dimension
+ Range tmp_range;
+ for (Range::iterator rit = solids.begin(); rit != solids.end(); rit++) {
+ rval = mb->get_entities_by_handle(*rit, tmp_range, true);
+ RR("Failed to get entities in solid.");
+ }
+ int dim = mb->dimension_from_handle(*tmp_range.rbegin());
+ assert(dim > 0 && dim < 4);
+
+ solid_elems = tmp_range.subset_by_dimension(dim);
+
+ setno_ptr = &FLUID_SETNO;
+ rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, fluids);
+ if (fluids.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any fluid sets.");
+
+ for (Range::iterator rit = fluids.begin(); rit != fluids.end(); rit++) {
+ rval = mb->get_entities_by_dimension(*rit, dim, fluid_elems, true);
+ RR("Failed to get entities in fluid.");
+ }
+ if (mb->dimension_from_handle(*fluid_elems.begin()) != dim) {
+ rval = MB_FAILURE;
+ RR("Fluid and solid elements must be same dimension.");
+ }
+
+ return MB_SUCCESS;
+}
diff --git a/examples/makefile b/examples/makefile
index ff5876c..87c0132 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -4,7 +4,7 @@ include ${MOAB_DIR}/lib/iMesh-Defs.inc
.SUFFIXES: .o .cpp .F90
-# MESH_DIR is the top-level MOAB source directory, used to locate mesh files that come with MOAB source
+# MESH_DIR is the directory containing mesh files that come with MOAB source
MESH_DIR="../MeshFiles/unittest"
EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles
@@ -48,6 +48,8 @@ TestExodusII: TestExodusII.o
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
point_in_elem_search: point_in_elem_search.o
+
+DeformMeshRemap: DeformMeshRemap.o
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
clean:
https://bitbucket.org/fathomteam/moab/commits/530265f9a39c/
Changeset: 530265f9a39c
Branch: None
User: tautges
Date: 2013-09-18 01:13:16
Summary: Merge branch 'a8f966f' into deformed-mesh-remap
Affected #: 2 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index ddb3f3a..55ac3cc 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -13,6 +13,7 @@
#include "moab/Core.hpp"
#include "moab/Range.hpp"
#include "moab/LloydSmoother.hpp"
+#include "moab/ProgOptions.hpp"
#include "MBTagConventions.hpp"
#include <iostream>
@@ -32,39 +33,48 @@ ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
ErrorCode write_to_coords(Range &elems, Tag tagh);
-string trimeshf = string(MESH_DIR) + string("/rodtri.g");
-string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
const int SOLID_SETNO = 100, FLUID_SETNO = 200;
Interface *mb;
#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
-int main(int /*argc*/, char **/*argv*/) {
+int main(int argc, char **argv) {
EntityHandle master, slave;
ErrorCode rval;
+ ProgOptions po("Deformed mesh options");
+ po.addOpt<std::string> ("master,m", "Specify the master meshfile name" );
+ po.addOpt<std::string> ("slave,s", "Specify the slave meshfile name" );
+ po.parseCommandLine(argc, argv);
+ std::string foo;
+ string masterf, slavef;
+ if(!po.getOpt("master", &masterf))
+ masterf = string(MESH_DIR) + string("/rodquad.g");
+ if(!po.getOpt("slave", &slavef))
+ slavef = string(MESH_DIR) + string("/rodtri.g");
+
mb = new Core();
// read master/slave files and get fluid/solid material sets
Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
- rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
- rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+ rval = read_file(masterf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+ rval = read_file(slavef, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
// deform the master's solid mesh, put results in a new tag
Tag xnew;
rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
- rval = write_to_coords(solid_elems[0], xnew); RR("");
- rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+ if (debug) write_and_save(solid_elems[0], master, xnew, "deformed.vtk");
// smooth the master mesh
LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
rval = ll->perform_smooth();
RR("Failed in lloyd smoothing.");
cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+ if (debug) write_and_save(fluid_elems[0], master, xnew, "smoothed.vtk");
- rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+ // map new locations to slave
delete ll;
delete mb;
@@ -72,6 +82,13 @@ int main(int /*argc*/, char **/*argv*/) {
return MB_SUCCESS;
}
+ErrorCode write_and_save(Range &ents, EntithHanlde seth, Tag tagh, const char *filename)
+{
+ rval = write_to_coords(ents, tagh); RR("");
+ rval = mb->write_file("deformed.vtk", NULL, NULL, &seth, 1); RR("");
+ return rval;
+}
+
ErrorCode write_to_coords(Range &elems, Tag tagh)
{
// write the tag to coordinates
diff --git a/examples/LloydRelaxation.cpp b/examples/LloydRelaxation.cpp
index 0cf9d1e..34002cd 100644
--- a/examples/LloydRelaxation.cpp
+++ b/examples/LloydRelaxation.cpp
@@ -14,8 +14,10 @@
* in the current directory (H5M format must be used since the file is written in parallel).
*/
-#include "moab/ParallelComm.hpp"
-#include "MBParallelConventions.h"
+#ifdef USE_MPI
+# include "moab/ParallelComm.hpp"
+# include "MBParallelConventions.h"
+#endif
#include "moab/Core.hpp"
#include "moab/Skinner.hpp"
#include "moab/CN.hpp"
@@ -30,7 +32,7 @@ string test_file_name = string(MESH_DIR) + string("/surfrandomtris-4part.h5m");
#define RC if (MB_SUCCESS != rval) return rval
-ErrorCode perform_lloyd_relaxation(ParallelComm *pc, Range &verts, Range &cells, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &cells, Tag fixed,
int num_its, int report_its);
int main(int argc, char **argv)
@@ -38,7 +40,9 @@ int main(int argc, char **argv)
int num_its = 10;
int report_its = 1;
+#ifdef USE_MPI
MPI_Init(&argc, &argv);
+#endif
// need option handling here for input filename
if (argc > 1){
@@ -49,10 +53,13 @@ int main(int argc, char **argv)
// get MOAB and ParallelComm instances
Interface *mb = new Core;
if (NULL == mb) return 1;
+ int nprocs = 1;
+
+#ifdef USE_MPI
// get the ParallelComm instance
- ParallelComm* pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
- int nprocs = pcomm->size();
-
+ ParallelComm *pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
+ nprocs = pcomm->size();
+#endif
string options;
if (nprocs > 1) // if reading in parallel, need to tell it how
options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
@@ -76,34 +83,45 @@ int main(int argc, char **argv)
// ok to mark non-owned skin vertices too, I won't move those anyway
// use MOAB's skinner class to find the skin
Skinner skinner(mb);
- rval = skinner.find_skin(faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
+ rval = skinner.find_skin(0, faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
std::vector<int> fix_tag(skin_verts.size(), 1); // initialized to 1 to indicate fixed
rval = mb->tag_set_data(fixed, skin_verts, &fix_tag[0]); RC;
// now perform the Lloyd relaxation
- rval = perform_lloyd_relaxation(pcomm, verts, faces, fixed, num_its, report_its); RC;
+ rval = perform_lloyd_relaxation(mb, verts, faces, fixed, num_its, report_its); RC;
// delete fixed tag, since we created it here
rval = mb->tag_delete(fixed); RC;
// output file, using parallel write
- rval = mb->write_file("lloydfinal.h5m", NULL, "PARALLEL=WRITE_PART"); RC;
+
+#ifdef USE_MPI
+ options = "PARALLEL=WRITE_PART";
+#endif
+
+ rval = mb->write_file("lloydfinal.h5m", NULL, options.c_str()); RC;
// delete MOAB instance
delete mb;
+#ifdef USE_MPI
MPI_Finalize();
+#endif
return 0;
}
-ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &faces, Tag fixed,
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &faces, Tag fixed,
int num_its, int report_its)
{
ErrorCode rval;
- Interface *mb = pcomm->get_moab();
- int nprocs = pcomm->size();
+ int nprocs = 1;
+
+#ifdef USE_MPI
+ ParallelComm *pcomm = ParallelComm::get_pcomm(mb, 0);
+ nprocs = pcomm->size();
+#endif
// perform Lloyd relaxation:
// 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
@@ -120,8 +138,10 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// filter verts down to owned ones and get fixed tag for them
Range owned_verts, shared_owned_verts;
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
if (rval != MB_SUCCESS) return rval;
+#endif
}
else
owned_verts = verts;
@@ -133,11 +153,13 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
vcentroids.resize(3*owned_verts.size());
rval = mb->tag_get_data(centroid, owned_verts, &vcentroids[0]); RC;
+#ifdef USE_MPI
// get shared owned verts, for exchanging tags
rval = pcomm->get_shared_entities(-1, shared_owned_verts, 0, false, true); RC;
// workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
// for all shared entities
if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+#endif
// some declarations for later iterations
std::vector<double> fcentroids(3*faces.size()); // fcentroids for face centroids
@@ -195,16 +217,22 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
// 2c. exchange tags on owned verts
if (nprocs > 1) {
+#ifdef USE_MPI
rval = pcomm->exchange_tags(centroid, shared_owned_verts); RC;
+#endif
}
if (!(nit%report_its)) {
// global reduce for maximum delta, then report it
double global_max = mxdelta;
+ int myrank = 0;
+#ifdef USE_MPI
if (nprocs > 1)
MPI_Reduce(&mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm());
- if (1 == nprocs || !pcomm->rank())
+ myrank = pcomm->rank();
+#endif
+ if (1 == nprocs || !myrank)
cout << "Max delta = " << global_max << endl;
}
}
https://bitbucket.org/fathomteam/moab/commits/e2e1ffac2d54/
Changeset: e2e1ffac2d54
Branch: deformed-mesh-remap
User: tautges
Date: 2013-10-02 18:38:22
Summary: Make a DeformMeshRemap class and implement through that.
Affected #: 1 file
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 55ac3cc..fac83ae 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -17,6 +17,7 @@
#include "MBTagConventions.hpp"
#include <iostream>
+#include <set>
#include <assert.h>
using namespace moab;
@@ -37,11 +38,212 @@ const int SOLID_SETNO = 100, FLUID_SETNO = 200;
Interface *mb;
#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
-
+
+const bool debug = true;
+
+class DeformMeshRemap
+{
+public:
+
+ //! enumerator for solid/fluid, master/slave
+ enum {MASTER=0, SLAVE, SOLID, FLUID};
+
+ //! constructor
+ DeformMeshRemap(Interface *impl) : mbImpl(impl), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew") {}
+
+ //! destructor
+ ~DeformMeshRemap();
+
+ //! execute the deformed mesh process
+ ErrorCode execute();
+
+ //! add a set number
+ ErrorCode add_set_no(int fluid_or_solid, int set_no);
+
+ //! remove a set number
+ ErrorCode remove_set_no(int fluid_or_solid, int set_no);
+
+ //! get the set numbers
+ ErrorCode get_set_nos(int fluid_or_solid, std::set<int> &set_nos) const;
+
+ //! get the xNew tag handle
+ inline Tag x_new() const {return xNew;}
+
+ //! get the tag name
+ std::string x_new_name() const {return xNewName;}
+
+ //! set the tag name
+ void x_new_name(const std::string &name) {xNewName = name;}
+
+ //! get/set the file name
+ std::string get_file_name(int m_or_s) const;
+
+ //! get/set the file name
+ void set_file_name(int m_or_s, const std::string &name);
+
+private:
+ //! apply a known deformation to the solid elements, putting the results in the xNew tag; also
+ //! write current coordinates to the xNew tag for fluid elements
+ ErrorCode deform_master(Range &fluid_elems, Range &solid_elems);
+
+ //! read a file and establish proper ranges
+ ErrorCode read_file(int m_or_s, string &fname, EntityHandle &seth);
+
+ //! write the input tag to the coordinates for the vertices in the input elems
+ ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+ //! write the tag to the vertices, then save to the specified file
+ ErrorCode write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename);
+
+ //! moab interface
+ Interface *mbImpl;
+
+ //! material set numbers for fluid materials
+ std::set<int> fluidSetNos;
+
+ //! material set numbers for solid materials
+ std::set<int> solidSetNos;
+
+ //! sets defining master/slave meshes
+ EntityHandle masterSet, slaveSet;
+
+ //! sets in master/slave meshes
+ Range fluidSets[2], solidSets[2];
+
+ //! elements in master/slave meshes
+ Range fluidElems[2], solidElems[2];
+
+ //! filenames for master/slave meshes
+ std::string masterFileName, slaveFileName;
+
+ //! tag used for new positions
+ Tag xNew;
+
+ //! tag name used for new positions
+ std::string xNewName;
+};
+
+ //! add a set number
+inline ErrorCode DeformMeshRemap::add_set_no(int f_or_s, int set_no)
+{
+ std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+
+ this_set->insert(set_no);
+
+ return MB_SUCCESS;
+}
+
+ //! remove a set number
+inline ErrorCode DeformMeshRemap::remove_set_no(int f_or_s, int set_no)
+{
+ std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+ std::set<int>::iterator sit = this_set->find(set_no);
+ if (sit != this_set->end()) {
+ this_set->erase(*sit);
+ return MB_SUCCESS;
+ }
+
+ return MB_FAILURE;
+}
+
+ //! get the set numbers
+inline ErrorCode DeformMeshRemap::get_set_nos(int f_or_s, std::set<int> &set_nos) const
+{
+ const std::set<int> *this_set;
+ switch (f_or_s) {
+ case FLUID:
+ this_set = &fluidSetNos; break;
+ case SOLID:
+ this_set = &solidSetNos; break;
+ default:
+ assert(false && "f_or_s should be FLUID or SOLID.");
+ return MB_FAILURE;
+ }
+
+ set_nos = *this_set;
+
+ return MB_SUCCESS;
+}
+
+ErrorCode DeformMeshRemap::execute()
+{
+ // read master/slave files and get fluid/solid material sets
+ ErrorCode rval = read_file(MASTER, masterFileName, masterSet);
+ if (MB_SUCCESS != rval) return rval;
+
+ rval = read_file(SLAVE, slaveFileName, slaveSet);
+ if (MB_SUCCESS != rval) return rval;
+
+ // deform the master's solid mesh, put results in a new tag
+ rval = deform_master(fluidElems[MASTER], solidElems[MASTER]); RR("");
+ if (debug) write_and_save(solidElems[MASTER], masterSet, xNew, "deformed.vtk");
+
+ // smooth the master mesh
+ LloydSmoother *ll = new LloydSmoother(mbImpl, NULL, fluidElems[MASTER], xNew);
+ rval = ll->perform_smooth();
+ RR("Failed in lloyd smoothing.");
+
+ cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+ if (debug) write_and_save(fluidElems[MASTER], masterSet, xNew, "smoothed.vtk");
+
+ // map new locations to slave
+
+ delete ll;
+
+ return MB_SUCCESS;
+}
+
+std::string DeformMeshRemap::get_file_name(int m_or_s) const
+{
+ switch (m_or_s) {
+ case MASTER:
+ return masterFileName;
+ case SLAVE:
+ return slaveFileName;
+ default:
+ assert(false && "m_or_s should be MASTER or SLAVE.");
+ return std::string();
+ }
+}
+
+void DeformMeshRemap::set_file_name(int m_or_s, const std::string &name)
+{
+ switch (m_or_s) {
+ case MASTER:
+ masterFileName = name; break;
+ case SLAVE:
+ slaveFileName = name; break;
+ default:
+ assert(false && "m_or_s should be MASTER or SLAVE.");
+ }
+}
+
+DeformMeshRemap::~DeformMeshRemap()
+{
+ // delete the tag
+ mbImpl->tag_delete(xNew);
+}
int main(int argc, char **argv) {
- EntityHandle master, slave;
ErrorCode rval;
ProgOptions po("Deformed mesh options");
@@ -56,49 +258,34 @@ int main(int argc, char **argv) {
slavef = string(MESH_DIR) + string("/rodtri.g");
mb = new Core();
-
- // read master/slave files and get fluid/solid material sets
- Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
- rval = read_file(masterf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
- rval = read_file(slavef, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
- // deform the master's solid mesh, put results in a new tag
- Tag xnew;
- rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
- if (debug) write_and_save(solid_elems[0], master, xnew, "deformed.vtk");
+ DeformMeshRemap dfr(mb);
+ dfr.set_file_name(DeformMeshRemap::MASTER, masterf);
+ dfr.set_file_name(DeformMeshRemap::SLAVE, slavef);
+ rval = dfr.add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+ rval = dfr.add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add fluid set no.");
- // smooth the master mesh
- LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
- rval = ll->perform_smooth();
- RR("Failed in lloyd smoothing.");
- cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
- if (debug) write_and_save(fluid_elems[0], master, xnew, "smoothed.vtk");
-
- // map new locations to slave
-
- delete ll;
- delete mb;
-
- return MB_SUCCESS;
+ rval = dfr.execute();
+ return rval;
}
-ErrorCode write_and_save(Range &ents, EntithHanlde seth, Tag tagh, const char *filename)
+ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename)
{
- rval = write_to_coords(ents, tagh); RR("");
- rval = mb->write_file("deformed.vtk", NULL, NULL, &seth, 1); RR("");
+ ErrorCode rval = write_to_coords(ents, tagh); RR("");
+ rval = mbImpl->write_file(filename, NULL, NULL, &seth, 1); RR("");
return rval;
}
-ErrorCode write_to_coords(Range &elems, Tag tagh)
+ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh)
{
// write the tag to coordinates
Range verts;
- ErrorCode rval = mb->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+ ErrorCode rval = mbImpl->get_adjacencies(elems, 0, false, verts, Interface::UNION);
RR("Failed to get adj vertices.");
std::vector<double> coords(3*verts.size());
- rval = mb->tag_get_data(tagh, verts, &coords[0]);
+ rval = mbImpl->tag_get_data(tagh, verts, &coords[0]);
RR("Failed to get tag data.");
- rval = mb->set_coords(verts, &coords[0]);
+ rval = mbImpl->set_coords(verts, &coords[0]);
RR("Failed to set coordinates.");
return MB_SUCCESS;
}
@@ -115,83 +302,92 @@ void deform_func(double *xold, double *xnew)
xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
}
-ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew)
+ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems)
{
// deform elements with an analytic function
// create the tag
- ErrorCode rval = mb->tag_get_handle("", 3, MB_TYPE_DOUBLE, xnew, MB_TAG_CREAT|MB_TAG_DENSE);
+ ErrorCode rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, xNew, MB_TAG_CREAT|MB_TAG_DENSE);
RR("Failed to create xnew tag.");
// get all the vertices and coords in the fluid, set xnew to them
Range verts;
- rval = mb->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+ rval = mbImpl->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
RR("Failed to get vertices.");
std::vector<double> coords(3*verts.size(), 0.0);
- rval = mb->get_coords(verts, &coords[0]);
+ rval = mbImpl->get_coords(verts, &coords[0]);
RR("Failed to get vertex coords.");
- rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
RR("Failed to set xnew tag on fluid verts.");
// get all the vertices and coords in the solid
verts.clear();
- rval = mb->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+ rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
RR("Failed to get vertices.");
coords.resize(3*verts.size(), 0.0);
- rval = mb->get_coords(verts, &coords[0]);
+ rval = mbImpl->get_coords(verts, &coords[0]);
RR("Failed to get vertex coords.");
unsigned int num_verts = verts.size();
for (unsigned int i = 0; i < num_verts; i++)
deform_func(&coords[3*i], &coords[3*i]);
// set the new tag to those coords
- rval = mb->tag_set_data(xnew, verts, &coords[0]);
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
RR("Failed to set tag data.");
return MB_SUCCESS;
}
-ErrorCode read_file(string &fname, EntityHandle &seth,
- Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems)
+ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &seth)
{
// create meshset
- ErrorCode rval = mb->create_meshset(0, seth);
+ ErrorCode rval = mbImpl->create_meshset(0, seth);
RR("Couldn't create master/slave set.");
- rval = mb->load_file(fname.c_str(), &seth);
+ rval = mbImpl->load_file(fname.c_str(), &seth);
RR("Couldn't load master/slave mesh.");
// get material sets for solid/fluid
Tag tagh;
- rval = mb->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
- const void *setno_ptr = &SOLID_SETNO;
- rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, solids);
- if (solids.empty()) rval = MB_FAILURE;
- RR("Couldn't get any solid sets.");
+ rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+ for (std::set<int>::iterator sit = solidSetNos.begin(); sit != solidSetNos.end(); sit++) {
+ Range sets;
+ int set_no = *sit;
+ const void *setno_ptr = &set_no;
+ rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+ if (sets.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any solid sets.");
+ solidSets[m_or_s].merge(sets);
+ }
// get solid entities, and dimension
Range tmp_range;
- for (Range::iterator rit = solids.begin(); rit != solids.end(); rit++) {
- rval = mb->get_entities_by_handle(*rit, tmp_range, true);
+ for (Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); rit++) {
+ rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
RR("Failed to get entities in solid.");
}
- int dim = mb->dimension_from_handle(*tmp_range.rbegin());
+ int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
assert(dim > 0 && dim < 4);
- solid_elems = tmp_range.subset_by_dimension(dim);
+ solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
- setno_ptr = &FLUID_SETNO;
- rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, fluids);
- if (fluids.empty()) rval = MB_FAILURE;
- RR("Couldn't get any fluid sets.");
-
- for (Range::iterator rit = fluids.begin(); rit != fluids.end(); rit++) {
- rval = mb->get_entities_by_dimension(*rit, dim, fluid_elems, true);
- RR("Failed to get entities in fluid.");
+ for (std::set<int>::iterator sit = fluidSetNos.begin(); sit != fluidSetNos.end(); sit++) {
+ Range sets;
+ int set_no = *sit;
+ const void *setno_ptr = &set_no;
+ rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+ if (sets.empty()) rval = MB_FAILURE;
+ RR("Couldn't get any fluid sets.");
+ fluidSets[m_or_s].merge(sets);
}
- if (mb->dimension_from_handle(*fluid_elems.begin()) != dim) {
- rval = MB_FAILURE;
- RR("Fluid and solid elements must be same dimension.");
+
+ // get fluid entities, and dimension
+ tmp_range.clear();
+ for (Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); rit++) {
+ rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+ RR("Failed to get entities in fluid.");
}
+ fluidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+
return MB_SUCCESS;
}
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.
More information about the moab-dev
mailing list