[MOAB-dev] commit/MOAB: iulian07: continue for evaluation errors for mbcoupler
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Apr 21 10:41:56 CDT 2014
1 new commit in MOAB:
https://bitbucket.org/fathomteam/moab/commits/746455ef8011/
Changeset: 746455ef8011
Branch: master
User: iulian07
Date: 2014-04-21 17:39:58
Summary: continue for evaluation errors for mbcoupler
the jacobian becomes negative for points outside the parametric element
be more lenient, print some debug info, and continue on
Affected #: 6 files
diff --git a/src/LocalDiscretization/SpectralQuad.cpp b/src/LocalDiscretization/SpectralQuad.cpp
index d7355e6..7b581cb 100644
--- a/src/LocalDiscretization/SpectralQuad.cpp
+++ b/src/LocalDiscretization/SpectralQuad.cpp
@@ -109,7 +109,10 @@ bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const
double dist = opt_findpt_2(&_data, (const double **)_xyz, x_star, r, &c);
// if it did not converge, get out with throw...
if (dist > 0.9e+30)
- throw Map::EvaluationError();
+ {
+ std::vector<CartVect> dummy;
+ throw Map::EvaluationError(CartVect(x_star), dummy);
+ }
//c tells us if we landed inside the element or exactly on a face, edge, or node
// also, dist shows the distance to the computed point.
//copy parametric coords back
diff --git a/src/moab/BoundBox.hpp b/src/moab/BoundBox.hpp
index c74707a..3f1f07f 100644
--- a/src/moab/BoundBox.hpp
+++ b/src/moab/BoundBox.hpp
@@ -14,6 +14,15 @@ namespace moab {
BoundBox(const CartVect &min, const CartVect &max) :
bMin(min), bMax(max) {}
BoundBox(const double *corners);
+ // constructor used in element maps
+ BoundBox(std::vector<CartVect> points): bMin(DBL_MAX), bMax(-DBL_MAX)
+ {
+ for (size_t i=0; i<points.size(); i++)
+ {
+ update_min( points[i].array() );
+ update_max( points[i].array() );
+ }
+ }
~BoundBox() {}
bool contains_point(const double *point, const double tol = 0.0) const;
diff --git a/tools/mbcoupler/Coupler.cpp b/tools/mbcoupler/Coupler.cpp
index afee8af..b868006 100644
--- a/tools/mbcoupler/Coupler.cpp
+++ b/tools/mbcoupler/Coupler.cpp
@@ -758,38 +758,53 @@ ErrorCode Coupler::nat_param(double xyz[3],
std::cout << "Problems getting coordinates of vertices\n";
return result;
}
-
+ CartVect pos(xyz);
if (etype == MBHEX) {
if (8==num_connect)
{
Element::LinearHex hexmap(coords_vert);
- tmp_nat_coords = hexmap.ievaluate(CartVect(xyz), epsilon);
- bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
- if (!inside)
+ if (!hexmap.inside_box( pos, epsilon))
+ continue;
+ try {
+ tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
+ bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
+ if (!inside) continue;
+ }
+ catch (Element::Map::EvaluationError) {
continue;
+ }
}
else if (27==num_connect)
{
Element::QuadraticHex hexmap(coords_vert);
- tmp_nat_coords = hexmap.ievaluate(CartVect(xyz), epsilon);
- bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
- if (!inside)
- continue;
+ if (!hexmap.inside_box( pos, epsilon))
+ continue;
+ try {
+ tmp_nat_coords = hexmap.ievaluate(pos, epsilon);
+ bool inside = hexmap.inside_nat_space(tmp_nat_coords, epsilon);
+ if (!inside) continue;
+ }
+ catch (Element::Map::EvaluationError) {
+ continue;
+ }
}
else // TODO this case not treated yet, no interpolation
continue;
}
else if (etype == MBTET){
Element::LinearTet tetmap(coords_vert);
- tmp_nat_coords = tetmap.ievaluate(CartVect(xyz));
+ // this is just a linear solve; unless degenerate, will not except
+ tmp_nat_coords = tetmap.ievaluate(pos);
bool inside = tetmap.inside_nat_space(tmp_nat_coords, epsilon);
if (!inside)
continue;
}
else if (etype == MBQUAD){
Element::LinearQuad quadmap(coords_vert);
+ if (!quadmap.inside_box( pos, epsilon))
+ continue;
try {
- tmp_nat_coords = quadmap.ievaluate(CartVect(xyz), epsilon);
+ tmp_nat_coords = quadmap.ievaluate(pos, epsilon);
bool inside = quadmap.inside_nat_space(tmp_nat_coords, epsilon);
if (!inside) continue;
}
diff --git a/tools/mbcoupler/ElemUtil.cpp b/tools/mbcoupler/ElemUtil.cpp
index cfc36e0..ee338c0 100644
--- a/tools/mbcoupler/ElemUtil.cpp
+++ b/tools/mbcoupler/ElemUtil.cpp
@@ -3,6 +3,7 @@
#include <assert.h>
#include "ElemUtil.hpp"
+#include "moab/BoundBox.hpp"
namespace moab {
namespace ElemUtil {
@@ -433,6 +434,15 @@ namespace Element {
this->vertex = v;
}
+ bool Map::inside_box(const CartVect & xi, double & tol) const
+ {
+ // bail out early, before doing an expensive NR iteration
+ // compute box
+ BoundBox box(this->vertex);
+ return box.contains_point(xi.array(), tol);
+
+ }
+
//
CartVect Map::ievaluate(const CartVect& x, double tol, const CartVect& x0) const {
// TODO: should differentiate between epsilons used for
@@ -448,12 +458,12 @@ namespace Element {
int iters=0;
while (delta % delta > error_tol_sqr) {
if(++iters>10)
- throw Map::EvaluationError();
+ throw Map::EvaluationError(x, vertex);
J = jacobian(xi);
det = J.determinant();
if (det < std::numeric_limits<double>::epsilon())
- throw Map::EvaluationError();
+ throw Map::EvaluationError(x, vertex);
xi -= J.inverse(1.0/det) * delta;
delta = evaluate( xi ) - x;
}
@@ -896,7 +906,10 @@ namespace Element {
real dist = opt_findpt_3(&_data, (const real **)_xyz, x_star, r, &c);
// if it did not converge, get out with throw...
if (dist > 0.9e+30)
- throw Map::EvaluationError();
+ {
+ std::vector<CartVect> dummy;
+ throw Map::EvaluationError(xyz, dummy);
+ }
//c tells us if we landed inside the element or exactly on a face, edge, or node
// also, dist shows the distance to the computed point.
//copy parametric coords back
@@ -1184,7 +1197,11 @@ namespace Element {
real dist = opt_findpt_2(&_data, (const real **)_xyz, x_star, r, &c);
// if it did not converge, get out with throw...
if (dist > 0.9e+30)
- throw Map::EvaluationError();
+ {
+ std::vector<CartVect> dummy;
+ throw Map::EvaluationError(xyz, dummy);
+ }
+
//c tells us if we landed inside the element or exactly on a face, edge, or node
// also, dist shows the distance to the computed point.
//copy parametric coords back
diff --git a/tools/mbcoupler/ElemUtil.hpp b/tools/mbcoupler/ElemUtil.hpp
index 219323b..b2a743a 100644
--- a/tools/mbcoupler/ElemUtil.hpp
+++ b/tools/mbcoupler/ElemUtil.hpp
@@ -103,10 +103,22 @@ namespace ElemUtil {
/**\brief Set vertices. */
virtual void set_vertices(const std::vector<CartVect>& v);
+ // will look at the box formed by vertex coordinates, and before doing any NR, bail out if necessary
+ virtual bool inside_box(const CartVect & xi, double & tol) const;
+
/* Exception thrown when an evaluation fails (e.g., ievaluate fails to converge). */
class EvaluationError {
public:
- EvaluationError(){};
+ EvaluationError(const CartVect & x, const std::vector<CartVect> & verts): p(x), vertices(verts){
+#ifndef NDEBUG
+ std::cout << "p:" << p << "\n vertices.size() " <<vertices.size() << "\n";
+ for (size_t i=0; i<vertices.size(); i++)
+ std::cout << vertices[i] << "\n";
+#endif
+ };
+ private:
+ CartVect p;
+ std::vector<CartVect> vertices;
};// class EvaluationError
/* Exception thrown when a bad argument is encountered. */
diff --git a/tools/mbcoupler/ElementTest.cpp b/tools/mbcoupler/ElementTest.cpp
index 4f37602..ca6e9e2 100644
--- a/tools/mbcoupler/ElementTest.cpp
+++ b/tools/mbcoupler/ElementTest.cpp
@@ -24,7 +24,58 @@ void test_tet() {
}// test_tet()
void test_hex() {
- moab::Element::LinearHex hex;
+ double positions[] =
+ {
+ 236.80706050970281, -139.55422526228017, 193.27999999999997,
+ 236.47511729348639, -141.33020962638582, 193.27999999999997,
+ 237.8457938295229, -142.57076074835663, 193.27999999999997,
+ 239.12702305519684, -139.96608933577852, 193.27999999999997,
+ 236.80841321361444, -139.55341321335499, 202.654,
+ 236.47655014713746, -141.32980272396816, 202.654,
+ 237.8477913707564, -142.57047282187165, 202.654,
+ 239.12865103844533, -139.96531051891105, 202.654
+ };
+ CartVect x(235.96518686964933, -142.43503000077749, 188.19999999999987);
+ std::vector<CartVect> vertices;
+ for (int i=0; i<8; i++)
+ vertices.push_back(CartVect(positions+3*i));
+
+ moab::Element::LinearHex hex(vertices);
+ double tol(0.0001);
+ if (hex.inside_box(x, tol))
+ {
+ CartVect nat_par = hex.ievaluate(x, 0.0001);
+ std::cout <<nat_par <<"\n";
+ }
+
+ double positions2[] =
+ { 49.890500000000024, -20.376134375374882, 312.72000000000003,
+ 52.015875000000044, -19.149048546996006, 312.72000000000003,
+ 48.430375821458099, -18.548796774572125, 312.72000000000003,
+ 47.717616239031223, -21.191360829777231, 312.72000000000003,
+ 49.890500000000024, -20.376134375374882, 322.88,
+ 52.015875000000044, -19.149048546996006, 322.88,
+ 48.429930354643275, -18.52828610485021, 322.88,
+ 47.720552036968819, -21.167591146685712, 322.88};
+
+ CartVect x2(51.469000000000015, -20.145482942833631, 317.80000000000001);
+
+ vertices.clear();
+ for (int i=0; i<8; i++)
+ vertices.push_back(CartVect(positions2+3*i));
+ moab::Element::LinearHex hex2(vertices);
+ if (hex2.inside_box(x2, tol))
+ {
+ try {
+ CartVect nat_par = hex.ievaluate(x, 0.0001);
+ std::cout <<nat_par <<"\n";
+ }
+ catch (Element::Map::EvaluationError) {
+ // nothing
+ }
+}
+
+
}// test_hex()
#include "moab/Core.hpp"
#include "moab/Range.hpp"
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