|
moab
|
Implement CAMAL geometry callbacks using smooth iMesh. More...
#include <SmoothFace.hpp>
Implement CAMAL geometry callbacks using smooth iMesh.
Definition at line 33 of file SmoothFace.hpp.
| moab::SmoothFace::SmoothFace | ( | Interface * | mb, |
| EntityHandle | surface_set, | ||
| GeomTopoTool * | gTool | ||
| ) |
Definition at line 45 of file SmoothFace.cpp.
: _mb(mb), _set(surface_set), _my_geomTopoTool(gTool), _evaluationsCounter(0) { //_smooth_face = NULL; //_mbOut->create_meshset(MESHSET_SET, _oSet); //will contain the // get also the obb_root if (_my_geomTopoTool) { _my_geomTopoTool->get_root(this->_set, _obb_root); if (debug_surf_eval1) _my_geomTopoTool->obb_tree()->stats(_obb_root, std::cout); } }
| moab::SmoothFace::~SmoothFace | ( | ) | [virtual] |
Definition at line 60 of file SmoothFace.cpp.
{
}
| void moab::SmoothFace::ac_at_edge | ( | CartVect & | fac, |
| CartVect & | eac, | ||
| int | edge_id | ||
| ) | [private] |
Definition at line 1533 of file SmoothFace.cpp.
{
double u, v, w;
switch (edge_id) {
case 0:
u = 0.0;
v = fac[1] / (fac[1] + fac[2]);//v = fac.y() / (fac.y() + fac.z());
w = 1.0 - v;
break;
case 1:
u = fac[0] / (fac[0] + fac[2]);// u = fac.x() / (fac.x() + fac.z());
v = 0.0;
w = 1.0 - u;
break;
case 2:
u = fac[0] / (fac[0] + fac[1]);//u = fac.x() / (fac.x() + fac.y());
v = 1.0 - u;
w = 0.0;
break;
default:
assert(0);
u = -1; // needed to eliminate warnings about used before set
v = -1; // needed to eliminate warnings about used before set
w = -1; // needed to eliminate warnings about used before set
break;
}
eac[0] = u;
eac[1] = v;
eac[2] = w; //= CartVect(u, v, w);
}
| void moab::SmoothFace::adjust_bounding_box | ( | CartVect & | vect | ) | [private] |
| void moab::SmoothFace::append_smooth_tags | ( | std::vector< Tag > & | smoothTags | ) |
Definition at line 93 of file SmoothFace.cpp.
{
// these are created locally, for each smooth face
smoothTags.push_back(_gradientTag);
smoothTags.push_back(_planeTag);
}
| double moab::SmoothFace::area | ( | ) | [virtual] |
Definition at line 64 of file SmoothFace.cpp.
{
// find the area of this entity
//assert(_smooth_face);
//double area1 = _smooth_face->area();
double totArea = 0.;
for (Range::iterator it = _triangles.begin(); it != _triangles.end(); it++)
{
EntityHandle tria = *it;
const EntityHandle * conn3;
int nnodes;
_mb->get_connectivity(tria, conn3, nnodes);
//
//double coords[9]; // store the coordinates for the nodes
//_mb->get_coords(conn3, 3, coords);
CartVect p[3];
_mb->get_coords(conn3, 3, (double*) &p[0]);
// need to compute the angles
// compute angles and the normal
//CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
CartVect AB(p[1] - p[0]);//(p2 - p1);
CartVect BC(p[2] - p[1]);//(p3 - p2);
CartVect normal = AB * BC;
totArea += normal.length() / 2;
}
return totArea;
}
| void moab::SmoothFace::bounding_box | ( | double | box_min[3], |
| double | box_max[3] | ||
| ) | [virtual] |
| ErrorCode moab::SmoothFace::compute_control_points_on_edges | ( | double | min_dot, |
| Tag | edgeCtrlTag, | ||
| Tag | markTag | ||
| ) |
Definition at line 157 of file SmoothFace.cpp.
{
_edgeCtrlTag = edgeCtrlTag;
_markTag = markTag;
// now, compute control points for all edges that are not marked already (they are no on the boundary!)
for (Range::iterator it = _edges.begin(); it != _edges.end(); it++)
{
EntityHandle edg = *it;
// is the edge marked? already computed
unsigned char tagVal = 0;
_mb->tag_get_data(_markTag, &edg, 1, &tagVal);
if (tagVal)
continue;
//double min_dot;
init_bezier_edge(edg, min_dot);
tagVal = 1;
_mb->tag_set_data(_markTag, &edg, 1, &tagVal);
}
return MB_SUCCESS;
}
| ErrorCode moab::SmoothFace::compute_internal_control_points_on_facets | ( | double | min_dot, |
| Tag | facetCtrlTag, | ||
| Tag | facetEdgeCtrlTag | ||
| ) |
Definition at line 483 of file SmoothFace.cpp.
{
// collect from each triangle the control points in order
//
_facetCtrlTag = facetCtrlTag;
_facetEdgeCtrlTag = facetEdgeCtrlTag;
for (Range::iterator it = _triangles.begin(); it != _triangles.end(); it++)
{
EntityHandle tri = *it;
// first get connectivity, and the edges
// we need a fast method to retrieve the adjacent edges to each triangle
const EntityHandle * conn3;
int nnodes;
ErrorCode rval = _mb->get_connectivity(tri, conn3, nnodes);
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
assert(3 == nnodes);
// would it be easier to do
CartVect vNode[3];// position at nodes
rval = _mb->get_coords(conn3, 3, (double*) &vNode[0]);
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
// get gradients (normal) at each node of triangle
CartVect NN[3];
rval = _mb->tag_get_data(_gradientTag, conn3, 3, &NN[0]);
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
EntityHandle edges[3];
int orient[3]; // + 1 or -1, if the edge is positive or negative within the face
rval = find_edges_orientations(edges, conn3, orient);// maybe we will set it?
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
// maybe we will store some tags with edges and their orientation with respect to
// a triangle;
CartVect P[3][5];
CartVect N[6], G[6];
// create the linear array for control points on edges, for storage (expensive !!!)
CartVect CP[9];
int index = 0;
// maybe store a tag / entity handle for edges?
for (int i = 0; i < 3; i++)
{
// populate P and N with the right vectors
int i1 = (i + 1) % 3; // the first node of the edge
int i2 = (i + 2) % 3; // the second node of the edge
N[2 * i] = NN[i1];
N[2 * i + 1] = NN[i2];
P[i][0] = vNode[i1];
rval = _mb->tag_get_data(_edgeCtrlTag, &edges[i], 1, &(P[i][1]));
// if sense is -1, swap 1 and 3 control points
if (orient[i] == -1)
{
CartVect tmp;
tmp = P[i][1];
P[i][1] = P[i][3];
P[i][3] = tmp;
}
P[i][4] = vNode[i2];
for (int j = 1; j < 4; j++)
CP[index++] = P[i][j];
// the first edge control points
}
// stat = facet->get_edge_control_points( P );
init_facet_control_points(N, P, G);
// what do we need to store in the tag control points?
rval = _mb->tag_set_data(_facetCtrlTag, &tri, 1, &G[0]);
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
// store here again the 9 control points on the edges
rval = _mb->tag_set_data(_facetEdgeCtrlTag, &tri, 1, &CP[0]);
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
// look at what we retrieve later
// adjust the bounding box
int j = 0;
for (j = 0; j < 3; j++)
adjust_bounding_box(vNode[j]);
// edge control points
for (j = 0; j < 9; j++)
adjust_bounding_box(CP[j]);
// internal facet control points
for (j = 0; j < 6; j++)
adjust_bounding_box(G[j]);
}
return MB_SUCCESS;
}
Definition at line 378 of file SmoothFace.cpp.
{
double defTangents[6] = { 0., 0., 0., 0., 0., 0. };
ErrorCode rval = _mb->tag_get_handle("TANGENTS", 6, MB_TYPE_DOUBLE,
_tangentsTag, MB_TAG_DENSE|MB_TAG_CREAT, &defTangents);
if (MB_SUCCESS != rval)
return MB_FAILURE;
// now, compute Tangents for all edges that are not on boundary, so they are not marked
for (Range::iterator it = _edges.begin(); it != _edges.end(); it++)
{
EntityHandle edg = *it;
int nnodes;
const EntityHandle * conn2;//
_mb->get_connectivity(edg, conn2, nnodes);
assert(nnodes == 2);
CartVect P[2]; // store the coordinates for the nodes
rval = _mb->get_coords(conn2, 2, (double *) &P[0]);
if (MB_SUCCESS != rval) return rval;
assert(rval==MB_SUCCESS);
CartVect T[2];
T[0] = P[1] - P[0];
T[0].normalize();
T[1] = T[0]; //
_mb->tag_set_data(_tangentsTag, &edg, 1, (double*) &T[0]);// set the tangents computed at every edge
}
return MB_SUCCESS;
}
Definition at line 639 of file SmoothFace.cpp.
{
// here, we will dump all control points from edges and facets (6 control points for each facet)
// we may also create some edges; maybe later...
// create a point3D file
// output a Point3D file (special visit format)
unsigned long setId = _mb->id_from_handle(_set);
char name[50] = { 0 };
sprintf(name, "%ldcontrol.Point3D", setId);// name should be something 2control.Point3D
std::ofstream point3DFile;
point3DFile.open(name);//("control.Point3D");
point3DFile << "# x y z \n";
std::ofstream point3DEdgeFile;
sprintf(name, "%ldcontrolEdge.Point3D", setId);//
point3DEdgeFile.open(name);//("controlEdge.Point3D");
point3DEdgeFile << "# x y z \n";
std::ofstream smoothPoints;
sprintf(name, "%ldsmooth.Point3D", setId);//
smoothPoints.open(name);//("smooth.Point3D");
smoothPoints << "# x y z \n";
CartVect controlPoints[3]; // edge control points
for (Range::iterator it = _edges.begin(); it != _edges.end(); it++)
{
EntityHandle edge = *it;
_mb->tag_get_data(_edgeCtrlTag, &edge, 1, (double*) &controlPoints[0]);
for (int i = 0; i < 3; i++)
{
CartVect & c = controlPoints[i];
point3DEdgeFile << std::setprecision(11) << c[0] << " " << c[1] << " "
<< c[2] << " \n";
}
}
CartVect controlTriPoints[6]; // triangle control points
CartVect P_facet[3];// result in 3 "mid" control points
for (Range::iterator it2 = _triangles.begin(); it2 != _triangles.end(); it2++)
{
EntityHandle tri = *it2;
_mb->tag_get_data(_facetCtrlTag, &tri, 1, (double*) &controlTriPoints[0]);
// draw a line of points between pairs of control points
int numPoints = 7;
for (int n = 0; n < numPoints; n++)
{
double a = 1. * n / (numPoints - 1);
double b = 1.0 - a;
P_facet[0] = a * controlTriPoints[3] + b * controlTriPoints[4];
//1,2,1
P_facet[1] = a * controlTriPoints[0] + b * controlTriPoints[5];
//1,1,2
P_facet[2] = a * controlTriPoints[1] + b * controlTriPoints[2];
for (int i = 0; i < 3; i++)
{
CartVect & c = P_facet[i];
point3DFile << std::setprecision(11) << c[0] << " " << c[1] << " "
<< c[2] << " \n";
}
}
// evaluate for each triangle a lattice of points
int N = 40;
for (int k = 0; k <= N; k++)
{
for (int m = 0; m <= N - k; m++)
{
int n = N - m - k;
CartVect areacoord(1. * k / N, 1. * m / N, 1. * n / N);
CartVect pt;
eval_bezier_patch(tri, areacoord, pt);
smoothPoints << std::setprecision(11) << pt[0] << " " << pt[1] << " "
<< pt[2] << " \n";
}
}
}
point3DFile.close();
smoothPoints.close();
point3DEdgeFile.close();
return;
}
| ErrorCode moab::SmoothFace::eval_bezier_patch | ( | EntityHandle | tri, |
| CartVect & | areacoord, | ||
| CartVect & | pt | ||
| ) |
Definition at line 773 of file SmoothFace.cpp.
{
//
// interpolate internal control points
CartVect gctrl_pts[6];
// get the control points facet->get_control_points( gctrl_pts );
//init_facet_control_points( N, P, G) ;
// what do we need to store in the tag control points?
ErrorCode rval = _mb->tag_get_data(_facetCtrlTag, &tri, 1, &gctrl_pts[0]);// get all 6 control points
assert(MB_SUCCESS == rval);
if (MB_SUCCESS != rval) return rval;
const EntityHandle * conn3;
int nnodes;
rval = _mb->get_connectivity(tri, conn3, nnodes);
assert(MB_SUCCESS == rval);
CartVect vN[3];
_mb->get_coords(conn3, 3, (double*) &vN[0]); // fill the coordinates of the vertices
if (fabs(areacoord[1] + areacoord[2]) < 1.0e-6)
{
pt = vN[0];
return MB_SUCCESS;
}
if (fabs(areacoord[0] + areacoord[2]) < 1.0e-6)
{
pt = vN[0];
return MB_SUCCESS;
}
if (fabs(areacoord[0] + areacoord[1]) < 1.0e-6)
{
pt = vN[0];
return MB_SUCCESS;
}
CartVect P_facet[3];
//2,1,1
P_facet[0] = (1.0e0 / (areacoord[1] + areacoord[2])) * (areacoord[1]
* gctrl_pts[3] + areacoord[2] * gctrl_pts[4]);
//1,2,1
P_facet[1] = (1.0e0 / (areacoord[0] + areacoord[2])) * (areacoord[0]
* gctrl_pts[0] + areacoord[2] * gctrl_pts[5]);
//1,1,2
P_facet[2] = (1.0e0 / (areacoord[0] + areacoord[1])) * (areacoord[0]
* gctrl_pts[1] + areacoord[1] * gctrl_pts[2]);
// sum the contribution from each of the control points
pt = CartVect(0.);// set all to 0, we start adding / accumulating different parts
// first edge is from node 0 to 1, index 2 in
// retrieve the points, in order, and the control points on edges
// store here again the 9 control points on the edges
CartVect CP[9];
rval = _mb->tag_get_data(_facetEdgeCtrlTag, &tri, 1, &CP[0]);
assert(MB_SUCCESS == rval);
//CubitFacetEdge *edge;
//edge = facet->edge(2);! start with edge 2, from 0-1
int k = 0;
CartVect ctrl_pts[5];
//edge->control_points(facet, ctrl_pts);
ctrl_pts[0] = vN[0]; //
for (k = 1; k < 4; k++)
ctrl_pts[k] = CP[k + 5]; // for edge index 2
ctrl_pts[4] = vN[1]; //
//i=4; j=0; k=0;
double B = quart(areacoord[0]);
pt += B * ctrl_pts[0];
//i=3; j=1; k=0;
B = 4.0 * cube(areacoord[0]) * areacoord[1];
pt += B * ctrl_pts[1];
//i=2; j=2; k=0;
B = 6.0 * sqr(areacoord[0]) * sqr(areacoord[1]);
pt += B * ctrl_pts[2];
//i=1; j=3; k=0;
B = 4.0 * areacoord[0] * cube(areacoord[1]);
pt += B * ctrl_pts[3];
//edge = facet->edge(0);
//edge->control_points(facet, ctrl_pts);
// edge index 0, from 1 to 2
ctrl_pts[0] = vN[1]; //
for (k = 1; k < 4; k++)
ctrl_pts[k] = CP[k - 1]; // for edge index 0
ctrl_pts[4] = vN[2]; //
//i=0; j=4; k=0;
B = quart(areacoord[1]);
pt += B * ctrl_pts[0];
//i=0; j=3; k=1;
B = 4.0 * cube(areacoord[1]) * areacoord[2];
pt += B * ctrl_pts[1];
//i=0; j=2; k=2;
B = 6.0 * sqr(areacoord[1]) * sqr(areacoord[2]);
pt += B * ctrl_pts[2];
//i=0; j=1; k=3;
B = 4.0 * areacoord[1] * cube(areacoord[2]);
pt += B * ctrl_pts[3];
//edge = facet->edge(1);
//edge->control_points(facet, ctrl_pts);
// edge index 1, from 2 to 0
ctrl_pts[0] = vN[2]; //
for (k = 1; k < 4; k++)
ctrl_pts[k] = CP[k + 2]; // for edge index 0
ctrl_pts[4] = vN[0]; //
//i=0; j=0; k=4;
B = quart(areacoord[2]);
pt += B * ctrl_pts[0];
//i=1; j=0; k=3;
B = 4.0 * areacoord[0] * cube(areacoord[2]);
pt += B * ctrl_pts[1];
//i=2; j=0; k=2;
B = 6.0 * sqr(areacoord[0]) * sqr(areacoord[2]);
pt += B * ctrl_pts[2];
//i=3; j=0; k=1;
B = 4.0 * cube(areacoord[0]) * areacoord[2];
pt += B * ctrl_pts[3];
//i=2; j=1; k=1;
B = 12.0 * sqr(areacoord[0]) * areacoord[1] * areacoord[2];
pt += B * P_facet[0];
//i=1; j=2; k=1;
B = 12.0 * areacoord[0] * sqr(areacoord[1]) * areacoord[2];
pt += B * P_facet[1];
//i=1; j=1; k=2;
B = 12.0 * areacoord[0] * areacoord[1] * sqr(areacoord[2]);
pt += B * P_facet[2];
return MB_SUCCESS;
}
| ErrorCode moab::SmoothFace::eval_bezier_patch_normal | ( | EntityHandle | facet, |
| CartVect & | areacoord, | ||
| CartVect & | normal | ||
| ) |
Definition at line 1738 of file SmoothFace.cpp.
{
// interpolate internal control points
CartVect gctrl_pts[6];
//facet->get_control_points( gctrl_pts );
ErrorCode rval = _mb->tag_get_data(_facetCtrlTag, &facet, 1, &gctrl_pts[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
// _gradientTag
// get normals at points
const EntityHandle * conn3;
int nnodes;
rval = _mb->get_connectivity(facet, conn3, nnodes);
if (MB_SUCCESS != rval) return rval;
CartVect NN[3];
rval = _mb->tag_get_data(_gradientTag, conn3, 3, &NN[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
if (fabs(areacoord[1] + areacoord[2]) < 1.0e-6)
{
normal = NN[0];
return MB_SUCCESS;
}
if (fabs(areacoord[0] + areacoord[2]) < 1.0e-6)
{
normal = NN[1];//facet->point(1)->normal(facet);
return MB_SUCCESS;
}
if (fabs(areacoord[0] + areacoord[1]) < 1.0e-6)
{
normal = NN[2]; //facet->point(2)->normal(facet);
return MB_SUCCESS;
}
// compute the hodograph of the quartic Gregory patch
CartVect Nijk[10];
//hodograph(facet,areacoord,Nijk);
// start copy from hodograph
//CubitVector gctrl_pts[6];
// facet->get_control_points( gctrl_pts );
CartVect P_facet[3];
//2,1,1
/*P_facet[0] = (1.0e0 / (areacoord.y() + areacoord.z())) *
(areacoord.y() * gctrl_pts[3] +
areacoord.z() * gctrl_pts[4]);*/
P_facet[0] = (1.0e0 / (areacoord[1] + areacoord[2])) * (areacoord[1]
* gctrl_pts[3] + areacoord[2] * gctrl_pts[4]);
//1,2,1
/*P_facet[1] = (1.0e0 / (areacoord.x() + areacoord.z())) *
(areacoord.x() * gctrl_pts[0] +
areacoord.z() * gctrl_pts[5]);*/
P_facet[1] = (1.0e0 / (areacoord[0] + areacoord[2])) * (areacoord[0]
* gctrl_pts[0] + areacoord[2] * gctrl_pts[5]);
//1,1,2
/*P_facet[2] = (1.0e0 / (areacoord.x() + areacoord.y())) *
(areacoord.x() * gctrl_pts[1] +
areacoord.y() * gctrl_pts[2]);*/
P_facet[2] = (1.0e0 / (areacoord[0] + areacoord[1])) * (areacoord[0]
* gctrl_pts[1] + areacoord[1] * gctrl_pts[2]);
// corner control points are just the normals at the points
//3, 0, 0
Nijk[0] = NN[0];
//0, 3, 0
Nijk[3] = NN[1];
//0, 0, 3
Nijk[9] = NN[2];//facet->point(2)->normal(facet);
// fill in the boundary control points. Define as the normal to the local
// triangle formed by the quartic control point lattice
// store here again the 9 control points on the edges
CartVect CP[9]; // 9 control points on the edges,
rval = _mb->tag_get_data(_facetEdgeCtrlTag, &facet, 1, &CP[0]);
if (MB_SUCCESS != rval) return rval;
// there are 3 CP for each edge, 0, 1, 2; first edge is 1-2
//CubitFacetEdge *edge;
//edge = facet->edge( 2 );
//CubitVector ctrl_pts[5];
//edge->control_points(facet, ctrl_pts);
//2, 1, 0
//Nijk[1] = (ctrl_pts[2] - ctrl_pts[1]) * (P_facet[0] - ctrl_pts[1]);
Nijk[1] = (CP[7] - CP[6]) * (P_facet[0] - CP[6]);
Nijk[1].normalize();
//1, 2, 0
//Nijk[2] = (ctrl_pts[3] - ctrl_pts[2]) * (P_facet[1] - ctrl_pts[2]);
Nijk[2] = (CP[8] - CP[7]) * (P_facet[1] - CP[7]);
Nijk[2].normalize();
//edge = facet->edge( 0 );
//edge->control_points(facet, ctrl_pts);
//0, 2, 1
//Nijk[6] = (ctrl_pts[1] - P_facet[1]) * (ctrl_pts[2] - P_facet[1]);
Nijk[6] = (CP[0] - P_facet[1]) * (CP[1] - P_facet[1]);
Nijk[6].normalize();
//0, 1, 2
//Nijk[8] = (ctrl_pts[2] - P_facet[2]) * (ctrl_pts[3] - P_facet[2]);
Nijk[8] = (CP[1] - P_facet[2]) * (CP[2] - P_facet[2]);
Nijk[8].normalize();
//edge = facet->edge( 1 );
//edge->control_points(facet, ctrl_pts);
//1, 0, 2
//Nijk[7] = (P_facet[2] - ctrl_pts[2]) * (ctrl_pts[1] - ctrl_pts[2]);
Nijk[7] = (P_facet[2] - CP[4]) * (CP[3] - CP[4]);
Nijk[7].normalize();
//2, 0, 1
//Nijk[4] = (P_facet[0] - ctrl_pts[3]) * (ctrl_pts[2] - ctrl_pts[3]);
Nijk[4] = (P_facet[0] - CP[5]) * (CP[4] - CP[5]);
Nijk[4].normalize();
//1, 1, 1
Nijk[5] = (P_facet[1] - P_facet[0]) * (P_facet[2] - P_facet[0]);
Nijk[5].normalize();
// end copy from hodograph
// sum the contribution from each of the control points
normal = CartVect(0.0e0, 0.0e0, 0.0e0);
//i=3; j=0; k=0;
double Bsum = 0.0;
double B = cube(areacoord[0]);
Bsum += B;
normal += B * Nijk[0];
//i=2; j=1; k=0;
B = 3.0 * sqr(areacoord[0]) * areacoord[1];
Bsum += B;
normal += B * Nijk[1];
//i=1; j=2; k=0;
B = 3.0 * areacoord[0] * sqr(areacoord[1]);
Bsum += B;
normal += B * Nijk[2];
//i=0; j=3; k=0;
B = cube(areacoord[1]);
Bsum += B;
normal += B * Nijk[3];
//i=2; j=0; k=1;
B = 3.0 * sqr(areacoord[0]) * areacoord[2];
Bsum += B;
normal += B * Nijk[4];
//i=1; j=1; k=1;
B = 6.0 * areacoord[0] * areacoord[1] * areacoord[2];
Bsum += B;
normal += B * Nijk[5];
//i=0; j=2; k=1;
B = 3.0 * sqr(areacoord[1]) * areacoord[2];
Bsum += B;
normal += B * Nijk[6];
//i=1; j=0; k=2;
B = 3.0 * areacoord[0] * sqr(areacoord[2]);
Bsum += B;
normal += B * Nijk[7];
//i=0; j=1; k=2;
B = 3.0 * areacoord[1] * sqr(areacoord[2]);
Bsum += B;
normal += B * Nijk[8];
//i=0; j=0; k=3;
B = cube(areacoord[2]);
Bsum += B;
normal += B * Nijk[9];
//assert(fabs(Bsum - 1.0) < 1e-9);
normal.normalize();
return MB_SUCCESS;
}
| int moab::SmoothFace::eval_counter | ( | ) | [inline] |
Definition at line 119 of file SmoothFace.hpp.
{
return _evaluationsCounter;
}
| ErrorCode moab::SmoothFace::evaluate_smooth_edge | ( | EntityHandle | eh, |
| double & | t, | ||
| CartVect & | outv | ||
| ) |
Definition at line 729 of file SmoothFace.cpp.
{
CartVect P[2]; // P0 and P1
CartVect controlPoints[3]; // edge control points
double t4, t3, t2, one_minus_t, one_minus_t2, one_minus_t3, one_minus_t4;
// project the position to the linear edge
// t is from 0 to 1 only!!
//double tt = (t + 1) * 0.5;
if (tt <= 0.0)
tt = 0.0;
if (tt >= 1.0)
tt = 1.0;
int nnodes;
const EntityHandle * conn2;
ErrorCode rval = _mb->get_connectivity(eh, conn2, nnodes);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
rval = _mb->get_coords(conn2, 2, (double*) &P[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
rval = _mb->tag_get_data(_edgeCtrlTag, &eh, 1, (double*) &controlPoints[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
t2 = tt * tt;
t3 = t2 * tt;
t4 = t3 * tt;
one_minus_t = 1. - tt;
one_minus_t2 = one_minus_t * one_minus_t;
one_minus_t3 = one_minus_t2 * one_minus_t;
one_minus_t4 = one_minus_t3 * one_minus_t;
outv = one_minus_t4 * P[0] + 4. * one_minus_t3 * tt * controlPoints[0] + 6.
* one_minus_t2 * t2 * controlPoints[1] + 4. * one_minus_t * t3
* controlPoints[2] + t4 * P[1];
return MB_SUCCESS;
}
| void moab::SmoothFace::facet_area_coordinate | ( | EntityHandle | facet, |
| CartVect & | pt_on_plane, | ||
| CartVect & | areacoord | ||
| ) |
Definition at line 954 of file SmoothFace.cpp.
{
const EntityHandle * conn3;
int nnodes;
ErrorCode rval = _mb->get_connectivity(facet, conn3, nnodes);
assert(MB_SUCCESS == rval);
if (rval) {} // empty statement to prevent compiler warning
//double coords[9]; // store the coordinates for the nodes
//_mb->get_coords(conn3, 3, coords);
CartVect p[3];
rval = _mb->get_coords(conn3, 3, (double*) &p[0]);
assert(MB_SUCCESS == rval);
if (rval) {} // empty statement to prevent compiler warning
double plane[4];
rval = _mb->tag_get_data(_planeTag, &facet, 1, plane);
assert(rval == MB_SUCCESS);
if (rval) {} // empty statement to prevent compiler warning
CartVect normal(&plane[0]);// just first 3 components are used
double area2;
double tol = GEOMETRY_RESABS * 1.e-5;// 1.e-11;
CartVect v1(p[1] - p[0]);
CartVect v2(p[2] - p[0]);
area2 = (v1 * v2).length_squared();// the same for CartVect
if (area2 < 100 * tol)
{
tol = .01 * area2;
}
CartVect absnorm(fabs(normal[0]), fabs(normal[1]), fabs(normal[2]));
// project to the closest coordinate plane so we only have to do this in 2D
if (absnorm[0] >= absnorm[1] && absnorm[0] >= absnorm[2])
{
area2 = determ3( p[0][1], p[0][2],
p[1][1], p[1][2],
p[2][1], p[2][2]);
if (fabs(area2) < tol)
{
areacoord = CartVect(-CUBIT_DBL_MAX);// .set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
}
else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
}
else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
}
else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
}
else
{
areacoord[0] = determ3(pt_on_plane[1], pt_on_plane[2],
p[1][1], p[1][2], p[2][1], p[2][2] ) / area2;
areacoord[1] = determ3( p[0][1], p[0][2],
pt_on_plane[1], pt_on_plane[2],
p[2][1], p[2][2]) / area2;
areacoord[2] = determ3( p[0][1], p[0][2], p[1][1], p[1][2],
pt_on_plane[1], pt_on_plane[2]) / area2;
}
}
else if (absnorm[1] >= absnorm[0] && absnorm[1] >= absnorm[2])
{
area2 = determ3(p[0][0], p[0][2],
p[1][0], p[1][2],
p[2][0], p[2][2]);
if (fabs(area2) < tol)
{
areacoord = CartVect(-CUBIT_DBL_MAX);//.set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
}
else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
}
else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
}
else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
}
else
{
areacoord[0] = determ3(pt_on_plane[0], pt_on_plane[2],
p[1][0], p[1][2], p[2][0], p[2][2] ) / area2;
areacoord[1] = determ3( p[0][0], p[0][2],
pt_on_plane[0], pt_on_plane[2],
p[2][0], p[2][2]) / area2;
areacoord[2] = determ3( p[0][0], p[0][2], p[1][0], p[1][2],
pt_on_plane[0], pt_on_plane[2]) / area2;
}
}
else
{
/*area2 = determ3(pt0->x(), pt0->y(),
pt1->x(), pt1->y(),
pt2->x(), pt2->y());*/
area2 = determ3( p[0][0], p[0][1],
p[1][0], p[1][1],
p[2][0], p[2][1]);
if (fabs(area2) < tol)
{
areacoord = CartVect(-CUBIT_DBL_MAX);//.set( -CUBIT_DBL_MAX, -CUBIT_DBL_MAX, -CUBIT_DBL_MAX );
}
else if (within_tolerance(p[0], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(1., 0., 0.);//.set( 1.0, 0.0, 0.0 );
}
else if (within_tolerance(p[1], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 1., 0.);//.set( 0.0, 1.0, 0.0 );
}
else if (within_tolerance(p[2], pt_on_plane, GEOMETRY_RESABS))
{
areacoord = CartVect(0., 0., 1.);//.set( 0.0, 0.0, 1.0 );
}
else
{
areacoord[0] = determ3(pt_on_plane[0], pt_on_plane[1],
p[1][0], p[1][1], p[2][0], p[2][1] ) / area2;
areacoord[1] = determ3( p[0][0], p[0][1],
pt_on_plane[0], pt_on_plane[1],
p[2][0], p[2][1]) / area2;
areacoord[2] = determ3( p[0][0], p[0][1], p[1][0], p[1][1],
pt_on_plane[0], pt_on_plane[1]) / area2;
}
}
}
| ErrorCode moab::SmoothFace::find_edges_orientations | ( | EntityHandle | edges[3], |
| const EntityHandle * | conn3, | ||
| int | orient[3] | ||
| ) | [private] |
Definition at line 449 of file SmoothFace.cpp.
{
// find the edge that is adjacent to 2 vertices at a time
for (int i = 0; i < 3; i++)
{
// edge 0 is 1-2, 1 is 3-1, 2 is 0-1
EntityHandle v[2];
v[0] = conn3[(i + 1) % 3];
v[1] = conn3[(i + 2) % 3];
std::vector<EntityHandle> adjacencies;
// generate all edges for these two hexes
ErrorCode rval = _mb->get_adjacencies(v, 2, 1, false, adjacencies,
Interface::INTERSECT);
if (MB_SUCCESS != rval) return rval;
// find the edge connected to both vertices, and then see its orientation
assert(adjacencies.size() == 1);
const EntityHandle * conn2;
int nnodes;
rval = _mb->get_connectivity(adjacencies[0], conn2, nnodes);
assert(rval == MB_SUCCESS);
assert(2 == nnodes);
edges[i] = adjacencies[0];
// what is the story morning glory?
if (conn2[0] == v[0] && conn2[1] == v[1])
orient[i] = 1;
else if (conn2[0] == v[1] && conn2[1] == v[0])
orient[i] = -1;
else
return MB_FAILURE;
}
return MB_SUCCESS;
}
| ErrorCode moab::SmoothFace::get_normals_for_vertices | ( | const EntityHandle * | conn2, |
| CartVect | N[2] | ||
| ) |
Definition at line 1929 of file SmoothFace.cpp.
{
//CartVect N[2];
//_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
ErrorCode rval = _mb->tag_get_data(_gradientTag, conn2, 2, (double*) &N[0]);
return rval;
}
| ErrorCode moab::SmoothFace::init_bezier_edge | ( | EntityHandle | edge, |
| double | min_dot | ||
| ) | [private] |
Definition at line 319 of file SmoothFace.cpp.
{
// min dot was used for angle here
//int stat = 0; // CUBIT_SUCCESS;
// all boundaries will be simple, initially
// we may complicate them afterwards
CartVect ctrl_pts[3];
int nnodes;
const EntityHandle * conn2;
ErrorCode rval = _mb->get_connectivity(edge, conn2, nnodes);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
assert(2 == nnodes);
//double coords[6]; // store the coordinates for the nodes
CartVect P[2];
//ErrorCode rval = _mb->get_coords(conn2, 2, coords);
rval = _mb->get_coords(conn2, 2, (double*) &P[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
//CartVect P0(&coords[0]);
//CartVect P3(&coords[3]);
//double normalVec[6];
CartVect N[2];
//_mb->tag_get_data(_gradientTag, conn2, 2, normalVec);
rval = _mb->tag_get_data(_gradientTag, conn2, 2, (double*) &N[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
CartVect T[2]; // T0, T3
rval = _mb->tag_get_data(_tangentsTag, &edge, 1, &T[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
rval = init_edge_control_points(P[0], P[1], N[0], N[1], T[0], T[1], ctrl_pts);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
rval = _mb->tag_set_data(_edgeCtrlTag, &edge, 1, &ctrl_pts[0]);
assert(rval == MB_SUCCESS);
if (MB_SUCCESS != rval) return rval;
if (debug_surf_eval1)
{
std::cout << "edge: " << _mb-> id_from_handle(edge) << " tangents: "
<< T[0] << T[1] << std::endl;
std::cout << " points: " << P[0] << " " << P[1] << std::endl;
std::cout << " normals: " << N[0] << " " << N[1] << std::endl;
std::cout << " Control points " << ctrl_pts[0] << " " << ctrl_pts[1]
<< " " << ctrl_pts[2] << std::endl;
}
return MB_SUCCESS;
}
| ErrorCode moab::SmoothFace::init_edge_control_points | ( | CartVect & | P0, |
| CartVect & | P3, | ||
| CartVect & | N0, | ||
| CartVect & | N3, | ||
| CartVect & | T0, | ||
| CartVect & | T3, | ||
| CartVect * | Pi | ||
| ) |
Definition at line 415 of file SmoothFace.cpp.
{
CartVect Vi[4];
Vi[0] = P0;
Vi[3] = P3;
CartVect P03(P3 - P0);
double di = P03.length();
double ai = N0 % N3;// this is the dot operator, the same as in cgm for CubitVector
double ai0 = N0 % T0;
double ai3 = N3 % T3;
double denom = 4 - ai * ai;
if (fabs(denom) < 1e-20)
{
return MB_FAILURE; // CUBIT_FAILURE;
}
double row = 6.0e0 * (2.0e0 * ai0 + ai * ai3) / denom;
double omega = 6.0e0 * (2.0e0 * ai3 + ai * ai0) / denom;
Vi[1] = Vi[0] + (di * (((6.0e0 * T0) - ((2.0e0 * row) * N0) + (omega * N3))
/ 18.0e0));
Vi[2] = Vi[3] - (di * (((6.0e0 * T3) + (row * N0) - ((2.0e0 * omega) * N3))
/ 18.0e0));
CartVect Wi[3];
Wi[0] = Vi[1] - Vi[0];
Wi[1] = Vi[2] - Vi[1];
Wi[2] = Vi[3] - Vi[2];
Pi[0] = 0.25 * Vi[0] + 0.75 * Vi[1];
Pi[1] = 0.50 * Vi[1] + 0.50 * Vi[2];
Pi[2] = 0.75 * Vi[2] + 0.25 * Vi[3];
return MB_SUCCESS;
}
| ErrorCode moab::SmoothFace::init_facet_control_points | ( | CartVect | N[6], |
| CartVect | P[3][5], | ||
| CartVect | G[6] | ||
| ) | [private] |
Definition at line 597 of file SmoothFace.cpp.
{
CartVect Di[4], Ai[3], N0, N3, Vi[4], Wi[3];
double denom;
double lambda[2], mu[2];
ErrorCode rval = MB_SUCCESS;
for (int i = 0; i < 3; i++)
{
N0 = N[i * 2];
N3 = N[i * 2 + 1];
Vi[0] = P[i][0];
Vi[1] = (P[i][1] - 0.25 * P[i][0]) / 0.75;
Vi[2] = (P[i][3] - 0.25 * P[i][4]) / 0.75;
Vi[3] = P[i][4];
Wi[0] = Vi[1] - Vi[0];
Wi[1] = Vi[2] - Vi[1];
Wi[2] = Vi[3] - Vi[2];
Di[0] = P[(i + 2) % 3][3] - 0.5 * (P[i][1] + P[i][0]);
Di[3] = P[(i + 1) % 3][1] - 0.5 * (P[i][4] + P[i][3]);
Ai[0] = (N0 * Wi[0]) / Wi[0].length();
Ai[2] = (N3 * Wi[2]) / Wi[2].length();
Ai[1] = Ai[0] + Ai[2];
denom = Ai[1].length();
Ai[1] /= denom;
lambda[0] = (Di[0] % Wi[0]) / (Wi[0] % Wi[0]);
lambda[1] = (Di[3] % Wi[2]) / (Wi[2] % Wi[2]);
mu[0] = (Di[0] % Ai[0]);
mu[1] = (Di[3] % Ai[2]);
G[i * 2] = 0.5 * (P[i][1] + P[i][2]) + 0.66666666666666 * lambda[0] * Wi[1]
+ 0.33333333333333 * lambda[1] * Wi[0] + 0.66666666666666 * mu[0]
* Ai[1] + 0.33333333333333 * mu[1] * Ai[0];
G[i * 2 + 1] = 0.5 * (P[i][2] + P[i][3]) + 0.33333333333333 * lambda[0]
* Wi[2] + 0.66666666666666 * lambda[1] * Wi[1] + 0.33333333333333
* mu[0] * Ai[2] + 0.66666666666666 * mu[1] * Ai[1];
}
return rval;
}
| int moab::SmoothFace::init_gradient | ( | ) |
Definition at line 181 of file SmoothFace.cpp.
{
// first, create a Tag for gradient (or normal)
// loop over all triangles in set, and modify the normals according to the angle as weight
// get all the edges from this subset
if (NULL == _mb)
return 1; //fail
_triangles.clear();
ErrorCode rval = _mb->get_entities_by_type(_set, MBTRI, _triangles);
if (MB_SUCCESS != rval)
return 1;
// get a new range of edges, and decide the loops from here
_edges.clear();
rval = _mb-> get_adjacencies(_triangles, 1, true, _edges, Interface::UNION);
assert(rval == MB_SUCCESS);
rval = _mb->get_adjacencies(_triangles, 0, false, _nodes, Interface::UNION);
assert(rval == MB_SUCCESS);
// initialize bounding box limits
CartVect vert1;
EntityHandle v1 = _nodes[0];// first vertex
_mb->get_coords(&v1, 1, (double*) &vert1);
_minim = vert1;
_maxim = vert1;
double defNormal[] = { 0., 0., 0. };
// look for a tag name here that is definitely unique. We do not want the tags to interfere with each other
// this normal will be for each node in a face
// some nodes have multiple normals, if they are at the feature edges
unsigned long setId = _mb->id_from_handle(_set);
char name[50] = { 0 };
sprintf(name, "GRADIENT%ld", setId);// name should be something like GRADIENT29, where 29 is the set ID of the face
rval = _mb->tag_get_handle(name, 3, MB_TYPE_DOUBLE, _gradientTag,
MB_TAG_DENSE|MB_TAG_CREAT, &defNormal);
double defPlane[4] = { 0., 0., 1., 0. };
// also define a plane tag ; this will be for each triangle
char namePlaneTag[50] = { 0 };
sprintf(namePlaneTag, "PLANE%ld", setId);
rval = _mb->tag_get_handle("PLANE", 4, MB_TYPE_DOUBLE, _planeTag,
MB_TAG_DENSE|MB_TAG_CREAT, &defPlane);
// the fourth double is for weight, accumulated at each vertex so far
// maybe not needed in the end
for (Range::iterator it = _triangles.begin(); it != _triangles.end(); it++)
{
EntityHandle tria = *it;
const EntityHandle * conn3;
int nnodes;
_mb->get_connectivity(tria, conn3, nnodes);
if (nnodes != 3)
return 1; // error
//double coords[9]; // store the coordinates for the nodes
//_mb->get_coords(conn3, 3, coords);
CartVect p[3];
_mb->get_coords(conn3, 3, (double*) &p[0]);
// need to compute the angles
// compute angles and the normal
//CartVect p1(&coords[0]), p2(&coords[3]), p3(&coords[6]);
CartVect AB(p[1] - p[0]);//(p2 - p1);
CartVect BC(p[2] - p[1]);//(p3 - p2);
CartVect CA(p[0] - p[2]);//(p1 - p3);
double a[3];
a[1] = angle(AB, -BC); // angle at B (p2), etc.
a[2] = angle(BC, -CA);
a[0] = angle(CA, -AB);
CartVect normal = -AB * CA;
normal.normalize();
double plane[4];
const double * coordNormal = normal.array();
plane[0] = coordNormal[0];
plane[1] = coordNormal[1];
plane[2] = coordNormal[2];
plane[3] = -normal % p[0]; // dot product
//set the plane
rval = _mb->tag_set_data(_planeTag, &tria, 1, plane);
assert(rval == MB_SUCCESS);
// add to each vertex the tag value the normal multiplied by the angle
double values[9];
_mb->tag_get_data(_gradientTag, conn3, 3, values);
for (int i = 0; i < 3; i++)
{
//values[4*i]+=a[i]; // first is the weight, which we do not really need
values[3 * i + 0] += a[i] * coordNormal[0];
values[3 * i + 1] += a[i] * coordNormal[1];
values[3 * i + 2] += a[i] * coordNormal[2];
}
// reset those values
_mb->tag_set_data(_gradientTag, conn3, 3, values);
}
// normalize the gradients at each node; maybe not needed here?
// no, we will do it, it is important
int numNodes = _nodes.size();
double * normalVal = new double[3 * numNodes];
_mb->tag_get_data(_gradientTag, _nodes, normalVal);// get all the normal values at the _nodes
for (int i = 0; i < numNodes; i++)
{
CartVect p1(&normalVal[3 * i]);
p1.normalize();
p1.get(&normalVal[3 * i]);
}
// reset the normal values after normalization
_mb->tag_set_data(_gradientTag, _nodes, normalVal);
// print the loops size and some other stuff
if (debug_surf_eval1)
{
std::cout << " normals at " << numNodes << " nodes" << std::endl;
int i = 0;
for (Range::iterator it = _nodes.begin(); it != _nodes.end(); it++, i++)
{
EntityHandle node = *it;
std::cout << " Node id " << _mb->id_from_handle(node) << " "
<< normalVal[3 * i] << " " << normalVal[3 * i + 1] << " "
<< normalVal[3 * i + 2] << std::endl;
}
}
delete[] normalVal;
return 0;
}
| bool moab::SmoothFace::is_at_vertex | ( | EntityHandle | facet, |
| CartVect & | pt, | ||
| CartVect & | ac, | ||
| double | compare_tol, | ||
| CartVect & | eval_pt, | ||
| CartVect * | eval_norm_ptr | ||
| ) |
Definition at line 1604 of file SmoothFace.cpp.
{
double dist;
CartVect vert_loc;
const double actol = 0.1;
// get coordinates get_coords
const EntityHandle * conn3;
int nnodes;
_mb->get_connectivity(facet, conn3, nnodes);
//
//double coords[9]; // store the coordinates for the nodes
//_mb->get_coords(conn3, 3, coords);
CartVect p[3];
_mb->get_coords(conn3, 3, (double*) &p[0]);
// also get the normals at nodes
CartVect NN[3];
_mb->tag_get_data(_gradientTag, conn3, 3, (double*) &NN[0]);
if (fabs(ac[0]) < actol && fabs(ac[1]) < actol)
{
vert_loc = p[2];
dist = (pt - vert_loc).length(); //pt.distance_between( vert_loc );
if (dist <= compare_tol)
{
eval_pt = vert_loc;
if (eval_norm_ptr)
{
*eval_norm_ptr = NN[2];
}
return true;
}
}
if (fabs(ac[0]) < actol && fabs(ac[2]) < actol)
{
vert_loc = p[1];
dist = (pt - vert_loc).length();//pt.distance_between( vert_loc );
if (dist <= compare_tol)
{
eval_pt = vert_loc;
if (eval_norm_ptr)
{
*eval_norm_ptr = NN[1];//facet->point(1)->normal( facet );
}
return true;
}
}
if (fabs(ac[1]) < actol && fabs(ac[2]) < actol)
{
vert_loc = p[0];
dist = (pt - vert_loc).length();//pt.distance_between( vert_loc );
if (dist <= compare_tol)
{
eval_pt = vert_loc;
if (eval_norm_ptr)
{
*eval_norm_ptr = NN[0];
}
return true;
}
}
return false;
}
| bool moab::SmoothFace::move_ac_inside | ( | CartVect & | ac, |
| double | tol | ||
| ) | [private] |
Definition at line 1682 of file SmoothFace.cpp.
{
int nout = 0;
if (ac[0] < -tol)
{
ac[0] = 0.0;
ac[1] = ac[1] / (ac[1] + ac[2]); //( ac.y() / (ac.y() + ac.z()) ;
ac[2] = 1. - ac[1]; //ac.z( 1.0 - ac.y() );
nout++;
}
if (ac[1] < -tol)
{
ac[1] = 0.;//ac.y( 0.0 );
ac[0] = ac[0] / (ac[0] + ac[2]);//ac.x( ac.x() / (ac.x() + ac.z()) );
ac[2] = 1. - ac[0];//ac.z( 1.0 - ac.x() );
nout++;
}
if (ac[2] < -tol)
{
ac[2] = 0.;// ac.z( 0.0 );
ac[0] = ac[0] / (ac[0] + ac[1]);//ac.x( ac.x() / (ac.x() + ac.y()) );
ac[1] = 1. - ac[0]; // ac.y( 1.0 - ac.x() );
nout++;
}
return (nout > 0) ? true : false;
}
| void moab::SmoothFace::move_to_surface | ( | double & | x, |
| double & | y, | ||
| double & | z | ||
| ) | [virtual] |
Definition at line 111 of file SmoothFace.cpp.
{
CartVect loc2(x, y, z);
bool trim = false;// is it needed?
bool outside = true;
CartVect closestPoint;
ErrorCode rval = project_to_facets_main(loc2, trim, outside, &closestPoint,
NULL);
if (MB_SUCCESS != rval) return;
assert(rval==MB_SUCCESS);
x = closestPoint[0];
y = closestPoint[1];
z = closestPoint[2];
}
| bool moab::SmoothFace::normal_at | ( | double | x, |
| double | y, | ||
| double | z, | ||
| double & | nx, | ||
| double & | ny, | ||
| double & | nz | ||
| ) | [virtual] |
Definition at line 136 of file SmoothFace.cpp.
{
CartVect loc2(x, y, z);
bool trim = false;// is it needed?
bool outside = true;
//CartVect closestPoint;// not needed
// not interested in normal
CartVect normal;
ErrorCode rval = project_to_facets_main(loc2, trim, outside, NULL, &normal);
if (MB_SUCCESS != rval) return false;
assert(rval==MB_SUCCESS);
nx = normal[0];
ny = normal[1];
nz = normal[2];
return true;
}
| ErrorCode moab::SmoothFace::project_to_facet | ( | EntityHandle | facet, |
| CartVect & | pt, | ||
| CartVect & | areacoord, | ||
| CartVect & | close_point, | ||
| bool & | outside_facet, | ||
| double | compare_tol | ||
| ) |
Definition at line 1573 of file SmoothFace.cpp.
{
ErrorCode stat = MB_SUCCESS;
const EntityHandle * conn3;
int nnodes;
_mb->get_connectivity(facet, conn3, nnodes);
//
//double coords[9]; // store the coordinates for the nodes
//_mb->get_coords(conn3, 3, coords);
CartVect p[3];
_mb->get_coords(conn3, 3, (double*) &p[0]);
int edge_id = -1;
stat = project_to_patch(facet, areacoord, pt, close_point, NULL,
outside_facet, compare_tol, edge_id);
/* }
break;
}*/
return stat;
}
| void moab::SmoothFace::project_to_facet_plane | ( | EntityHandle | tri, |
| CartVect & | pt, | ||
| CartVect & | point_on_plane, | ||
| double & | dist_to_plane | ||
| ) |
Definition at line 929 of file SmoothFace.cpp.
{
double plane[4];
ErrorCode rval = _mb->tag_get_data(_planeTag, &tri, 1, plane);
if (MB_SUCCESS != rval) return;
assert(rval == MB_SUCCESS);
// _planeTag
CartVect normal(&plane[0]);// just first 3 components are used
double dist = normal % pt + plane[3]; // coeff d is saved!!!
dist_to_plane = fabs(dist);
point_on_plane = pt - dist * normal;
return;
}
| ErrorCode moab::SmoothFace::project_to_facets | ( | std::vector< EntityHandle > & | facet_list, |
| EntityHandle & | lastFacet, | ||
| int | interpOrder, | ||
| double | compareTol, | ||
| CartVect & | this_point, | ||
| bool | trim, | ||
| bool & | outside, | ||
| CartVect * | closest_point_ptr, | ||
| CartVect * | normal_ptr | ||
| ) |
best_facet
Definition at line 1127 of file SmoothFace.cpp.
{
bool outside_facet, best_outside_facet = true;
double mindist = 1.e20;
CartVect close_point, best_point(mindist, mindist, mindist), best_areacoord;
EntityHandle best_facet = 0L;// no best facet found yet
EntityHandle facet;
assert(facet_list.size() > 0);
double big_dist = compareTol * 1.0e3;
// from the list of close facets, determine the closest point
for (size_t i = 0; i < facet_list.size(); i++)
{
facet = facet_list[i];
CartVect pt_on_plane;
double dist_to_plane;
project_to_facet_plane(facet, this_point, pt_on_plane, dist_to_plane);
CartVect areacoord;
//CartVect close_point;
facet_area_coordinate(facet, pt_on_plane, areacoord);
if (interpOrder != 0)
{
// modify the areacoord - project to the bezier patch- snaps to the
// edge of the patch if necessary
if (project_to_facet(facet, this_point, areacoord, close_point,
outside_facet, compareTol) != MB_SUCCESS)
{
return MB_FAILURE;
}
//if (closest_point_ptr)
//*closest_point_ptr = close_point;
}
// keep track of the minimum distance
double dist = (close_point - this_point).length();//close_point.distance_between(this_point);
if ((best_outside_facet == outside_facet && dist < mindist)
|| (best_outside_facet && !outside_facet && (dist < big_dist
|| best_facet == 0L)))
{
mindist = dist;
best_point = close_point;
best_facet = facet;
best_areacoord = areacoord;
best_outside_facet = outside_facet;
if (dist < compareTol)
{
break;
}
big_dist = 10.0 * mindist;
}
//facet->marked(1);
//used_facet_list.append(facet);
}
if (normal_ptr)
{
CartVect normal;
if (eval_bezier_patch_normal(best_facet, best_areacoord, normal)
!= MB_SUCCESS)
{
return MB_FAILURE;
}
*normal_ptr = normal;
}
if (closest_point_ptr)
{
*closest_point_ptr = best_point;
}
outside = best_outside_facet;
lastFacet = best_facet;
return MB_SUCCESS;
//end copy
}
| ErrorCode moab::SmoothFace::project_to_facets_main | ( | CartVect & | this_point, |
| bool | trim, | ||
| bool & | outside, | ||
| CartVect * | closest_point_ptr = NULL, |
||
| CartVect * | normal_ptr = NULL |
||
| ) |
Definition at line 1105 of file SmoothFace.cpp.
{
// if there are a lot of facets on this surface - use the OBB search first
// to narrow the selection
_evaluationsCounter++;
double tolerance = 1.e-5;
std::vector<EntityHandle> facets_out;
// we will start with a list of facets anyway, the best among them wins
ErrorCode rval = _my_geomTopoTool->obb_tree()->closest_to_location(
(double*) &this_point, _obb_root, tolerance, facets_out);
int interpOrder = 4;
double compareTol = 1.e-5;
EntityHandle lastFacet = facets_out.front();
rval = project_to_facets(facets_out, lastFacet, interpOrder, compareTol,
this_point, trim, outside, closest_point_ptr, normal_ptr);
return rval;
}
| ErrorCode moab::SmoothFace::project_to_patch | ( | EntityHandle | facet, |
| CartVect & | ac, | ||
| CartVect & | pt, | ||
| CartVect & | eval_pt, | ||
| CartVect * | eval_norm, | ||
| bool & | outside, | ||
| double | compare_tol, | ||
| int | edge_id | ||
| ) |
Definition at line 1224 of file SmoothFace.cpp.
{
ErrorCode status = MB_SUCCESS;
// see if we are at a vertex
#define INCR 0.01
const double tol = compare_tol;
if (is_at_vertex(facet, pt, ac, compare_tol, eval_pt, eval_norm))
{
outside = false;
return MB_SUCCESS;
}
// check if the start ac is inside the patch -if not, then move it there
int nout = 0;
const double atol = 0.001;
if (move_ac_inside(ac, atol))
nout++;
int diverge = 0;
int iter = 0;
CartVect newpt;
eval_bezier_patch(facet, ac, newpt);
CartVect move = pt - newpt;
double lastdist = move.length();
double bestdist = lastdist;
CartVect bestac = ac;
CartVect bestpt = newpt;
CartVect bestnorm(0, 0, 0);
// If we are already close enough, then return now
if (lastdist <= tol && !eval_norm && nout == 0)
{
eval_pt = pt;
outside = false;
return status;
}
double ratio, mag, umove, vmove, det, distnew, movedist;
CartVect lastpt = newpt;
CartVect lastac = ac;
CartVect norm;
CartVect xpt, ypt, zpt, xac, yac, zac, xvec, yvec, zvec;
CartVect du, dv, newac;
bool done = false;
while (!done)
{
// We will be locating the projected point within the u,v,w coordinate
// system of the triangular bezier patch. Since u+v+w=1, the components
// are not linearly independent. We will choose only two of the
// coordinates to use and compute the third.
int system;
if (lastac[0] >= lastac[1] && lastac[0] >= lastac[2])
{
system = 0;
}
else if (lastac[1] >= lastac[2])
{
system = 1;
}
else
{
system = 2;
}
// compute the surface derivatives with respect to each
// of the barycentric coordinates
if (system == 1 || system == 2)
{
xac[0] = lastac[0] + INCR; // xac.x( lastac.x() + INCR );
if (lastac[1] + lastac[2] == 0.0)
return MB_FAILURE;
ratio = lastac[2] / (lastac[1] + lastac[2]);//ratio = lastac.z() / (lastac.y() + lastac.z());
xac[1] = (1.0 - xac[0]) * (1.0 - ratio);//xac.y( (1.0 - xac.x()) * (1.0 - ratio) );
xac[2] = 1.0 - xac[0] - xac[1]; // xac.z( 1.0 - xac.x() - xac.y() );
eval_bezier_patch(facet, xac, xpt);
xvec = xpt - lastpt;
xvec /= INCR;
}
if (system == 0 || system == 2)
{
yac[1] = (lastac[1] + INCR);//yac.y( lastac.y() + INCR );
if (lastac[0] + lastac[2] == 0.0)//if (lastac.x() + lastac.z() == 0.0)
return MB_FAILURE;
ratio = lastac[2] / (lastac[0] + lastac[2]);//ratio = lastac.z() / (lastac.x() + lastac.z());
yac[0] = ((1.0 - yac[1]) * (1.0 - ratio));//yac.x( (1.0 - yac.y()) * (1.0 - ratio) );
yac[2] = (1.0 - yac[0] - yac[1]);//yac.z( 1.0 - yac.x() - yac.y() );
eval_bezier_patch(facet, yac, ypt);
yvec = ypt - lastpt;
yvec /= INCR;
}
if (system == 0 || system == 1)
{
zac[2] = (lastac[2] + INCR);//zac.z( lastac.z() + INCR );
if (lastac[0] + lastac[1] == 0.0)//if (lastac.x() + lastac.y() == 0.0)
return MB_FAILURE;
ratio = lastac[1] / (lastac[0] + lastac[1]);//ratio = lastac.y() / (lastac.x() + lastac.y());
zac[0] = ((1.0 - zac[2]) * (1.0 - ratio));//zac.x( (1.0 - zac.z()) * (1.0 - ratio) );
zac[1] = (1.0 - zac[0] - zac[2]);//zac.y( 1.0 - zac.x() - zac.z() );
eval_bezier_patch(facet, zac, zpt);
zvec = zpt - lastpt;
zvec /= INCR;
}
// compute the surface normal
switch (system) {
case 0:
du = yvec;
dv = zvec;
break;
case 1:
du = zvec;
dv = xvec;
break;
case 2:
du = xvec;
dv = yvec;
break;
}
norm = du * dv;
mag = norm.length();
if (mag < DBL_EPSILON)
{
return MB_FAILURE;
// do something else here (it is likely a flat triangle -
// so try evaluating just an edge of the bezier patch)
}
norm /= mag;
if (iter == 0)
bestnorm = norm;
// project the move vector to the tangent plane
move = (norm * move) * norm;
// compute an equivalent u-v-w vector
CartVect absnorm(fabs(norm[0]), fabs(norm[1]), fabs(norm[2]));
if (absnorm[2] >= absnorm[1] && absnorm[2] >= absnorm[0])
{
det = du[0] * dv[1] - dv[0] * du[1];
if (fabs(det) <= DBL_EPSILON)
{
return MB_FAILURE; // do something else here
}
umove = (move[0] * dv[1] - dv[0] * move[1]) / det;
vmove = (du[0] * move[1] - move[0] * du[1]) / det;
}
else if (absnorm[1] >= absnorm[2] && absnorm[1] >= absnorm[0])
{
det = du[0] * dv[2] - dv[0] * du[2];
if (fabs(det) <= DBL_EPSILON)
{
return MB_FAILURE;
}
umove = (move[0] * dv[2] - dv[0] * move[2]) / det;
vmove = (du[0] * move[2] - move[0] * du[2]) / det;
}
else
{
det = du[1] * dv[2] - dv[1] * du[2];
if (fabs(det) <= DBL_EPSILON)
{
return MB_FAILURE;
}
umove = (move[1] * dv[2] - dv[1] * move[2]) / det;
vmove = (du[1] * move[2] - move[1] * du[2]) / det;
}
/* === compute the new u-v coords and evaluate surface at new location */
switch (system) {
case 0:
newac[1] = (lastac[1] + umove);//newac.y( lastac.y() + umove );
newac[2] = (lastac[2] + vmove);//newac.z( lastac.z() + vmove );
newac[0] = (1.0 - newac[1] - newac[2]);//newac.x( 1.0 - newac.y() - newac.z() );
break;
case 1:
newac[2] = (lastac[2] + umove);//newac.z( lastac.z() + umove );
newac[0] = (lastac[0] + vmove);//newac.x( lastac.x() + vmove );
newac[1] = (1.0 - newac[2] - newac[0]);//newac.y( 1.0 - newac.z() - newac.x() );
break;
case 2:
newac[0] = (lastac[0] + umove);//newac.x( lastac.x() + umove );
newac[1] = (lastac[1] + vmove);//newac.y( lastac.y() + vmove );
newac[2] = (1.0 - newac[0] - newac[1]);//newac.z( 1.0 - newac.x() - newac.y() );
break;
}
// Keep it inside the patch
if (newac[0] >= -atol && newac[1] >= -atol && newac[2] >= -atol)
{
nout = 0;
}
else
{
if (move_ac_inside(newac, atol))
nout++;
}
// Evaluate at the new location
if (edge_id != -1)
ac_at_edge(newac, newac, edge_id); // move to edge first
eval_bezier_patch(facet, newac, newpt);
// Check for convergence
distnew = (pt - newpt).length();//pt.distance_between(newpt);
move = newpt - lastpt;
movedist = move.length();
if (movedist < tol || distnew < tol)
{
done = true;
if (distnew < bestdist)
{
bestdist = distnew;
bestac = newac;
bestpt = newpt;
bestnorm = norm;
}
}
else
{
// don't allow more than 30 iterations
iter++;
if (iter > 30)
{
//if (movedist > tol * 100.0) nout=1;
done = true;
}
// Check for divergence - don't allow more than 5 divergent
// iterations
if (distnew > lastdist)
{
diverge++;
if (diverge > 10)
{
done = true;
//if (movedist > tol * 100.0) nout=1;
}
}
// Check if we are continuing to project outside the facet.
// If so, then stop now
if (nout > 3)
{
done = true;
}
// set up for next iteration
if (!done)
{
if (distnew < bestdist)
{
bestdist = distnew;
bestac = newac;
bestpt = newpt;
bestnorm = norm;
}
lastdist = distnew;
lastpt = newpt;
lastac = newac;
move = pt - lastpt;
}
}
}
eval_pt = bestpt;
if (eval_norm)
{
*eval_norm = bestnorm;
}
outside = (nout > 0) ? true : false;
ac = bestac;
return status;
}
| ErrorCode moab::SmoothFace::ray_intersection_correct | ( | EntityHandle | facet, |
| CartVect & | pt, | ||
| CartVect & | ray, | ||
| CartVect & | eval_pt, | ||
| double & | distance, | ||
| bool & | outside | ||
| ) |
Definition at line 1939 of file SmoothFace.cpp.
{
// find a point on the smooth surface
CartVect currentPoint = eval_pt;
int numIter = 0;
double improvement = 1.e20;
CartVect diff;
while (numIter++ < 5 && improvement > 0.01)
{
CartVect newPos;
bool trim = false;// is it needed?
outside = true;
CartVect closestPoint;
CartVect normal;
ErrorCode rval = project_to_facets_main(currentPoint, trim, outside,
&newPos, &normal);
if (MB_SUCCESS != rval) return rval;
assert(rval==MB_SUCCESS);
diff = newPos - currentPoint;
improvement = diff.length();
// ( pt + t * ray - closest ) % normal = 0;
// intersect tangent plane that goes through closest point with the direction
// t = normal%(closest-pt) / normal%ray;
double dot = normal % ray; // if it is 0, get out while we can
if (dot < 0.00001)
{
// bad convergence, get out, do not modify anything
return MB_SUCCESS;
}
double t = ((newPos - pt) % normal) / (dot);
currentPoint = pt + t * ray;
}
eval_pt = currentPoint;
diff = currentPoint - pt;
distance = diff.length();
return MB_SUCCESS;
}
Tag moab::SmoothFace::_edgeCtrlTag [private] |
Definition at line 205 of file SmoothFace.hpp.
Range moab::SmoothFace::_edges [private] |
Definition at line 174 of file SmoothFace.hpp.
long moab::SmoothFace::_evaluationsCounter [private] |
Definition at line 231 of file SmoothFace.hpp.
Tag moab::SmoothFace::_facetCtrlTag [private] |
Definition at line 211 of file SmoothFace.hpp.
Tag moab::SmoothFace::_facetEdgeCtrlTag [private] |
Definition at line 219 of file SmoothFace.hpp.
Tag moab::SmoothFace::_gradientTag [private] |
Definition at line 193 of file SmoothFace.hpp.
Tag moab::SmoothFace::_markTag [private] |
Definition at line 188 of file SmoothFace.hpp.
CartVect moab::SmoothFace::_maxim [private] |
Definition at line 171 of file SmoothFace.hpp.
Interface* moab::SmoothFace::_mb [private] |
Definition at line 226 of file SmoothFace.hpp.
CartVect moab::SmoothFace::_minim [private] |
Definition at line 170 of file SmoothFace.hpp.
GeomTopoTool* moab::SmoothFace::_my_geomTopoTool [private] |
Definition at line 228 of file SmoothFace.hpp.
Range moab::SmoothFace::_nodes [private] |
Definition at line 175 of file SmoothFace.hpp.
EntityHandle moab::SmoothFace::_obb_root [private] |
Definition at line 229 of file SmoothFace.hpp.
Tag moab::SmoothFace::_planeTag [private] |
Definition at line 224 of file SmoothFace.hpp.
EntityHandle moab::SmoothFace::_set [private] |
Definition at line 227 of file SmoothFace.hpp.
Tag moab::SmoothFace::_tangentsTag [private] |
Definition at line 198 of file SmoothFace.hpp.
Range moab::SmoothFace::_triangles [private] |
Definition at line 173 of file SmoothFace.hpp.
1.7.6.1