[MOAB-dev] commit/MOAB: tautges: Merge branch 'deformed-mesh-remap' of git at bitbucket.org:fathomteam/moab.git
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Jan 7 14:05:17 CST 2014
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/37765e24606e/
Changeset: 37765e24606e
Branch: master
User: tautges
Date: 2014-01-07 21:04:26
Summary: Merge branch 'deformed-mesh-remap' of git at bitbucket.org:fathomteam/moab.git
Affected #: 47 files
diff --git a/.gitignore b/.gitignore
index e824a6b..d842c2b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -165,6 +165,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/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
new file mode 100644
index 0000000..5442d1d
--- /dev/null
+++ b/examples/DeformMeshRemap.cpp
@@ -0,0 +1,488 @@
+/** @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 "moab/ProgOptions.hpp"
+#include "moab/BoundBox.hpp"
+#include "moab/SpatialLocator.hpp"
+#include "MBTagConventions.hpp"
+#include "DataCoupler.hpp"
+
+#ifdef USE_MPI
+# include "moab/ParallelComm.hpp"
+#endif
+
+#include <iostream>
+#include <set>
+#include <sstream>
+#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(const BoundBox &bbox, 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);
+
+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
+ //! if master is NULL, the MOAB part is run in serial;
+ //! if slave is NULL but the master isn't, the slave is copied from the master
+ //! Create communicators using moab::ParallelComm::get_pcomm
+ DeformMeshRemap(Interface *impl, ParallelComm *master = NULL, ParallelComm *slave = NULL);
+
+ //! 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, const char *tag_name = NULL);
+
+ //! 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;
+
+#ifdef USE_MPI
+ //! ParallelComm for master, slave meshes
+ ParallelComm *pcMaster, *pcSlave;
+#endif
+
+ //! 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;
+
+ Range src_elems = solidElems[MASTER];
+ src_elems.merge(fluidElems[MASTER]);
+ // locate slave vertices in master, orig coords; do this with a data coupler, so you can
+ // later interpolate
+ Range tgt_verts, tmp_range = solidElems[SLAVE];
+ tmp_range.merge(fluidElems[SLAVE]);
+ rval = mbImpl->get_adjacencies(tmp_range, 0, false, tgt_verts, Interface::UNION);
+ RR("Failed to get target verts.");
+
+
+ // initialize data coupler on source elements
+ DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
+
+ // locate slave vertices, caching results in dc
+ rval = dc_master.locate_points(tgt_verts); RR("Point location of tgt verts failed.");
+ int num_located = dc_master.spatial_locator()->local_num_located();
+ if (num_located != (int)tgt_verts.size()) {
+ rval = MB_FAILURE;
+ std::cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located." << std::endl;
+ return rval;
+ }
+
+ // deform the master's solid mesh, put results in a new tag
+ rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
+
+ { // to isolate the lloyd smoother & delete when done
+
+ // smooth the master mesh
+ LloydSmoother ll(mbImpl, NULL, fluidElems[MASTER], xNew);
+ rval = ll.perform_smooth();
+ RR("Failed in lloyd smoothing.");
+ cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
+ }
+
+ // map new locations to slave
+ // interpolate xNew to slave points
+ rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
+
+ // transfer xNew to coords, for master and slave
+ rval = write_to_coords(fluidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
+ rval = write_to_coords(tgt_verts, xNew); RR("Failed writing tag to slave verts.");
+
+ if (debug) {
+ std::string str;
+#ifdef USE_MPI
+ if (pcMaster && pcMaster->size() > 1)
+ str = "PARALLEL=WRITE_PART";
+#endif
+ rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
+#ifdef USE_MPI
+ str.clear();
+ if (pcSlave && pcSlave->size() > 1)
+ str = "PARALLEL=WRITE_PART";
+#endif
+ rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
+ }
+
+ 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(Interface *impl, ParallelComm *master, ParallelComm *slave)
+ : mbImpl(impl), pcMaster(master), pcSlave(slave), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew")
+{
+ if (!pcSlave && pcMaster)
+ pcSlave = pcMaster;
+}
+
+DeformMeshRemap::~DeformMeshRemap()
+{
+ // delete the tag
+ mbImpl->tag_delete(xNew);
+}
+
+int main(int argc, char **argv) {
+
+ 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();
+
+ DeformMeshRemap *dfr;
+#ifdef USE_MPI
+ ParallelComm *pc = new ParallelComm(mb, MPI_COMM_WORLD);
+ dfr = new DeformMeshRemap(mb, pc);
+#else
+ dfr = new DeformMeshRemap(mb);
+#endif
+ 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.");
+
+ rval = dfr->execute();
+
+ delete dfr;
+ delete mb;
+
+ return rval;
+}
+
+ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename)
+{
+ ErrorCode rval = write_to_coords(ents, tagh); RR("");
+ rval = mbImpl->write_file(filename, NULL, NULL, &seth, 1); RR("");
+ return rval;
+}
+
+ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh)
+{
+ // write the tag to coordinates
+ Range verts;
+ 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 = mbImpl->tag_get_data(tagh, verts, &coords[0]);
+ RR("Failed to get tag data.");
+ rval = mbImpl->set_coords(verts, &coords[0]);
+ RR("Failed to set coordinates.");
+ return MB_SUCCESS;
+}
+
+void deform_func(const BoundBox &bbox, double *xold, double *xnew)
+{
+/* Deformation function based on max delx and dely at top of rod
+ 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;
+*/
+
+/* Deformation function based on fraction of bounding box dimension in each direction */
+ double frac = 0.01; // taken from approximate relative deformation from LLNL Diablo of XX09 assys
+ CartVect *xo = reinterpret_cast<CartVect*>(xold), *xn = reinterpret_cast<CartVect*>(xnew);
+ CartVect disp = frac * (*xo - bbox.bMin);
+ *xn = *xo + disp;
+}
+
+ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name)
+{
+ // deform elements with an analytic function
+
+ // create the tag
+ ErrorCode rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 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 = 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 = mbImpl->get_coords(verts, &coords[0]);
+ RR("Failed to get vertex coords.");
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+ RR("Failed to set xnew tag on fluid verts.");
+
+ // get the bounding box of the solid mesh
+ BoundBox bbox;
+ bbox.update(*mbImpl, solid_elems);
+
+ // get all the vertices and coords in the solid
+ verts.clear();
+ rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+ RR("Failed to get vertices.");
+ coords.resize(3*verts.size(), 0.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(bbox, &coords[3*i], &coords[3*i]);
+
+ // set the new tag to those coords
+ rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+ RR("Failed to set tag data.");
+
+ return MB_SUCCESS;
+}
+
+ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &seth)
+{
+ // create meshset
+ ErrorCode rval = mbImpl->create_meshset(0, seth);
+ RR("Couldn't create master/slave set.");
+ ostringstream ostr;
+#ifdef USE_MPI
+ ParallelComm *pc = (m_or_s == MASTER ? pcMaster : pcSlave);
+ if (pc && pc->size() > 1) {
+ if (debug) ostr << "DEBUG_IO=1;CPUTIME;";
+ ostr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
+ << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
+ }
+#endif
+ rval = mbImpl->load_file(fname.c_str(), &seth, ostr.str().c_str());
+ RR("Couldn't load master/slave mesh.");
+
+ // get material sets for solid/fluid
+ Tag tagh;
+ 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 = 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 = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+ assert(dim > 0 && dim < 4);
+
+ solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+
+ 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);
+ }
+
+ // 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;
+}
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;
}
}
diff --git a/examples/makefile b/examples/makefile
index 5ae5af0..a6e3d1f 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
@@ -14,45 +14,48 @@ EXOIIEXAMPLES = TestExodusII
default: ${EXAMPLES}
-HelloMOAB : HelloMOAB.o
+HelloMOAB : HelloMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-GetEntities: GetEntities.o
+GetEntities: GetEntities.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-SetsNTags: SetsNTags.o
+SetsNTags: SetsNTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-LloydRelaxation: LloydRelaxation.o
+LloydRelaxation: LloydRelaxation.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-StructuredMeshSimple : StructuredMeshSimple.o
+StructuredMeshSimple : StructuredMeshSimple.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessWithHoles: DirectAccessWithHoles.o
+DirectAccessWithHoles: DirectAccessWithHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHoles: DirectAccessNoHoles.o
+DirectAccessNoHoles: DirectAccessNoHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHolesF90: DirectAccessNoHolesF90.o
+DirectAccessNoHolesF90: DirectAccessNoHolesF90.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${IMESH_LIBS}
-ReduceExchangeTags : ReduceExchangeTags.o
+ReduceExchangeTags : ReduceExchangeTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-HelloParMOAB: HelloParMOAB.o
+HelloParMOAB: HelloParMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-TestExodusII: TestExodusII.o
+TestExodusII: TestExodusII.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-point_in_elem_search: point_in_elem_search.o
+point_in_elem_search: point_in_elem_search.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
PushParMeshIntoMoabF90: PushParMeshIntoMoabF90.o
${MOAB_CXX} -o $@ $< ${IMESH_LIBS} -lgfortran -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl
+DeformMeshRemap: DeformMeshRemap.o ${MOAB_LIBDIR}/libMOAB.la
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} -lmbcoupler ${MOAB_LIBS_LINK}
+
clean:
rm -rf *.o ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES}
diff --git a/src/AdaptiveKDTree.cpp b/src/AdaptiveKDTree.cpp
index f66acbc..1324c30 100644
--- a/src/AdaptiveKDTree.cpp
+++ b/src/AdaptiveKDTree.cpp
@@ -1303,7 +1303,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -1322,7 +1323,7 @@ namespace moab {
treeStats.nodesVisited++;
ErrorCode rval = get_bounding_box(box, &node);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) return MB_SUCCESS;
+ if (!box.contains_point(point, iter_tol)) return MB_SUCCESS;
rval = moab()->get_child_meshsets( node, children );
if (MB_SUCCESS != rval)
@@ -1346,7 +1347,7 @@ namespace moab {
treeStats.leavesVisited++;
if (myEval && params) {
- rval = myEval->find_containing_entity(node, point, tol,
+ rval = myEval->find_containing_entity(node, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
}
@@ -1358,7 +1359,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
AdaptiveKDTreeIter& leaf_it,
- double tol,
+ const double iter_tol,
+ const double /*inside_tol*/,
bool *multiple_leaves,
EntityHandle *start_node)
{
@@ -1372,7 +1374,7 @@ namespace moab {
leaf_it.mBox[1] = boundBox.bMax;
// test that point is inside tree
- if (!boundBox.contains_point(point, tol)) {
+ if (!boundBox.contains_point(point, iter_tol)) {
treeStats.nodesVisited++;
return MB_ENTITY_NOT_FOUND;
}
@@ -1423,7 +1425,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -1445,7 +1448,7 @@ namespace moab {
// (zero if inside box)
BoundBox box;
rval = get_bounding_box(box);
- if (MB_SUCCESS == rval && !box.contains_point(from_point, tol)) {
+ if (MB_SUCCESS == rval && !box.contains_point(from_point, iter_tol)) {
treeStats.nodesVisited++;
return MB_SUCCESS;
}
@@ -1484,7 +1487,7 @@ namespace moab {
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(node.handle, from_point, tol,
+ rval = myEval->find_containing_entity(node.handle, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
@@ -1802,7 +1805,7 @@ namespace moab {
// the same distance from the input point as the current closest
// point is.
CartVect diff = closest_pt - from;
- rval = distance_search(from_coords, sqrt(diff%diff), leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(from_coords, sqrt(diff%diff), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// Check any close leaves to see if they contain triangles that
@@ -1834,7 +1837,7 @@ namespace moab {
// get leaves of tree that intersect sphere
assert(tree_root);
- rval = distance_search(center, radius, leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(center, radius, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// search each leaf for triangles intersecting sphere
diff --git a/src/BVHTree.cpp b/src/BVHTree.cpp
index 8903070..b65eae4 100644
--- a/src/BVHTree.cpp
+++ b/src/BVHTree.cpp
@@ -435,7 +435,8 @@ namespace moab
ErrorCode BVHTree::find_point(const std::vector<double> &point,
const unsigned int &index,
- const double tol,
+ const double iter_tol,
+ const double inside_tol,
std::pair<EntityHandle, CartVect> &result)
{
if (index == 0) treeStats.numTraversals++;
@@ -453,7 +454,7 @@ namespace moab
for(Range::iterator i = entities.begin(); i != entities.end(); i++) {
treeStats.leafObjectTests++;
myEval->set_ent_handle(*i);
- myEval->reverse_eval(&point[0], tol, params.array(), &is_inside);
+ myEval->reverse_eval(&point[0], iter_tol, inside_tol, params.array(), &is_inside);
if (is_inside) {
result.first = *i;
result.second = params;
@@ -469,11 +470,11 @@ namespace moab
rval = mbImpl->get_child_meshsets(startSetHandle+index, children);
if (MB_SUCCESS != rval || children.size() != 2) return rval;
- if((node.Lmax+tol) < (node.Rmin-tol)){
- if(point[node.dim] <= (node.Lmax + tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
- else if(point[ node.dim] >= (node.Rmin - tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ if((node.Lmax+iter_tol) < (node.Rmin-iter_tol)){
+ if(point[node.dim] <= (node.Lmax + iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
+ else if(point[ node.dim] >= (node.Rmin - iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
result = std::make_pair(0, CartVect(&point[0])); //point lies in empty space.
return MB_SUCCESS;
}
@@ -482,11 +483,11 @@ namespace moab
//left of Rmin, you must be on the left
//we can't be sure about the boundaries since the boxes overlap
//this was a typo in the paper which caused pain.
- if(point[node.dim] < (node.Rmin - tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
+ if(point[node.dim] < (node.Rmin - iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
//if you are on the right Lmax, you must be on the right
- else if(point[ node.dim] > (node.Lmax+tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ else if(point[ node.dim] > (node.Lmax+iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
/* pg5 of paper
* However, instead of always traversing either subtree
@@ -501,22 +502,22 @@ namespace moab
//branch predicted..
//bool dir = (point[ node.dim] - node.Rmin) <=
// (node.Lmax - point[ node.dim]);
- find_point(point, children[0]-startSetHandle, tol, result);
+ find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
if(result.first == 0){
- return find_point(point, children[1]-startSetHandle, tol, result);
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
}
return MB_SUCCESS;
}
- EntityHandle BVHTree::bruteforce_find(const double *point, const double tol) {
+ EntityHandle BVHTree::bruteforce_find(const double *point, const double iter_tol, const double inside_tol) {
treeStats.numTraversals++;
CartVect params;
for(unsigned int i = 0; i < myTree.size(); i++) {
- if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, tol)) continue;
+ if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, iter_tol)) continue;
if (myEval) {
EntityHandle entity = 0;
treeStats.leavesVisited++;
- ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, tol,
+ ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, iter_tol, inside_tol,
entity, params.array(), &treeStats.leafObjectTests);
if (entity) return entity;
else if (MB_SUCCESS != rval) return 0;
@@ -542,7 +543,8 @@ namespace moab
ErrorCode BVHTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -572,7 +574,7 @@ namespace moab
// test box of this node
ErrorCode rval = get_bounding_box(box, &this_set);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) continue;
+ if (!box.contains_point(point, iter_tol)) continue;
// else if not a leaf, test children & put on list
else if (myTree[ind].dim != 3) {
@@ -581,7 +583,7 @@ namespace moab
continue;
}
else if (myTree[ind].dim == 3 && myEval && params) {
- rval = myEval->find_containing_entity(startSetHandle+ind, point, tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (leaf_out || MB_SUCCESS != rval) return rval;
}
@@ -599,7 +601,8 @@ namespace moab
ErrorCode BVHTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double params_tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -653,7 +656,7 @@ namespace moab
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(startSetHandle+ind, from_point, params_tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
new file mode 100644
index 0000000..b1614d2
--- /dev/null
+++ b/src/LloydSmoother.cpp
@@ -0,0 +1,232 @@
+#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"
+#include "MBParallelConventions.h"
+#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 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);
+ 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());
+ 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);
+ 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 && 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 && myPcomm->size() > 1) {
+ rval = myPcomm->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 && myPcomm->size() > 1)
+ MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
+ if (!myPcomm || !myPcomm->rank())
+#endif
+ std::cout << "Max residual = " << global_max << std::endl;
+ }
+
+ } // end while
+
+ // write the tag back onto vertex coordinates
+ 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;
+}
+
+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/LocalDiscretization/ElemEvaluator.cpp b/src/LocalDiscretization/ElemEvaluator.cpp
index b132da9..4ca3943 100644
--- a/src/LocalDiscretization/ElemEvaluator.cpp
+++ b/src/LocalDiscretization/ElemEvaluator.cpp
@@ -6,6 +6,7 @@
// need to include eval set types here to support get_eval_set; alternative would be to have some
// type of registration, but we'd still need static registration for the built-in types
+#include "moab/LinearTri.hpp"
#include "moab/LinearQuad.hpp"
#include "moab/LinearTet.hpp"
#include "moab/LinearHex.hpp"
@@ -16,12 +17,12 @@
namespace moab {
ErrorCode EvalSet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -30,41 +31,48 @@ namespace moab {
CartVect new_pos;
// evaluate that first guess to get a new position
- ErrorCode rval = (*eval)(cvparams->array(), verts, ndim, ndim, work, new_pos.array());
+ ErrorCode rval = (*eval)(cvparams->array(), verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
// residual is diff between old and new pos; need to minimize that
CartVect res = new_pos - *cvposn;
Matrix3 J;
+ bool dum, *tmp_inside = (inside ? inside : &dum);
int iters=0;
// while |res| larger than tol
while (res % res > error_tol_sqr) {
if(++iters>10) {
- if (inside) {
- // if we haven't converged but we're outside, that's defined as success
- *inside = (*inside_f)(params, ndim, tol);
- if (!(*inside)) return MB_SUCCESS;
- }
- return MB_FAILURE;
+ // if we haven't converged but we're outside, that's defined as success
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
}
// get jacobian at current params
rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
double det = J.determinant();
- assert(det > std::numeric_limits<double>::epsilon());
+ if (det < std::numeric_limits<double>::epsilon()) {
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
+ }
// new params tries to eliminate residual
*cvparams -= J.inverse(1.0/det) * res;
// get the new forward-evaluated position, and its difference from the target pt
- rval = (*eval)(params, verts, ndim, ndim, work, new_pos.array());
+ rval = (*eval)(params, verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
res = new_pos - *cvposn;
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
@@ -85,6 +93,7 @@ namespace moab {
case MBEDGE:
break;
case MBTRI:
+ if (LinearTri::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
break;
case MBQUAD:
if (LinearQuad::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
@@ -105,9 +114,9 @@ namespace moab {
return MB_NOT_IMPLEMENTED;
}
- ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, double tol,
- EntityHandle &containing_ent, double *params,
- unsigned int *num_evals)
+ ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, const double iter_tol,
+ const double inside_tol, EntityHandle &containing_ent,
+ double *params, unsigned int *num_evals)
{
bool is_inside;
ErrorCode rval = MB_SUCCESS;
@@ -116,7 +125,7 @@ namespace moab {
for(i = entities.begin(); i != entities.end(); i++) {
nevals++;
set_ent_handle(*i);
- rval = reverse_eval(point, tol, params, &is_inside);
+ rval = reverse_eval(point, iter_tol, inside_tol, params, &is_inside);
if (MB_SUCCESS != rval) return rval;
if (is_inside) break;
}
diff --git a/src/LocalDiscretization/LinearHex.cpp b/src/LocalDiscretization/LinearHex.cpp
index 14de2d9..5e730af 100644
--- a/src/LocalDiscretization/LinearHex.cpp
+++ b/src/LocalDiscretization/LinearHex.cpp
@@ -96,10 +96,12 @@ namespace moab
ErrorCode LinearHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
}
bool LinearHex::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/LinearQuad.cpp b/src/LocalDiscretization/LinearQuad.cpp
index ddd4918..63c9959 100644
--- a/src/LocalDiscretization/LinearQuad.cpp
+++ b/src/LocalDiscretization/LinearQuad.cpp
@@ -16,23 +16,23 @@ namespace moab
*/
const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
- ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int ndim,
+ ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int /*ndim*/,
double *, double *result)
{
Matrix3 *J = reinterpret_cast<Matrix3*>(result);
*J = Matrix3(0.0);
for (unsigned i = 0; i < 4; ++i) {
- const double params_p = 1 + params[0]*corner[i][0];
+ const double xi_p = 1 + params[0]*corner[i][0];
const double eta_p = 1 + params[1]*corner[i][1];
- const double dNi_dparams = corner[i][0] * eta_p;
- const double dNi_deta = corner[i][1] * params_p;
- (*J)(0,0) += dNi_dparams * verts[i*ndim+0];
- (*J)(1,0) += dNi_dparams * verts[i*ndim+1];
- (*J)(0,1) += dNi_deta * verts[i*ndim+0];
- (*J)(1,1) += dNi_deta * verts[i*ndim+1];
+ const double dNi_dxi = corner[i][0] * eta_p;
+ const double dNi_deta = corner[i][1] * xi_p;
+ (*J)(0,0) += dNi_dxi * verts[i*3+0];
+ (*J)(1,0) += dNi_dxi * verts[i*3+1];
+ (*J)(0,1) += dNi_deta * verts[i*3+0];
+ (*J)(1,1) += dNi_deta * verts[i*3+1];
}
- (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
(*J) *= 0.25;
+ (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
return MB_SUCCESS;
}// LinearQuad::jacobian()
@@ -62,7 +62,7 @@ namespace moab
for(unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2) {
x[1] = LinearQuad::gauss[j2][1];
double w2 = LinearQuad::gauss[j2][0];
- rval = evalFcn(x.array(),field, ndim, num_tuples, NULL, tmp_result);
+ rval = evalFcn(x.array(), field, ndim, num_tuples, NULL, tmp_result);
if (MB_SUCCESS != rval) return rval;
rval = jacobianFcn(x.array(), verts, nverts, ndim, work, J[0]);
if (MB_SUCCESS != rval) return rval;
@@ -75,9 +75,11 @@ namespace moab
ErrorCode LinearQuad::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
}
bool LinearQuad::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/LinearTet.cpp b/src/LocalDiscretization/LinearTet.cpp
index 705dbee..4965cfc 100644
--- a/src/LocalDiscretization/LinearTet.cpp
+++ b/src/LocalDiscretization/LinearTet.cpp
@@ -74,27 +74,29 @@ namespace moab
ErrorCode LinearTet::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
+ work, params, is_inside);
}
bool LinearTet::insideFcn(const double *params, const int , const double tol)
{
return (params[0] >= -1.0-tol && params[1] >= -1.0-tol && params[2] >= -1.0-tol &&
- params[0] + params[1] + params[2] <= -1.0+tol);
+ params[0] + params[1] + params[2] <= 1.0+tol);
}
ErrorCode LinearTet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -138,7 +140,7 @@ namespace moab
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
diff --git a/src/LocalDiscretization/LinearTri.cpp b/src/LocalDiscretization/LinearTri.cpp
new file mode 100644
index 0000000..8fe9d23
--- /dev/null
+++ b/src/LocalDiscretization/LinearTri.cpp
@@ -0,0 +1,145 @@
+#include "moab/LinearTri.hpp"
+#include "moab/Forward.hpp"
+#include <algorithm>
+
+namespace moab
+{
+
+ const double LinearTri::corner[3][2] = { {0,0},
+ {1,0},
+ {0,1}};
+
+ ErrorCode LinearTri::initFcn(const double *verts, const int /*nverts*/, double *&work) {
+ // allocate work array as:
+ // work[0..8] = T
+ // work[9..17] = Tinv
+ // work[18] = detT
+ // work[19] = detTinv
+ assert(!work && verts);
+ work = new double[20];
+ Matrix3 *T = reinterpret_cast<Matrix3*>(work),
+ *Tinv = reinterpret_cast<Matrix3*>(work+9);
+ double *detT = work+18, *detTinv = work+19;
+
+ *T = Matrix3(verts[1*3+0]-verts[0*3+0],verts[2*3+0]-verts[0*3+0],0.0,
+ verts[1*3+1]-verts[0*3+1],verts[2*3+1]-verts[0*3+1],0.0,
+ verts[1*3+2]-verts[0*3+2],verts[2*3+2]-verts[0*3+2],1.0);
+ *T *= 0.5;
+ (*T)(2,2) = 1.0;
+
+ *Tinv = T->inverse();
+ *detT = T->determinant();
+ *detTinv = (0.0 == *detT ? HUGE : 1.0 / *detT);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples,
+ double */*work*/, double *result) {
+ assert(params && field && num_tuples > 0);
+ // convert to [0,1]
+ double p1 = 0.5 * (1.0 + params[0]),
+ p2 = 0.5 * (1.0 + params[1]),
+ p0 = 1.0 - p1 - p2;
+
+ for (int j = 0; j < num_tuples; j++)
+ result[j] = p0 * field[0*num_tuples+j] + p1 * field[1*num_tuples+j] + p2 * field[2*num_tuples+j];
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::integrateFcn(const double *field, const double */*verts*/, const int /*nverts*/, const int /*ndim*/, const int num_tuples,
+ double *work, double *result)
+ {
+ assert(field && num_tuples > 0);
+ double tmp = work[18];
+
+ for (int i = 0; i < num_tuples; i++)
+ result[i] = tmp * (field[num_tuples+i] + field[2*num_tuples+i]);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::jacobianFcn(const double *, const double *, const int, const int ,
+ double *work, double *result)
+ {
+ // jacobian is cached in work array
+ assert(work);
+ std::copy(work, work+9, result);
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
+ const double *posn, const double *verts, const int nverts, const int ndim,
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
+ {
+ assert(posn && verts);
+ return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
+ params, is_inside);
+ }
+
+ bool LinearTri::insideFcn(const double *params, const int , const double tol)
+ {
+ return (params[0] >= -1.0-tol && params[1] >= -1.0-tol &&
+ params[0] + params[1] <= 1.0+tol);
+
+ }
+
+ ErrorCode LinearTri::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
+ const double *posn, const double *verts, const int nverts,
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside) {
+ // TODO: should differentiate between epsilons used for
+ // Newton Raphson iteration, and epsilons used for curved boundary geometry errors
+ // right now, fix the tolerance used for NR
+ const double error_tol_sqr = iter_tol*iter_tol;
+ CartVect *cvparams = reinterpret_cast<CartVect*>(params);
+ const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
+
+ // find best initial guess to improve convergence
+ CartVect tmp_params[] = {CartVect(-1,-1,-1), CartVect(1,-1,-1), CartVect(-1,1,-1)};
+ double resl = HUGE;
+ CartVect new_pos, tmp_pos;
+ ErrorCode rval;
+ for (unsigned int i = 0; i < 3; i++) {
+ rval = (*eval)(tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ double tmp_resl = (tmp_pos-*cvposn).length_squared();
+ if (tmp_resl < resl) {
+ *cvparams = tmp_params[i];
+ new_pos = tmp_pos;
+ resl = tmp_resl;
+ }
+ }
+
+ // residual is diff between old and new pos; need to minimize that
+ CartVect res = new_pos - *cvposn;
+ Matrix3 J;
+ rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
+ double det = J.determinant();
+ assert(det > std::numeric_limits<double>::epsilon());
+ Matrix3 Ji = J.inverse(1.0/det);
+
+ int iters=0;
+ // while |res| larger than tol
+ while (res % res > error_tol_sqr) {
+ if(++iters>25)
+ return MB_FAILURE;
+
+ // new params tries to eliminate residual
+ *cvparams -= Ji * res;
+
+ // get the new forward-evaluated position, and its difference from the target pt
+ rval = (*eval)(params, verts, ndim, 3, work, new_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ res = new_pos - *cvposn;
+ }
+
+ if (inside)
+ *inside = (*inside_f)(params, ndim, inside_tol);
+
+ return MB_SUCCESS;
+ }// Map::evaluate_reverse()
+
+} // namespace moab
diff --git a/src/LocalDiscretization/Makefile.am b/src/LocalDiscretization/Makefile.am
index 7309aef..cb34071 100644
--- a/src/LocalDiscretization/Makefile.am
+++ b/src/LocalDiscretization/Makefile.am
@@ -17,6 +17,7 @@ MOAB_LOCALDISCR_SRCS = \
LinearHex.cpp \
LinearQuad.cpp \
LinearTet.cpp \
+ LinearTri.cpp \
QuadraticHex.cpp
MOAB_LOCALDISCR_HDRS = \
@@ -24,6 +25,7 @@ MOAB_LOCALDISCR_HDRS = \
moab/LinearHex.hpp \
moab/LinearQuad.hpp \
moab/LinearTet.hpp \
+ moab/LinearTri.hpp \
moab/QuadraticHex.hpp
# The list of source files, and any header files that do not need to be installed
diff --git a/src/LocalDiscretization/QuadraticHex.cpp b/src/LocalDiscretization/QuadraticHex.cpp
index ca4758b..5c621f6 100644
--- a/src/LocalDiscretization/QuadraticHex.cpp
+++ b/src/LocalDiscretization/QuadraticHex.cpp
@@ -107,10 +107,12 @@ namespace moab
ErrorCode QuadraticHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
+ work, params, is_inside);
}
bool QuadraticHex::insideFcn(const double *params, const int ndim, const double tol)
diff --git a/src/LocalDiscretization/SpectralHex.cpp b/src/LocalDiscretization/SpectralHex.cpp
index 0e40eba..31f5391 100644
--- a/src/LocalDiscretization/SpectralHex.cpp
+++ b/src/LocalDiscretization/SpectralHex.cpp
@@ -73,7 +73,8 @@ CartVect SpectralHex::evaluate( const CartVect& params ) const
return result;
}
// replicate the functionality of hex_findpt
-bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double tol, const CartVect &init) const
+ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double iter_tol, const double inside_tol,
+ const CartVect &init) const
{
params = init;
@@ -92,7 +93,7 @@ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, doubl
//copy parametric coords back
params = r;
- return is_inside(params, tol);
+ return is_inside(params, inside_tol);
}
Matrix3 SpectralHex::jacobian(const CartVect& params) const
{
diff --git a/src/LocalDiscretization/SpectralQuad.cpp b/src/LocalDiscretization/SpectralQuad.cpp
index cbd1251..5c11d13 100644
--- a/src/LocalDiscretization/SpectralQuad.cpp
+++ b/src/LocalDiscretization/SpectralQuad.cpp
@@ -95,7 +95,8 @@ CartVect SpectralQuad::evalFcn(const double *params, const double *field, const
}
// replicate the functionality of hex_findpt
bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
params = init;
@@ -114,7 +115,7 @@ bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const
//copy parametric coords back
params = r;
- return insideFcn(params, 2, tol);
+ return insideFcn(params, 2, inside_tol);
}
diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index 48d34ac..29d1051 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -25,7 +25,8 @@ namespace moab {
typedef ErrorCode (*ReverseEvalFcn)(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *is_inside);
class EvalSet
{
@@ -76,8 +77,8 @@ namespace moab {
/** \brief Common function to do reverse evaluation based on evaluation and jacobian functions */
static ErrorCode evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
- const int ndim, const double tol, double *work, double *params,
- bool *inside);
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *inside);
/** \brief Common function that returns true if params is in [-1,1]^ndims */
static bool inside_function(const double *params, const int ndims, const double tol);
};
@@ -119,9 +120,9 @@ namespace moab {
* \param impl MOAB instance
* \param ent Entity handle to cache on the evaluator
* \param tag Tag to cache on the evaluator
- * \param tag_dim Tag dimension to cache on the evaluator
+ * \param tagged_ent_dim Dimension of entities to be tagged to cache on the evaluator
*/
- ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tag_dim = -1);
+ ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1);
/** \brief Evaluate cached tag at a given parametric location within the cached entity
* If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
@@ -133,12 +134,14 @@ namespace moab {
/** \brief Reverse-evaluate the cached entity at a given physical position
* \param posn Position at which to evaluate parameters
- * \param tol Tolerance of reverse evaluation, usually 10^-6 or so
+ * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
+ * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
* \param params Result of evaluation
* \param is_inside If non-NULL, returns true of resulting parameters place the point inside the element
* (in most cases, within [-1]*(dim)
*/
- ErrorCode reverse_eval(const double *posn, double tol, double *params, bool *is_inside = NULL) const;
+ ErrorCode reverse_eval(const double *posn, double iter_tol, double inside_tol, double *params,
+ bool *is_inside = NULL) const;
/** \brief Evaluate the jacobian of the cached entity at a given parametric location
* \param params Parameters at which to evaluate jacobian
@@ -166,14 +169,16 @@ namespace moab {
* object is changed.
* \param entities Entities tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(Range &entities, const double *point, double tol,
+ ErrorCode find_containing_entity(Range &entities, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -186,14 +191,16 @@ namespace moab {
* object is changed.
* \param ent_set Entity set containing the entities to be tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ ErrorCode find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -218,8 +225,12 @@ namespace moab {
inline ErrorCode set_ent_handle(EntityHandle ent);
/** \brief Get entity handle for this ElemEval */
- inline EntityHandle get_ent_handle() const {return entHandle;};
+ inline EntityHandle get_ent_handle() const {return entHandle;}
+ /* \brief Get vertex positions cached on this EE
+ */
+ inline double *get_vert_pos() {return vertPos[0].array();}
+
/* \brief Get the vertex handles cached here */
inline const EntityHandle *get_vert_handles() const {return vertHandles;}
@@ -233,24 +244,29 @@ namespace moab {
* To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
* and a tag dimension of 0.
* \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag_handle(Tag tag, int tag_dim = -1);
+ inline ErrorCode set_tag_handle(Tag tag, int tagged_ent_dim = -1);
/* \brief Set the name of the tag to cache on this evaluator
* To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
* and a tag dimension of 0.
- * \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tag_name Tag handle to cache, or 0 to cache vertex positions
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag(const char *tag_name, int ent_dim = -1);
+ inline ErrorCode set_tag(const char *tag_name, int tagged_ent_dim = -1);
/* \brief Get the dimension of the entities on which tag is set */
- inline int get_tag_dim() const {return tagDim;};
+ inline int get_tagged_ent_dim() const {return taggedEntDim;};
/* \brief Set the dimension of entities having the tag */
- inline ErrorCode set_tag_dim(int dim);
+ inline ErrorCode set_tagged_ent_dim(int dim);
+ /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as tags on entities
+ * Can't be const because most of the functions (evaluate, integrate, etc.) take non-const work space *'s
+ */
+ inline double *get_work_space() {return workSpace;}
+
/* \brief MOAB interface cached on this evaluator */
inline Interface *get_moab() {return mbImpl;}
@@ -287,7 +303,7 @@ namespace moab {
int numTuples;
/** \brief Dimension of entities from which to grab tag */
- int tagDim;
+ int taggedEntDim;
/** \brief Tag space */
std::vector<unsigned char> tagSpace;
@@ -300,13 +316,13 @@ namespace moab {
}; // class ElemEvaluator
- inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tag_dim)
+ inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tagged_ent_dim)
: mbImpl(impl), entHandle(0), entType(MBMAXTYPE), entDim(-1), numVerts(0),
vertHandles(NULL), tagHandle(0), tagCoords(false), numTuples(0),
- tagDim(0), workSpace(NULL)
+ taggedEntDim(0), workSpace(NULL)
{
if (ent) set_ent_handle(ent);
- if (tag) set_tag_handle(tag, tag_dim);
+ if (tag) set_tag_handle(tag, tagged_ent_dim);
}
inline ErrorCode ElemEvaluator::set_ent_handle(EntityHandle ent)
@@ -329,18 +345,17 @@ namespace moab {
rval = set_tag_handle(tagHandle);
if (MB_SUCCESS != rval) return rval;
}
-
if (evalSets[entType].initFcn) return (*evalSets[entType].initFcn)(vertPos[0].array(), numVerts, workSpace);
return MB_SUCCESS;
}
- inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
- if (!tag && !tag_dim) {
+ if (!tag && !tagged_ent_dim) {
tagCoords = true;
numTuples = 3;
- tagDim = 0;
+ taggedEntDim = 0;
tagHandle = 0;
return rval;
}
@@ -355,14 +370,14 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? 0 : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? 0 : tagged_ent_dim);
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -371,14 +386,14 @@ namespace moab {
return rval;
}
- inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
if (!tag_name || strlen(tag_name) == 0) return MB_FAILURE;
Tag tag;
if (!strcmp(tag_name, "COORDS")) {
tagCoords = true;
- tagDim = 0;
+ taggedEntDim = 0;
numTuples = 3;
tagHandle = 0;
// can return here, because vertex coords already cached when entity handle set
@@ -399,15 +414,15 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? entDim : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? entDim : tagged_ent_dim);
}
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -431,16 +446,17 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].evalFcn)(params,
(tagCoords ? (const double*) vertPos[0].array(): (const double*)&tagSpace[0]),
- entDim, (-1 == num_tuples ? numTuples : num_tuples),
- workSpace, result);
+ entDim,
+ (-1 == num_tuples ? numTuples : num_tuples), workSpace, result);
}
- inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double tol, double *params, bool *ins) const
+ inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double iter_tol, const double inside_tol,
+ double *params, bool *ins) const
{
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].reverseEvalFcn)(evalSets[entType].evalFcn, evalSets[entType].jacobianFcn, evalSets[entType].insideFcn,
- posn, vertPos[0].array(), numVerts, entDim, tol, workSpace,
- params, ins);
+ posn, vertPos[0].array(), numVerts,
+ entDim, iter_tol, inside_tol, workSpace, params, ins);
}
/** \brief Evaluate the jacobian of the cached entity at a given parametric location */
@@ -456,7 +472,7 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType && (tagCoords || tagHandle));
ErrorCode rval = MB_SUCCESS;
if (!tagCoords) {
- if (0 == tagDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
+ if (0 == taggedEntDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
else rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, (void*)&tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -465,7 +481,8 @@ namespace moab {
workSpace, result);
}
- inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals)
{
@@ -473,7 +490,7 @@ namespace moab {
Range entities;
ErrorCode rval = mbImpl->get_entities_by_handle(ent_set, entities);
if (MB_SUCCESS != rval) return rval;
- else return find_containing_entity(entities, point, tol, containing_ent, params, num_evals);
+ else return find_containing_entity(entities, point, iter_tol, inside_tol, containing_ent, params, num_evals);
}
inline bool ElemEvaluator::inside(const double *params, const double tol) const
diff --git a/src/LocalDiscretization/moab/LinearHex.hpp b/src/LocalDiscretization/moab/LinearHex.hpp
index e7ca3f4..03ed33f 100644
--- a/src/LocalDiscretization/moab/LinearHex.hpp
+++ b/src/LocalDiscretization/moab/LinearHex.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
This diff is so big that we needed to truncate the remainder.
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