[MOAB-dev] commit/MOAB: 3 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Wed Jun 19 10:28:39 CDT 2013
3 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/23fc4da626cd/
Changeset: 23fc4da626cd
Branch: None
User: iulian07
Date: 2013-06-14 01:27:17
Summary: more checks for convexity, and different oriented angle formula; break the concave quad in 2 triangles; bail out if more than 2 angles are bigger than 180 degrees; it is a signal for really bad quads, self intersecting, probably
Affected #: 3 files
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 61d1bf0..610c39d 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -12,6 +12,11 @@
#include "ElemUtil.hpp"
#include "moab/MergeMesh.hpp"
+// this is for sstream
+#include <sstream>
+
+#include <queue>
+
namespace moab {
// vec utilities that could be common between quads on a plane or sphere
double dist2(double * a, double * b)
@@ -636,7 +641,7 @@ double oriented_spherical_angle(double * A, double * B, double * C)
CartVect a(A), b(B), c(C);
CartVect normalOAB = a * b;
CartVect normalOCB = c * b;
- CartVect orient = (b-a)*(c-a);
+ CartVect orient = (c-b)*(a-b);
double ang = angle(normalOAB, normalOCB); // this is between 0 and M_PI
if (ang!=ang)
{
@@ -644,8 +649,8 @@ double oriented_spherical_angle(double * A, double * B, double * C)
std::cout << a << " " << b << " " << c <<"\n";
std::cout << ang << "\n";
}
- if (orient%a < 0)
- return (2*M_PI-ang);// the other angle
+ if (orient%b < 0)
+ return (2*M_PI-ang);// the other angle, supplement
return ang;
@@ -861,4 +866,133 @@ void departure_point_case1(CartVect & arrival_point, double t, double delta_t, C
departure_point = spherical_to_cart(sph_dep);
return;
}
+// break the nonconvex quads into triangles; remove the quad from the set? yes.
+// maybe radius is not needed;
+//
+ErrorCode enforce_convexity(Interface * mb, EntityHandle lset)
+{
+ // look at each quad; compute all 4 angles; if one is reflex, break along that diagonal
+ // replace it with 2 triangles, and remove from set;
+ // it should work for all polygons / tested first for case 1, with dt 0.5 (too much deformation)
+ // get all entities of dimension 2
+ // then get the connectivity, etc
+ Range inputRange;
+ ErrorCode rval = mb->get_entities_by_dimension(lset, 2, inputRange);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ std::vector<double> coords;
+ coords.resize(3*10); // at most 10 vertices per polygon
+ // we should create a queue with new polygons that need processing for reflex angles
+ // (obtuse)
+ std::queue<EntityHandle> newPolys;
+ int brokenQuads=0;
+ for (Range::iterator eit = inputRange.begin(); eit != inputRange.end()
+ || !newPolys.empty() ; eit++)
+ {
+ EntityHandle eh;
+ if (eit != inputRange.end())
+ eh = *eit;
+ else
+ {
+ eh = newPolys.front();
+ newPolys.pop();
+ }
+ // get the nodes, then the coordinates
+ const EntityHandle * verts;
+ int num_nodes;
+ rval = mb->get_connectivity(eh, verts, num_nodes);
+ if (MB_SUCCESS != rval)
+ return rval;
+ coords.resize(3 * num_nodes);
+ if (num_nodes < 4)
+ continue; // if already triangles, don't bother
+ // get coordinates
+ rval = mb->get_coords(verts, num_nodes, &coords[0]);
+ if (MB_SUCCESS != rval)
+ return rval;
+ // compute each angle
+ bool alreadyBroken = false;
+
+ for (int i=0; i<num_nodes; i++)
+ {
+ double * A = &coords[3*i];
+ double * B = &coords[3*((i+1)%num_nodes)];
+ double * C = &coords[3*((i+2)%num_nodes)];
+ double angle = oriented_spherical_angle(A, B, C);
+ if (angle-M_PI > 0.) // even almost reflex is bad; break it!
+ {
+ if (alreadyBroken)
+ {
+ mb->list_entities(&eh, 1);
+ mb->list_entities(verts, num_nodes);
+ double * D = &coords[3*((i+3)%num_nodes)];
+ std::cout<< "ABC: " << angle << " \n";
+ std::cout<< "BCD: " << oriented_spherical_angle( B, C, D) << " \n";
+ std::cout<< "CDA: " << oriented_spherical_angle( C, D, A) << " \n";
+ std::cout<< "DAB: " << oriented_spherical_angle( D, A, B)<< " \n";
+ std::cout << " this quad has at least 2 angles > 180, it has serious issues\n";
+
+ return MB_FAILURE;
+ }
+ // the bad angle is at i+1;
+ // create 1 triangle and one polygon; add the polygon to the input range, so
+ // it will be processed too
+ // also, add both to the set :) and remove the original polygon from the set
+ // break the next triangle, even though not optimal
+ // so create the triangle i+1, i+2, i+3; remove i+2 from original list
+ // even though not optimal in general, it is good enough.
+ EntityHandle conn3[3]={ verts[ (i+1)%num_nodes],
+ verts[ (i+2)%num_nodes],
+ verts[ (i+3)%num_nodes] };
+ // create a polygon with num_nodes-1 vertices, and connectivity
+ // verts[i+1], verts[i+3], (all except i+2)
+ std::vector<EntityHandle> conn(num_nodes-1);
+ for (int j=1; j<num_nodes; j++)
+ {
+ conn[j-1]=verts[(i+j+2)%num_nodes];
+ }
+ EntityHandle newElement;
+ rval = mb->create_element(MBTRI, conn3, 3, newElement);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = mb->add_entities(lset, &newElement, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (num_nodes == 4)
+ {
+ // create another triangle
+ rval = mb->create_element(MBTRI, &conn[0], 3, newElement);
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+ else
+ {
+ // create another polygon, and add it to the inputRange
+ rval = mb->create_element(MBPOLYGON, &conn[0], num_nodes-1, newElement);
+ if (MB_SUCCESS != rval)
+ return rval;
+ newPolys.push(newElement); // because it has less number of edges, the
+ // reverse should work to find it.
+ }
+ rval = mb->add_entities(lset, &newElement, 1);
+ if (MB_SUCCESS != rval)
+ return rval;
+ mb->remove_entities(lset, &eh, 1);
+ brokenQuads++;
+ /*std::cout<<"remove: " ;
+ mb->list_entities(&eh, 1);
+
+ std::stringstream fff;
+ fff << "file0" << brokenQuads<< ".vtk";
+ mb->write_file(fff.str().c_str(), 0, 0, &lset, 1);*/
+ alreadyBroken=true; // get out of the loop, element is broken
+ }
+ }
+ }
+ std::cout << brokenQuads << " quads were decomposed in triangles \n";
+ return MB_SUCCESS;
+}
+
} //namespace moab
diff --git a/tools/mbcslam/CslamUtils.hpp b/tools/mbcslam/CslamUtils.hpp
index e244138..1373a4f 100644
--- a/tools/mbcslam/CslamUtils.hpp
+++ b/tools/mbcslam/CslamUtils.hpp
@@ -122,5 +122,9 @@ double distance_on_great_circle(CartVect & p1, CartVect & p2);
void departure_point_case1(CartVect & arrival_point, double t, double delta_t, CartVect & departure_point);
+// break the nonconvex quads into triangles; remove the quad from the set? yes.
+// maybe radius is not needed;
+//
+ErrorCode enforce_convexity(Interface * mb, EntityHandle set);
}
#endif /* CSLAMUTILS_HPP_ */
diff --git a/tools/mbcslam/case1_test.cpp b/tools/mbcslam/case1_test.cpp
index 045cc1e..c15e3c7 100644
--- a/tools/mbcslam/case1_test.cpp
+++ b/tools/mbcslam/case1_test.cpp
@@ -166,6 +166,13 @@ int main(int argc, char **argv)
rval = manufacture_lagrange_mesh_on_sphere(&mb, euler_set, lagrange_set);
if (MB_SUCCESS != rval)
return 1;
+ rval = mb.write_file("lagrIni.h5m", 0, 0, &lagrange_set, 1);
+ if (MB_SUCCESS != rval)
+ std::cout << "can't write lagr set\n";
+
+ rval = enforce_convexity(&mb, lagrange_set);
+ if (MB_SUCCESS != rval)
+ return 1;
rval = mb.write_file("lagr.h5m", 0, 0, &lagrange_set, 1);
if (MB_SUCCESS != rval)
@@ -174,6 +181,7 @@ int main(int argc, char **argv)
Intx2MeshOnSphere worker(&mb);
double radius = CubeSide / 2 * sqrt(3.); // input
+
worker.SetRadius(radius);
worker.SetEntityType(MBQUAD);
https://bitbucket.org/fathomteam/moab/commits/a58f50c32a8c/
Changeset: a58f50c32a8c
Branch: None
User: iulian07
Date: 2013-06-19 05:48:09
Summary: change the intersection algorithms to accept any polygons
(right now, use case1_test with dt=0.43; it will have 6 quads concave, that
need to be decomposed in triangles, after manufacturing)
remove "set entity type" method, it is not be needed anymore
it will start with manifold sets that cover the same domain
change accordingly CslamUtils methods; edge intersect will take
convex polygons from now on.
Affected #: 14 files
diff --git a/test/parallel/par_intx_sph.cpp b/test/parallel/par_intx_sph.cpp
index de30f8b..98d0f7a 100644
--- a/test/parallel/par_intx_sph.cpp
+++ b/test/parallel/par_intx_sph.cpp
@@ -180,7 +180,7 @@ void test_intx_in_parallel()
double radius= CubeSide/2 * sqrt(3.) ; // input
worker.SetRadius(radius);
worker.set_box_error(EPS1);//
- worker.SetEntityType(MBQUAD);
+ //worker.SetEntityType(MBQUAD);
worker.SetErrorTolerance(radius*1.e-8);
worker.locate_departure_points(euler_set);
@@ -247,7 +247,7 @@ void test_intx_in_parallel_elem_based()
double radius= CubeSide/2 * sqrt(3.) ; // input
worker.SetRadius(radius);
worker.set_box_error(EPS1);//
- worker.SetEntityType(MBQUAD);
+ //worker.SetEntityType(MBQUAD);
worker.SetErrorTolerance(radius*1.e-8);
std::cout << "error tolerance epsilon_1="<< radius*1.e-8 << "\n";
diff --git a/tools/mbcslam/CslamUtils.cpp b/tools/mbcslam/CslamUtils.cpp
index 610c39d..06c3e18 100644
--- a/tools/mbcslam/CslamUtils.cpp
+++ b/tools/mbcslam/CslamUtils.cpp
@@ -29,24 +29,26 @@ double area2D(double *a, double *b, double *c)
// (b-a)x(c-a) / 2
return ((b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0])) / 2;
}
-int borderPointsOfXinY2(double * X, double * Y, int nsides, double * P, int side[4])
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES])
{
// 2 triangles, 3 corners, is the corner of X in Y?
// Y must have a positive area
/*
*/
int extraPoint = 0;
- for (int i = 0; i < nsides; i++)
+ for (int i = 0; i < nX; i++)
{
- // compute twice the area of all 4 triangles formed by a side of Y and a corner of X; if one is negative, stop
+ // compute twice the area of all nY triangles formed by a side of Y and a corner of X; if one is negative, stop
+ // (negative means it is outside; X and Y are all oriented such that they are positive oriented;
+ // if one area is negative, it means it is outside the convex region, for sure)
double * A = X + 2 * i;
int inside = 1;
- for (int j = 0; j < nsides; j++)
+ for (int j = 0; j < nY; j++)
{
double * B = Y + 2 * j;
- int j1 = (j + 1) % nsides;
+ int j1 = (j + 1) % nY;
double * C = Y + 2 * j1; // no copy of data
double area2 = (B[0] - A[0]) * (C[1] - A[1])
@@ -59,7 +61,8 @@ int borderPointsOfXinY2(double * X, double * Y, int nsides, double * P, int side
}
if (inside)
{
- side[i] = 1;
+ side[i] = 1;// so vertex i of X is inside the convex region formed by Y
+ // so side has nX dimension (first array)
P[extraPoint * 2] = A[0];
P[extraPoint * 2 + 1] = A[1];
extraPoint++;
@@ -89,7 +92,8 @@ int SortAndRemoveDoubles2(double * P, int & nP, double epsilon_1)
}
c[0] /= nP;
c[1] /= nP;
- double angle[24]; // could be at most 24 points; much less usually
+ // how many
+ std::vector<double> angle(nP); // could be at most nP points
for (k = 0; k < nP; k++)
{
double x = P[2 * k] - c[0], y = P[2 * k + 1] - c[1];
@@ -113,7 +117,7 @@ int SortAndRemoveDoubles2(double * P, int & nP, double epsilon_1)
if (angle[k] > angle[k + 1])
{
sorted = 0;
- swap2(angle + k, angle + k + 1);
+ swap2(&angle[k], &angle[k+1]);
swap2(P + (2 * k), P + (2 * k + 2));
swap2(P + (2 * k + 1), P + (2 * k + 3));
}
@@ -151,8 +155,8 @@ int SortAndRemoveDoubles2(double * P, int & nP, double epsilon_1)
// the marks will show what edges of blue intersect the red
-int EdgeIntersections2(double * blue, double * red, int nsides, int markb[4], int markr[4],
- double * points, int & nPoints)
+int EdgeIntersections2(double * blue, int nsBlue, double * red, int nsRed,
+ int markb[MAXEDGES], int markr[MAXEDGES], double * points, int & nPoints)
{
/* EDGEINTERSECTIONS computes edge intersections of two elements
[P,n]=EdgeIntersections(X,Y) computes for the two given elements * red
@@ -163,30 +167,20 @@ int EdgeIntersections2(double * blue, double * red, int nsides, int markb[4], in
with blue are given.
*/
- // points is an array with 48 slots (24 * 2 doubles)
+ // points is an array with enough slots (24 * 2 doubles)
nPoints = 0;
- markb[0] = markb[1] = markb[2] = markb[3] = 0; // no neighbors of red involved yet
- markr[0] = markr[1] = markr[2] = markr[3] = 0;
- /*for i=1:3 % find all intersections of edges
- for j=1:3
- b=Y(:,j)-X(:,i);
- A=[X(:,mod(i,3)+1)-X(:,i) -Y(:,mod(j,3)+1)+Y(:,j)];
- if rank(A)==2 % edges not parallel
- r=A\b;
- if r(1)>=0 & r(1)<=1 & r(2)>=0 & r(2)<=1, % intersection found
- k=k+1; P(:,k)=X(:,i)+r(1)*(X(:,mod(i,3)+1)-X(:,i)); n(i)=1;
- end;
- end;
- end;
- end;*/
- for (int i = 0; i < nsides; i++)
- {
- for (int j = 0; j < nsides; j++)
+ for (int i=0; i<MAXEDGES; i++){
+ markb[i]=markr[i] = 0;
+ }
+
+ for (int i = 0; i < nsBlue; i++)
+ {
+ for (int j = 0; j < nsRed; j++)
{
double b[2];
double a[2][2]; // 2*2
- int iPlus1 = (i + 1) % nsides;
- int jPlus1 = (j + 1) % nsides;
+ int iPlus1 = (i + 1) % nsBlue;
+ int jPlus1 = (j + 1) % nsRed;
for (int k = 0; k < 2; k++)
{
b[k] = red[2 * j + k] - blue[2 * i + k];
diff --git a/tools/mbcslam/CslamUtils.hpp b/tools/mbcslam/CslamUtils.hpp
index 1373a4f..c32a264 100644
--- a/tools/mbcslam/CslamUtils.hpp
+++ b/tools/mbcslam/CslamUtils.hpp
@@ -11,16 +11,20 @@
#include "moab/Core.hpp"
#include "moab/Interface.hpp"
+// maximum number of edges on each convex polygon of interest
+#define MAXEDGES 10
+#define MAXEDGES2 20 // used for coordinates in plane
+
namespace moab
{
double dist2(double * a, double * b);
double area2D(double *a, double *b, double *c);
-int borderPointsOfXinY2(double * X, double * Y, int nsides, double * P, int side[4]);
+int borderPointsOfXinY2(double * X, int nX, double * Y, int nY, double * P, int side[MAXEDGES]);
int SortAndRemoveDoubles2(double * P, int & nP, double epsilon);
// the marks will show what edges of blue intersect the red
-int EdgeIntersections2(double * blue, double * red, int nsides, int markb[4], int markr[4],
- double * points, int & nPoints);
+int EdgeIntersections2(double * blue, int nsBlue, double * red, int nsRed,
+ int markb[MAXEDGES], int markr[MAXEDGES], double * points, int & nPoints);
// vec utils related to gnomonic projection on a sphere
diff --git a/tools/mbcslam/Intx2Mesh.cpp b/tools/mbcslam/Intx2Mesh.cpp
index e7f5f26..a2d52e6 100644
--- a/tools/mbcslam/Intx2Mesh.cpp
+++ b/tools/mbcslam/Intx2Mesh.cpp
@@ -48,14 +48,14 @@ void Intx2Mesh::createTags()
ERRORV(rval, "can't get red flag tag");
// assume that the edges are on the red triangles
- Range redQuads;
+ Range redElements;
//Range redEdges;
- rval = mb->get_entities_by_type(mbs2, type, redQuads, false);
- ERRORV(rval, "can't get ents by type");
+ rval = mb->get_entities_by_dimension(mbs2, 2, redElements, false); // so all tri, quad, poly
+ ERRORV(rval, "can't get ents by dimension");
// create red edges if they do not exist yet; so when they are looked upon, they are found
// this is the only call that is potentially NlogN, in the whole method
- rval = mb->get_adjacencies(redQuads, 1, true, RedEdges, Interface::UNION);
+ rval = mb->get_adjacencies(redElements, 1, true, RedEdges, Interface::UNION);
ERRORV(rval, "can't get adjacent red edges");
// now, create a map from each edge to a list of potential new nodes on a red edge
@@ -91,18 +91,17 @@ void Intx2Mesh::createTags()
// specify also desired set; we are interested only in neighbors in the set!
-// we should always get manifold mesh, each edge is adjacent to 2 quads
-ErrorCode Intx2Mesh::GetOrderedNeighbors(EntityHandle set, EntityHandle quad,
- EntityHandle neighbors[4])
+// we should always get manifold mesh, each edge is adjacent to 2 cell
+// maybe we should check that first, just in case
+ErrorCode Intx2Mesh::GetOrderedNeighbors(EntityHandle set, EntityHandle cell,
+ EntityHandle neighbors[MAXEDGES])
{
- int nsides = 3;
- if (type == MBQUAD)
- nsides = 4;
- // will get the 4 ordered neighbors;
- // first quad is for nodes 0, 1, second to 1, 2, third to 2, 3,
- int nnodes;
+ int nnodes = 3;
+ // will get the nnodes ordered neighbors;
+ // first cell is for nodes 0, 1, second to 1, 2, third to 2, 3, last to nnodes-1,
const EntityHandle * conn4;
- ErrorCode rval = mb->get_connectivity(quad, conn4, nnodes);
+ ErrorCode rval = mb->get_connectivity(cell, conn4, nnodes);
+ int nsides = nnodes; // just keep it for historical purposes; it is indeed nnodes
ERRORR(rval, "can't get connectivity on an element");
for (int i = 0; i < nsides; i++)
{
@@ -110,45 +109,44 @@ ErrorCode Intx2Mesh::GetOrderedNeighbors(EntityHandle set, EntityHandle quad,
v[0] = conn4[i];
v[1] = conn4[(i + 1) % nsides];
// get quads adjacent to vertices
- std::vector<EntityHandle> quads;
- std::vector<EntityHandle> quadsInSet;
- rval = mb->get_adjacencies(v, 2, 2, false, quads, Interface::INTERSECT);
+ std::vector<EntityHandle> cells;
+ std::vector<EntityHandle> cellsInSet;
+ rval = mb->get_adjacencies(v, 2, 2, false, cells, Interface::INTERSECT);
ERRORR(rval, "can't get adjacencies on 2 nodes");
- size_t siz = quads.size();
+ size_t siz = cells.size();
for (size_t j = 0; j < siz; j++)
- if (mb->contains_entities(set, &(quads[j]), 1))
- quadsInSet.push_back(quads[j]);
- siz = quadsInSet.size();
+ if (mb->contains_entities(set, &(cells[j]), 1))
+ cellsInSet.push_back(cells[j]);
+ siz = cellsInSet.size();
if (siz > 2)
{
std::cout << "non manifold mesh, error"
- << mb->list_entities(&(quadsInSet[0]), quadsInSet.size()) << "\n";
+ << mb->list_entities(&(cellsInSet[0]), cellsInSet.size()) << "\n";
return MB_FAILURE; // non-manifold
}
if (siz == 1)
{
// it must be the border,
neighbors[i] = 0; // we are guaranteed that ids are !=0; this is marking a border
+ // borders do not appear for a sphere in serial, but they do appear for
+ // parallel processing anyway
continue;
}
// here siz ==2, it is either the first or second
- if (quad == quads[0])
- neighbors[i] = quads[1];
+ if (cell == cellsInSet[0])
+ neighbors[i] = cellsInSet[1];
else
- neighbors[i] = quads[0];
+ neighbors[i] = cellsInSet[0];
}
return MB_SUCCESS;
}
// main interface; this will do the advancing front trick
-// some are triangles, some are quads
+// some are triangles, some are quads, some are polygons ...
ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
EntityHandle & outputSet)
{
- int nsides=3;
- if (type == MBQUAD)
- nsides=4;
ErrorCode rval;
mbs1 = mbset1; // set 1 is departure, and it is completely covering the euler set on proc
mbs2 = mbset2;
@@ -160,8 +158,8 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
Range rs1;
Range rs2;
- mb->get_entities_by_type(mbs1, type, rs1);
- mb->get_entities_by_type(mbs2, type, rs2);
+ mb->get_entities_by_dimension(mbs1, 2, rs1);
+ mb->get_entities_by_dimension(mbs2, 2, rs2);
while (!rs2.empty())
{
for (Range::iterator it = rs1.begin(); it != rs1.end(); it++)
@@ -173,11 +171,13 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
startRed = *it2;
double area = 0;
// if area is > 0 , we have intersections
- double P[48]; // max 8 intx points + 8 more in the polygon
- // the red quad is convex, always, while the blue can be concave
+ double P[10*MAXEDGES]; // max 8 intx points + 8 more in the polygon
+ //
int nP = 0;
- int nb[4], nr[4]; // sides 3 or 4? also, check boxes first
- computeIntersectionBetweenRedAndBlue(startRed, startBlue, P, nP, area, nb, nr, true);
+ int nb[MAXEDGES], nr[MAXEDGES]; // sides 3 or 4? also, check boxes first
+ int nsRed, nsBlue;
+ computeIntersectionBetweenRedAndBlue(startRed, startBlue, P, nP, area, nb, nr,
+ nsBlue, nsRed, true);
if (area > 0)
{
found = 1;
@@ -206,9 +206,11 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
{
// flags for the side : 0 means a blue quad not found on side
// a paired blue not found yet for the neighbors of red
- EntityHandle n[4] = { EntityHandle(0) };
+ EntityHandle n[MAXEDGES] = { EntityHandle(0) };
+ int nsidesRed; // will be initialized later, when we compute intersection
EntityHandle currentRed = redQueue.front();
+
redQueue.pop();
// for (k=0; k<m_numPos; k++)
// redFlag[k] = 0;
@@ -222,6 +224,12 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
ERRORR(rval, "can't set blue unused tag");
}
//rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused);
+ /*if (dbg_1)
+ {
+ std::cout << "reset blues: ";
+ mb->list_entities(toResetBlues);
+ }*/
+ //rval = mb2->tag_set_data(RedFlagTag, toResetReds, &unused);
if (dbg_1)
{
std::cout << "reset blues: ";
@@ -247,19 +255,32 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
//
EntityHandle blueT = localBlue.front();
localBlue.pop();
- double P[48], area; // area is in 2d, points are in 3d (on a sphere), back-projected
- int nP = 0; // intersection points (could include the vertices of initial quads)
- int nb[4] = { 0, 0, 0, 0 }; // means no intersection on the side (markers)
- int nr[4] = { 0, 0, 0, 0 }; // means no intersection on the side (markers)
- // nc [j] = 1 means that the side j (from j to j+1) of blue quad intersects the
- // red quad. A potential next quad is the red quad that is adjacent to this side
+ double P[10*MAXEDGES], area; //
+ int nP = 0;
+ int nb[MAXEDGES] = {0};
+ int nr[MAXEDGES] = {0};
+
+ int nsidesBlue; // nsides red is the same, nsidesBlue changes for next element
+ // area is in 2d, points are in 3d (on a sphere), back-projected
+ // intersection points (could include the vertices of initial quads)
+ // means no intersection on the side (markers)
+ // means no intersection on the side (markers)
+ // nc [j] = 1 means that the side j (from j to j+1) of blue poly intersects the
+ // red poly. A potential next poly is the red poly that is adjacent to this side
computeIntersectionBetweenRedAndBlue(/* red */currentRed, blueT, P, nP,
- area, nb, nr);
+ area, nb, nr, nsidesBlue, nsidesRed);
if (nP > 0)
{
+ if (dbg_1)
+ {
+ for (int k=0; k<3; k++)
+ {
+ std::cout << " nb, nr: " << k << " " << nb[k] << " " << nr[k] << "\n";
+ }
+ }
// intersection found: output P and original triangles if nP > 2
- EntityHandle neighbors[4];
+ EntityHandle neighbors[MAXEDGES];
rval = GetOrderedNeighbors(mbs1, blueT, neighbors);
if (rval != MB_SUCCESS)
{
@@ -269,7 +290,7 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
}
// add neighbors to the localBlue queue, if they are not marked
- for (int nn = 0; nn < nsides; nn++)
+ for (int nn = 0; nn < nsidesBlue; nn++)
{
EntityHandle neighbor = neighbors[nn];
if (neighbor > 0 && nb[nn]>0) // advance across blue boundary n
@@ -280,11 +301,17 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
if (status == 0)
{
localBlue.push(neighbor);
+ /*if (dbg_1)
+ {
+ std::cout << " local blue elem "
+ << mb->list_entities(&neighbor, 1) << "\n for red:"
+ << mb->list_entities(¤tRed, 1) << "\n";
+ }*/
if (dbg_1)
{
std::cout << " local blue elem " << mb->id_from_handle(neighbor)
- << " for red:" << mb->id_from_handle(currentRed) << "\n"
- << mb->list_entities(&neighbor, 1) << "\n";
+ << " for red:" << mb->id_from_handle(currentRed) << "\n";
+ mb->list_entities(&neighbor, 1);
}
rval = mb->tag_set_data(BlueFlagTag, &neighbor, 1, &used);
//redFlag[neighbor] = 1; // flag it to not be added anymore
@@ -297,9 +324,14 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
}
if (nP > 1) // this will also construct triangles/polygons in the new mesh, if needed
- findNodes(currentRed, blueT, P, nP);
+ findNodes(currentRed, nsidesRed, blueT, nsidesBlue, P, nP);
}
else if (dbg_1)
+ /*{
+ std::cout << " red, blue, do not intersect: "
+ << mb->list_entities(¤tRed, 1) << " "
+ << mb->list_entities(&blueT, 1) << "\n";
+ }*/
{
std::cout << " red, blue, do not intersect: "
<< mb->id_from_handle(currentRed) << " "
@@ -311,13 +343,13 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
rs2.erase(currentRed);
// also, look at its neighbors, and add to the seeds a next one
- EntityHandle redNeighbors[4];
+ EntityHandle redNeighbors[MAXEDGES];
rval = GetOrderedNeighbors(mbs2, currentRed, redNeighbors);
ERRORR(rval, "can't get neighbors");
if (dbg_1)
{
std::cout << "Next: neighbors for current red ";
- for (int kk = 0; kk < nsides; kk++)
+ for (int kk = 0; kk < nsidesRed; kk++)
{
if (redNeighbors[kk] > 0)
std::cout << mb->id_from_handle(redNeighbors[kk]) << " ";
@@ -326,7 +358,7 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
}
std::cout << std::endl;
}
- for (int j = 0; j < nsides; j++)
+ for (int j = 0; j < nsidesRed; j++)
{
EntityHandle redNeigh = redNeighbors[j];
unsigned char status = 1;
@@ -339,7 +371,7 @@ ErrorCode Intx2Mesh::intersect_meshes(EntityHandle mbset1, EntityHandle mbset2,
redQueue.push(redNeigh);
blueQueue.push(n[j]);
if (dbg_1)
- std::cout << "new quads pushed: blue, red:"
+ std::cout << "new polys pushed: blue, red:"
<< mb->id_from_handle(redNeigh) << " "
<< mb->id_from_handle(n[j]) << std::endl;
mb->tag_set_data(RedFlagTag, &redNeigh, 1, &used);
diff --git a/tools/mbcslam/Intx2Mesh.hpp b/tools/mbcslam/Intx2Mesh.hpp
index 1e4fe04..b0424d9 100644
--- a/tools/mbcslam/Intx2Mesh.hpp
+++ b/tools/mbcslam/Intx2Mesh.hpp
@@ -39,24 +39,31 @@ public:
ErrorCode intersect_meshes(EntityHandle mbs1, EntityHandle mbs2,
EntityHandle & outputSet);
- // mark could be 3 or 4, depending on type
- // this is pure abstract, this need s to be implemented by
+ // mark could be (3 or 4, depending on type: ) no, it could go to 10
+ // no, it will be MAXEDGES = 10
+ // this is pure abstract, this needs to be implemented by
// all derivations
+ // the max number of intersection points could be 2*MAXEDGES
+ // so P must be dimensioned to 4*MAXEDGES (2*2*MAXEDGES)
+ // so, if you intersect 2 convex polygons with MAXEDGES , you will get a convex polygon
+ // with 2*MAXEDGES, at most
+ // will also return the number of nodes of red and blue elements
virtual int computeIntersectionBetweenRedAndBlue(EntityHandle red,
EntityHandle blue, double * P, int & nP, double & area,
- int markb[4], int markr[4], bool check_boxes_first=false)=0;
+ int markb[MAXEDGES], int markr[MAXEDGES], int & nsidesBlue,
+ int & nsidesRed, bool check_boxes_first=false)=0;
// this is also abstract
- virtual int findNodes(EntityHandle red, EntityHandle blue,
+ virtual int findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue,
double * iP, int nP)=0;
virtual void createTags();
ErrorCode GetOrderedNeighbors(EntityHandle set, EntityHandle quad,
- EntityHandle neighbors[4]);
+ EntityHandle neighbors[MAXEDGES]);
void SetErrorTolerance(double eps) { epsilon_1=eps;}
- void SetEntityType (EntityType tp) { type=tp;}
+ //void SetEntityType (EntityType tp) { type=tp;}
// clean some memory allocated
void clean();
@@ -81,6 +88,8 @@ public:
void correct_polygon(EntityHandle * foundIds, int & nP);
ErrorCode correct_intersection_points_positions();
+
+ void enable_debug() {dbg_1=1;};
protected: // so it can be accessed in derived classes, InPlane and OnSphere
Interface * mb;
@@ -102,20 +111,21 @@ protected: // so it can be accessed in derived classes, InPlane and OnSphere
Tag blueParentTag;
Tag countTag;
- EntityType type;
+ //EntityType type; // this will be tri, quad or MBPOLYGON...
const EntityHandle * redConn;
const EntityHandle * blueConn;
- CartVect redCoords[4];
- CartVect blueCoords[4];
- double redQuad[8]; // these are in plane
- double blueQuad[8]; // these are in plane
+ CartVect redCoords[MAXEDGES];
+ CartVect blueCoords[MAXEDGES];
+ double redCoords2D[MAXEDGES2]; // these are in plane
+ double blueCoords2D[MAXEDGES2]; // these are in plane
std::ofstream mout_1[6]; // some debug files
int dbg_1;
// for each red edge, we keep a vector of extra nodes, coming from intersections
// use the index in RedEdges range, instead of a map, as before
// std::map<EntityHandle, std::vector<EntityHandle> *> extraNodesMap;
+ // so the extra nodes on each red edge are kept track of
std::vector<std::vector<EntityHandle> *> extraNodesVec;
double epsilon_1;
diff --git a/tools/mbcslam/Intx2MeshInPlane.cpp b/tools/mbcslam/Intx2MeshInPlane.cpp
index 718c99a..4978f91 100644
--- a/tools/mbcslam/Intx2MeshInPlane.cpp
+++ b/tools/mbcslam/Intx2MeshInPlane.cpp
@@ -17,20 +17,18 @@ Intx2MeshInPlane::~Intx2MeshInPlane() {
}
int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, EntityHandle blue,
- double * P, int & nP, double & area, int markb[4], int markr[4], bool check_boxes_first)
+ double * P, int & nP, double & area, int markb[MAXEDGES], int markr[MAXEDGES],
+ int & nsBlue, int & nsRed, bool check_boxes_first)
{
- // the points will be at most 9; they will describe a convex patch, after the points will be ordered and
+ // the points will be at most ?; they will describe a convex patch, after the points will be ordered and
// collapsed (eliminate doubles)
// the area is not really required
- int nsides = 3;
- if (type == MBQUAD)
- nsides = 4;
+
int num_nodes;
ErrorCode rval = mb->get_connectivity(red, redConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != nsides)
- return 1;
+ nsRed = num_nodes;
//CartVect coords[4];
rval = mb->get_coords(redConn, num_nodes, &(redCoords[0][0]));
@@ -39,8 +37,9 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
rval = mb->get_connectivity(blue, blueConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != nsides)
+ if (MB_SUCCESS != rval)
return 1;
+ nsBlue = num_nodes;
rval = mb->get_coords(blueConn, num_nodes, &(blueCoords[0][0]));
if (MB_SUCCESS != rval)
return 1;
@@ -48,12 +47,12 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
if (dbg_1)
{
std::cout << "red " << mb->id_from_handle(red) << "\n";
- for (int j = 0; j < num_nodes; j++)
+ for (int j = 0; j < nsRed; j++)
{
std::cout << redCoords[j] << "\n";
}
std::cout << "blue " << mb->id_from_handle(blue) << "\n";
- for (int j = 0; j < num_nodes; j++)
+ for (int j = 0; j < nsBlue; j++)
{
std::cout << blueCoords[j] << "\n";
}
@@ -65,80 +64,112 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
if (check_boxes_first)
{
// look at the boxes formed with vertices; if they are far away, return false early
- if (!GeomUtil::bounding_boxes_overlap(redCoords, num_nodes, blueCoords, num_nodes, box_error))
+ if (!GeomUtil::bounding_boxes_overlap(redCoords, nsRed, blueCoords, nsBlue, box_error))
return 0; // no error, but no intersection, decide early to get out
}
- for (int j = 0; j < nsides; j++)
+ for (int j = 0; j < nsRed; j++)
{
// populate coords in the plane for intersection
// they should be oriented correctly, positively
- redQuad[2 * j]=redCoords[j][0]; // x coordinate,
- redQuad[2 * j + 1] = redCoords[j][1]; // y coordinate
- blueQuad[2 * j]=blueCoords[j][0]; // x coordinate,
- blueQuad[2 * j + 1] = blueCoords[j][1]; // y coordinate
+ redCoords2D[2 * j]=redCoords[j][0]; // x coordinate,
+ redCoords2D[2 * j + 1] = redCoords[j][1]; // y coordinate
}
- if (dbg_1)
+ for (int j=0; j<nsBlue; j++)
{
- //std::cout << "gnomonic plane: " << plane << "\n";
- std::cout << " red blue\n";
- for (int j = 0; j < 4; j++)
- {
- std::cout << redQuad[2 * j] << " " << redQuad[2 * j + 1] << " "
- << blueQuad[2 * j] << " " << blueQuad[2 * j + 1] << "\n";
- }
+ blueCoords2D[2 * j]=blueCoords[j][0]; // x coordinate,
+ blueCoords2D[2 * j + 1] = blueCoords[j][1]; // y coordinate
}
+ if (dbg_1)
+ {
+ //std::cout << "gnomonic plane: " << plane << "\n";
+ std::cout << " red \n";
+ for (int j = 0; j < nsRed; j++)
+ {
+ std::cout << redCoords2D[2 * j] << " " << redCoords2D[2 * j + 1] << "\n ";
+ }
+ std::cout << " blue\n";
+ for (int j = 0; j < nsBlue; j++)
+ {
+ std::cout << blueCoords2D[2 * j] << " " << blueCoords2D[2 * j + 1] << "\n";
+ }
+ }
- nP = 0; // number of intersection points
- int ret = EdgeIntersections2(blueQuad, redQuad, nsides, markb,
- markr, P, nP);
- if (ret != 0)
- return 1;// some unforeseen error
-
- // start copy
- int side[4] = { 0 };
- int extraPoints = borderPointsOfXinY2(blueQuad, redQuad, nsides,
- &(P[2 * nP]), side);
- if (extraPoints >= 1)
- {
- for (int k = 0; k < nsides; k++)
- {
- if (side[k])
- {
- markb[k] = 1;
- markb[(k + nsides-1) % nsides] = 1; // it is the previous edge, actually, but instead of doing -1, it is
- // better to do modulo +3 (modulo 4)
- // null side b for next call
- side[k]=0;
- }
- }
- }
- nP += extraPoints;
+ int ret = EdgeIntersections2(blueCoords2D, nsBlue, redCoords2D, nsRed, markb, markr, P, nP);
+ if (ret != 0)
+ return 1; // some unforeseen error
+ if (dbg_1)
+ {
+ for (int k=0; k<3; k++)
+ {
+ std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
+ }
+ }
- extraPoints = borderPointsOfXinY2(redQuad, blueQuad, nsides, &(P[2 * nP]), side);
- if (extraPoints >= 1)
- {
- for (int k = 0; k < 4; k++)
- {
- if (side[k])
- {
- markr[k] = 1;
- markr[(k + nsides-1) % nsides] = 1; // it is the previous edge, actually, but instead of doing -1, it is
- // better to do modulo +3 (modulo 4)
- // null side b for next call
- }
- }
- }
- nP += extraPoints;
+ int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
+ if (extraPoints >= 1)
+ {
+ for (int k = 0; k < nsBlue; k++)
+ {
+ if (side[k])
+ {
+ // this means that vertex k of blue is inside convex red; mark edges k-1 and k in blue,
+ // as being "intersected" by red; (even though they might not be intersected by other edges,
+ // the fact that their apex is inside, is good enough)
+ markb[k] = 1;
+ markb[(k + nsBlue-1) % nsBlue] = 1; // it is the previous edge, actually, but instead of doing -1, it is
+ // better to do modulo +3 (modulo 4)
+ // null side b for next call
+ side[k]=0;
+ }
+ }
+ }
+ if (dbg_1)
+ {
+ for (int k=0; k<3; k++)
+ {
+ std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
+ }
+ }
+ nP += extraPoints;
- //
- SortAndRemoveDoubles2(P, nP, epsilon_1); // nP should be at most 6 in the end
- // if there are more than 3 points, some area will be positive
- area = 0.;
- if (nP >= 3) {
- for (int k = 1; k < nP - 1; k++)
- area += area2D(P, &P[2 * k], &P[2 * k + 2]);
- }
- return 0; // no error
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
+ if (extraPoints >= 1)
+ {
+ for (int k = 0; k < nsRed; k++)
+ {
+ if (side[k])
+ {
+ // this is to mark that red edges k-1 and k are intersecting blue
+ markr[k] = 1;
+ markr[(k + nsRed-1) % nsRed] = 1; // it is the previous edge, actually, but instead of doing -1, it is
+ // better to do modulo +3 (modulo 4)
+ // null side b for next call
+ }
+ }
+ }
+ if (dbg_1)
+ {
+ for (int k=0; k<3; k++)
+ {
+ std::cout << " markb, markr: " << k << " " << markb[k] << " " << markr[k] << "\n";
+ }
+ }
+ nP += extraPoints;
+
+ // now sort and orient the points in P, such that they are forming a convex polygon
+ // this will be the foundation of our new mesh
+ // this works if the polygons are convex
+ SortAndRemoveDoubles2(P, nP, epsilon_1); // nP should be at most 8 in the end ?
+ // if there are more than 3 points, some area will be positive
+
+ if (nP >= 3)
+ {
+ for (int k = 1; k < nP - 1; k++)
+ area += area2D(P, &P[2 * k], &P[2 * k + 2]);
+ }
+
+ return 0; // no error
}
// this method will also construct the triangles/polygons in the new mesh
@@ -146,8 +177,11 @@ int Intx2MeshInPlane::computeIntersectionBetweenRedAndBlue(EntityHandle red, Ent
// also, we could just create new vertices every time, and merge only in the end;
// could be too expensive, and the tolerance for merging could be an
// interesting topic
-int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP, int nP)
+int Intx2MeshInPlane::findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue,
+ double * iP, int nP)
{
+ // except for gnomonic projection, everything is the same as spherical intx
+ // start copy
// first of all, check against red and blue vertices
//
if (dbg_1)
@@ -159,16 +193,13 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
}
- int nsides = 3;
- if (type == MBQUAD)
- nsides = 4;
// get the edges for the red triangle; the extra points will be on those edges, saved as
// lists (unordered)
- EntityHandle redEdges[4];
+ std::vector<EntityHandle> redEdges(nsRed);//
int i = 0;
- for (i = 0; i < nsides; i++)
+ for (i = 0; i < nsRed; i++)
{
- EntityHandle v[2] = { redConn[i], redConn[(i + 1) % nsides] };
+ EntityHandle v[2] = { redConn[i], redConn[(i + 1) % nsRed] };
std::vector<EntityHandle> adj_entities;
ErrorCode rval = mb->get_adjacencies(v, 2, 1, false, adj_entities,
Interface::INTERSECT);
@@ -183,18 +214,17 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
for (i = 0; i < nP; i++)
{
double * pp = &iP[2 * i]; // iP+2*i
- // project the point back on the sphere
- CartVect pos(pp[0], pp[1], 0.); // constructor with 3 doubles
- //reverse_gnomonic_projection(pp[0], pp[1], R, plane, pos);
+ //
+ CartVect pos(pp[0], pp[1], 0.);
int found = 0;
// first, are they on vertices from red or blue?
// priority is the red mesh (mb2?)
int j = 0;
EntityHandle outNode = (EntityHandle) 0;
- for (j = 0; j < nsides && !found; j++)
+ for (j = 0; j < nsRed && !found; j++)
{
//int node = redTri.v[j];
- double d2 = dist2(pp, &redQuad[2 * j]);
+ double d2 = dist2(pp, &redCoords2D[2 * j]);
if (d2 < epsilon_1)
{
@@ -202,15 +232,15 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
found = 1;
if (dbg_1)
std::cout << " red node j:" << j << " id:"
- << mb->id_from_handle(redConn[j]) << " 2d coords:" << redCoords[2 * j] << " "
- << redCoords[2 * j + 1] << " d2: " << d2 << " \n";
+ << mb->id_from_handle(redConn[j]) << " 2d coords:" << redCoords2D[2 * j] << " "
+ << redCoords2D[2 * j + 1] << " d2: " << d2 << " \n";
}
}
- for (j = 0; j < nsides && !found; j++)
+ for (j = 0; j < nsBlue && !found; j++)
{
//int node = blueTri.v[j];
- double d2 = dist2(pp, &blueQuad[2 * j]);
+ double d2 = dist2(pp, &blueCoords2D[2 * j]);
if (d2 < epsilon_1)
{
// suspect is blueConn[j] corresponding in mbOut
@@ -225,17 +255,17 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
}
if (!found)
{
- // find the edge it belongs, first, on the red quad
+ // find the edge it belongs, first, on the red element
//
- for (j = 0; j < nsides; j++)
+ for (j = 0; j < nsRed; j++)
{
- int j1 = (j + 1) % nsides;
- double area = area2D(&redQuad[2 * j], &redQuad[2 * j1], pp);
+ int j1 = (j + 1) % nsRed;
+ double area = area2D(&redCoords2D[2 * j], &redCoords2D[2 * j1], pp);
if (dbg_1)
std::cout << " edge " << j << ": "
<< mb->id_from_handle(redEdges[j]) << " " << redConn[j] << " "
<< redConn[j1] << " area : " << area << "\n";
- if (fabs(area) < epsilon_1)
+ if (fabs(area) < epsilon_1/2)
{
// found the edge; now find if there is a point in the list here
//std::vector<EntityHandle> * expts = extraNodesMap[redEdges[j]];
@@ -261,8 +291,6 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
std::cout << " found node:" << foundIds[i] << std::endl;
}
}
- delete[] coords1;
-
if (!found)
{
// create a new point in 2d (at the intersection)
@@ -277,29 +305,37 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
if (dbg_1)
std::cout << " new node: " << outNode << std::endl;
}
-
+ delete[] coords1;
}
}
}
if (!found)
{
- std::cout << " red quad: ";
- for (int j1 = 0; j1 < 4; j1++)
+ std::cout << " red polygon: ";
+ for (int j1 = 0; j1 < nsRed; j1++)
{
- std::cout << redQuad[2 * j1] << " " << redQuad[2 * j1 + 1] << "\n";
+ std::cout << redCoords2D[2 * j1] << " " << redCoords2D[2 * j1 + 1] << "\n";
}
- std::cout << " a point pp is not on a red quad " << *pp << " " << pp[1]
- << " red quad " << mb->id_from_handle(red) << " \n";
+ std::cout << " a point pp is not on a red polygon " << *pp << " " << pp[1]
+ << " red polygon " << mb->id_from_handle(red) << " \n";
return 1;
}
}
- // nodes might still be duplicated
+ if (dbg_1)
+ {
+ std::cout << " candidate polygon: nP " << nP << "\n";
+ for (int i1 = 0; i1 < nP; i1++)
+ std::cout << iP[2 * i1] << " " << iP[2 * i1 + 1] << " " << foundIds[i1] << "\n";
+ }
+ // first, find out if we have nodes collapsed; shrink them
+ // we may have to reduce nP
+ // it is possible that some nodes are collapsed after intersection only
+ // nodes will always be in order (convex intersection)
correct_polygon(foundIds, nP);
// now we can build the triangles, from P array, with foundIds
// we will put them in the out set
if (nP >= 3)
{
-
EntityHandle polyNew;
mb->create_element(MBPOLYGON, foundIds, nP, polyNew);
mb->add_entities(outSet, &polyNew, 1);
@@ -321,7 +357,7 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
std::cout << " polygon " << mb->id_from_handle(polyNew) << " nodes: " << nP << " :";
for (int i1 = 0; i1 < nP; i1++)
std::cout << " " << mb->id_from_handle(foundIds[i1]);
- //std::cout << " plane: " << plane << "\n";
+ std::cout << "\n";
std::vector<CartVect> posi(nP);
mb->get_coords(foundIds, nP, &(posi[0][0]));
for (int i1 = 0; i1 < nP; i1++)
@@ -332,20 +368,18 @@ int Intx2MeshInPlane::findNodes(EntityHandle red, EntityHandle blue, double * iP
mb->write_mesh(fff.str().c_str(), &outSet, 1);
}
-
-
-
}
delete[] foundIds;
foundIds = NULL;
return 0;
+ // end copy
}
bool Intx2MeshInPlane::is_inside_element(double xyz[3], EntityHandle eh)
{
int num_nodes;
ErrorCode rval = mb->get_connectivity(eh, redConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != 4 || num_nodes!=3)
+ if (MB_SUCCESS != rval)
return false;
int nsides = num_nodes;
@@ -358,14 +392,14 @@ bool Intx2MeshInPlane::is_inside_element(double xyz[3], EntityHandle eh)
{
// populate coords in the plane for decision making
// they should be oriented correctly, positively
- redQuad[2 * j] = redCoords[j][0];
- redQuad[2 * j + 1] = redCoords[j][1];
+ redCoords2D[2 * j] = redCoords[j][0];
+ redCoords2D[2 * j + 1] = redCoords[j][1];
}
double pt[2]={xyz[0], xyz[1]};// xy plane only
// now, is the projected point inside the red quad?
// cslam utils
- if (point_in_interior_of_convex_polygon (redQuad, num_nodes, pt))
+ if (point_in_interior_of_convex_polygon (redCoords2D, num_nodes, pt))
return true;
return false;
}
diff --git a/tools/mbcslam/Intx2MeshInPlane.hpp b/tools/mbcslam/Intx2MeshInPlane.hpp
index 1f58ddf..fcba728 100644
--- a/tools/mbcslam/Intx2MeshInPlane.hpp
+++ b/tools/mbcslam/Intx2MeshInPlane.hpp
@@ -17,9 +17,10 @@ public:
virtual ~Intx2MeshInPlane();
int computeIntersectionBetweenRedAndBlue(EntityHandle red, EntityHandle blue,
- double * P, int & nP, double & area, int markb[4], int markr[4], bool check_boxes_first=false );
+ double * P, int & nP, double & area, int markb[MAXEDGES], int markr[MAXEDGES],
+ int & nsBlue, int & nsRed, bool check_boxes_first=false);
- int findNodes(EntityHandle red, EntityHandle blue, double * iP, int nP);
+ int findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue, double * iP, int nP);
bool is_inside_element(double xyz[3], EntityHandle eh);
};
diff --git a/tools/mbcslam/Intx2MeshOnSphere.cpp b/tools/mbcslam/Intx2MeshOnSphere.cpp
index 2ea349f..c8590ad 100644
--- a/tools/mbcslam/Intx2MeshOnSphere.cpp
+++ b/tools/mbcslam/Intx2MeshOnSphere.cpp
@@ -22,16 +22,14 @@ Intx2MeshOnSphere::~Intx2MeshOnSphere()
// TODO Auto-generated destructor stub
}
-/* the red quad is convex for sure, intersect it with the blue quad
- * if the blue is convex, intersection is convex, it is pretty easy to order them and eliminate doubles
- * first decide the projection plane, then do a gnomonic projection of both,
+/* the elements are convex for sure, then do a gnomonic projection of both,
* compute intersection in the plane, then go back to the sphere for the points
- * For the time being, blue is convex too, but later on, we should separate the case of concave blue
- */
+ * */
int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, EntityHandle blue,
- double * P, int & nP, double & area, int markb[4], int markr[4], bool check_boxes_first)
+ double * P, int & nP, double & area, int markb[MAXEDGES], int markr[MAXEDGES],
+ int & nsBlue, int & nsRed, bool check_boxes_first)
{
- // the points will be at most 16; they will describe a convex patch, after the points will be ordered and
+ // the points will be at most 40; they will describe a convex patch, after the points will be ordered and
// collapsed (eliminate doubles)
// the area is not really required, except to see if it is greater than 0
@@ -42,35 +40,38 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
int num_nodes;
ErrorCode rval = mb->get_connectivity(red, redConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != 4)
+ if (MB_SUCCESS != rval )
return 1;
- int nsides = num_nodes;
+ nsRed = num_nodes;
//CartVect coords[4];
rval = mb->get_coords(redConn, num_nodes, &(redCoords[0][0]));
if (MB_SUCCESS != rval)
return 1;
- CartVect middle = 0.25
- * (redCoords[0] + redCoords[1] + redCoords[2] + redCoords[3]);
+ CartVect middle = redCoords[0];
+ for (int i=1; i<nsRed; i++)
+ middle += redCoords[i];
+ middle = 1./nsRed * middle;
decide_gnomonic_plane(middle, plane);// output the plane
//CartVect bluecoords[4];
rval = mb->get_connectivity(blue, blueConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != 4)
+ if (MB_SUCCESS != rval )
return 1;
- rval = mb->get_coords(blueConn, num_nodes, &(blueCoords[0][0]));
+ nsBlue = num_nodes;
+ rval = mb->get_coords(blueConn, nsBlue, &(blueCoords[0][0]));
if (MB_SUCCESS != rval)
return 1;
if (dbg_1)
{
std::cout << "red " << mb->id_from_handle(red) << "\n";
- for (int j = 0; j < num_nodes; j++)
+ for (int j = 0; j < nsRed; j++)
{
std::cout << redCoords[j] << "\n";
}
std::cout << "blue " << mb->id_from_handle(blue) << "\n";
- for (int j = 0; j < num_nodes; j++)
+ for (int j = 0; j < nsBlue; j++)
{
std::cout << blueCoords[j] << "\n";
}
@@ -83,19 +84,22 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
if (check_boxes_first)
{
// look at the boxes formed with vertices; if they are far away, return false early
- if (!GeomUtil::bounding_boxes_overlap(redCoords, num_nodes, blueCoords, num_nodes, box_error))
+ if (!GeomUtil::bounding_boxes_overlap(redCoords, nsRed, blueCoords, nsBlue, box_error))
return 0; // no error, but no intersection, decide early to get out
}
- for (int j = 0; j < nsides; j++)
+ for (int j = 0; j < nsRed; j++)
{
// populate coords in the plane for intersection
// they should be oriented correctly, positively
- int rc = gnomonic_projection(redCoords[j], R, plane, redQuad[2 * j],
- redQuad[2 * j + 1]);
+ int rc = gnomonic_projection(redCoords[j], R, plane, redCoords2D[2 * j],
+ redCoords2D[2 * j + 1]);
if (rc != 0)
return 1;
- rc = gnomonic_projection(blueCoords[j], R, plane, blueQuad[2 * j],
- blueQuad[2 * j + 1]);
+ }
+ for (int j=0; j<nsBlue; j++)
+ {
+ int rc = gnomonic_projection(blueCoords[j], R, plane, blueCoords2D[2 * j],
+ blueCoords2D[2 * j + 1]);
if (rc != 0)
return 1;
}
@@ -103,27 +107,33 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
{
std::cout << "gnomonic plane: " << plane << "\n";
std::cout << " red blue\n";
- for (int j = 0; j < 4; j++)
+ for (int j = 0; j < nsRed; j++)
{
- std::cout << redQuad[2 * j] << " " << redQuad[2 * j + 1] << " "
- << blueQuad[2 * j] << " " << blueQuad[2 * j + 1] << "\n";
+ std::cout << redCoords2D[2 * j] << " " << redCoords2D[2 * j + 1] << "\n";
+ }
+ for (int j = 0; j < nsBlue; j++)
+ {
+ std::cout << blueCoords2D[2 * j] << " " << blueCoords2D[2 * j + 1] << "\n";
}
}
- int ret = EdgeIntersections2(blueQuad, redQuad, nsides, markb, markr, P, nP);
+ int ret = EdgeIntersections2(blueCoords2D, nsBlue, redCoords2D, nsRed, markb, markr, P, nP);
if (ret != 0)
return 1; // some unforeseen error
- int side[4] = { 0 };
- int extraPoints = borderPointsOfXinY2(blueQuad, redQuad, nsides, &(P[2 * nP]), side);
+ int side[MAXEDGES] = { 0 };// this refers to what side? blue or red?
+ int extraPoints = borderPointsOfXinY2(blueCoords2D, nsBlue, redCoords2D, nsRed, &(P[2 * nP]), side);
if (extraPoints >= 1)
{
- for (int k = 0; k < nsides; k++)
+ for (int k = 0; k < nsBlue; k++)
{
if (side[k])
{
+ // this means that vertex k of blue is inside convex red; mark edges k-1 and k in blue,
+ // as being "intersected" by red; (even though they might not be intersected by other edges,
+ // the fact that their apex is inside, is good enough)
markb[k] = 1;
- markb[(k + nsides-1) % nsides] = 1; // it is the previous edge, actually, but instead of doing -1, it is
+ markb[(k + nsBlue-1) % nsBlue] = 1; // it is the previous edge, actually, but instead of doing -1, it is
// better to do modulo +3 (modulo 4)
// null side b for next call
side[k]=0;
@@ -132,15 +142,16 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
}
nP += extraPoints;
- extraPoints = borderPointsOfXinY2(redQuad, blueQuad, nsides, &(P[2 * nP]), side);
+ extraPoints = borderPointsOfXinY2(redCoords2D, nsRed, blueCoords2D, nsBlue, &(P[2 * nP]), side);
if (extraPoints >= 1)
{
- for (int k = 0; k < 4; k++)
+ for (int k = 0; k < nsRed; k++)
{
if (side[k])
{
+ // this is to mark that red edges k-1 and k are intersecting blue
markr[k] = 1;
- markr[(k + nsides-1) % nsides] = 1; // it is the previous edge, actually, but instead of doing -1, it is
+ markr[(k + nsRed-1) % nsRed] = 1; // it is the previous edge, actually, but instead of doing -1, it is
// better to do modulo +3 (modulo 4)
// null side b for next call
}
@@ -169,13 +180,11 @@ int Intx2MeshOnSphere::computeIntersectionBetweenRedAndBlue(EntityHandle red, En
// also, we could just create new vertices every time, and merge only in the end;
// could be too expensive, and the tolerance for merging could be an
// interesting topic
-int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * iP, int nP)
+int Intx2MeshOnSphere::findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue,
+ double * iP, int nP)
{
// first of all, check against red and blue vertices
//
- int nsides = 4;
- if (type == MBTRI)
- nsides=3;
if (dbg_1)
{
std::cout << "red, blue, nP, P " << mb->id_from_handle(red) << " "
@@ -187,11 +196,11 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * i
// get the edges for the red triangle; the extra points will be on those edges, saved as
// lists (unordered)
- EntityHandle redEdges[4];// at most 4
+ std::vector<EntityHandle> redEdges(nsRed);//
int i = 0;
- for (i = 0; i < nsides; i++)
+ for (i = 0; i < nsRed; i++)
{
- EntityHandle v[2] = { redConn[i], redConn[(i + 1) % nsides] };
+ EntityHandle v[2] = { redConn[i], redConn[(i + 1) % nsRed] };
std::vector<EntityHandle> adj_entities;
ErrorCode rval = mb->get_adjacencies(v, 2, 1, false, adj_entities,
Interface::INTERSECT);
@@ -214,10 +223,10 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * i
// priority is the red mesh (mb2?)
int j = 0;
EntityHandle outNode = (EntityHandle) 0;
- for (j = 0; j < nsides && !found; j++)
+ for (j = 0; j < nsRed && !found; j++)
{
//int node = redTri.v[j];
- double d2 = dist2(pp, &redQuad[2 * j]);
+ double d2 = dist2(pp, &redCoords2D[2 * j]);
if (d2 < epsilon_1)
{
@@ -225,15 +234,15 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * i
found = 1;
if (dbg_1)
std::cout << " red node j:" << j << " id:"
- << mb->id_from_handle(redConn[j]) << " 2d coords:" << redCoords[2 * j] << " "
- << redCoords[2 * j + 1] << " d2: " << d2 << " \n";
+ << mb->id_from_handle(redConn[j]) << " 2d coords:" << redCoords2D[2 * j] << " "
+ << redCoords2D[2 * j + 1] << " d2: " << d2 << " \n";
}
}
- for (j = 0; j < nsides && !found; j++)
+ for (j = 0; j < nsBlue && !found; j++)
{
//int node = blueTri.v[j];
- double d2 = dist2(pp, &blueQuad[2 * j]);
+ double d2 = dist2(pp, &blueCoords2D[2 * j]);
if (d2 < epsilon_1)
{
// suspect is blueConn[j] corresponding in mbOut
@@ -248,12 +257,12 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * i
}
if (!found)
{
- // find the edge it belongs, first, on the red quad
+ // find the edge it belongs, first, on the red element
//
- for (j = 0; j < nsides; j++)
+ for (j = 0; j < nsRed; j++)
{
- int j1 = (j + 1) % nsides;
- double area = area2D(&redQuad[2 * j], &redQuad[2 * j1], pp);
+ int j1 = (j + 1) % nsRed;
+ double area = area2D(&redCoords2D[2 * j], &redCoords2D[2 * j1], pp);
if (dbg_1)
std::cout << " edge " << j << ": "
<< mb->id_from_handle(redEdges[j]) << " " << redConn[j] << " "
@@ -305,9 +314,9 @@ int Intx2MeshOnSphere::findNodes(EntityHandle red, EntityHandle blue, double * i
if (!found)
{
std::cout << " red quad: ";
- for (int j1 = 0; j1 < nsides; j1++)
+ for (int j1 = 0; j1 < nsRed; j1++)
{
- std::cout << redQuad[2 * j1] << " " << redQuad[2 * j1 + 1] << "\n";
+ std::cout << redCoords2D[2 * j1] << " " << redCoords2D[2 * j1 + 1] << "\n";
}
std::cout << " a point pp is not on a red quad " << *pp << " " << pp[1]
<< " red quad " << mb->id_from_handle(red) << " \n";
@@ -371,9 +380,9 @@ bool Intx2MeshOnSphere::is_inside_element(double xyz[3], EntityHandle eh)
int num_nodes;
ErrorCode rval = mb->get_connectivity(eh, redConn, num_nodes);
- if (MB_SUCCESS != rval || num_nodes != 4)
+ if (MB_SUCCESS != rval)
return false;
- int nsides = num_nodes;
+ int nsRed = num_nodes;
//CartVect coords[4];
rval = mb->get_coords(redConn, num_nodes, &(redCoords[0][0]));
@@ -384,12 +393,12 @@ bool Intx2MeshOnSphere::is_inside_element(double xyz[3], EntityHandle eh)
center += redCoords[k];
center = 1./num_nodes*center;
decide_gnomonic_plane(center, plane);// output the plane
- for (int j = 0; j < nsides; j++)
+ for (int j = 0; j < nsRed; j++)
{
// populate coords in the plane for decision making
// they should be oriented correctly, positively
- int rc = gnomonic_projection(redCoords[j], R, plane, redQuad[2 * j],
- redQuad[2 * j + 1]);
+ int rc = gnomonic_projection(redCoords[j], R, plane, redCoords2D[2 * j],
+ redCoords2D[2 * j + 1]);
if (rc != 0)
return false;
}
@@ -401,7 +410,7 @@ bool Intx2MeshOnSphere::is_inside_element(double xyz[3], EntityHandle eh)
// now, is the projected point inside the red quad?
// cslam utils
- if (point_in_interior_of_convex_polygon (redQuad, 4, pt))
+ if (point_in_interior_of_convex_polygon (redCoords2D, nsRed, pt))
return true;
return false;
}
diff --git a/tools/mbcslam/Intx2MeshOnSphere.hpp b/tools/mbcslam/Intx2MeshOnSphere.hpp
index fdc5b20..77cd830 100644
--- a/tools/mbcslam/Intx2MeshOnSphere.hpp
+++ b/tools/mbcslam/Intx2MeshOnSphere.hpp
@@ -23,9 +23,11 @@ public:
int computeIntersectionBetweenRedAndBlue(EntityHandle red, EntityHandle blue,
- double * P, int & nP, double & area, int markb[4], int markr[4], bool check_boxes_first=false);
+ double * P, int & nP, double & area, int markb[MAXEDGES], int markr[MAXEDGES],
+ int & nsBlue, int & nsRed, bool check_boxes_first=false);
- int findNodes(EntityHandle red, EntityHandle blue, double * iP, int nP);
+ int findNodes(EntityHandle red, int nsRed, EntityHandle blue, int nsBlue,
+ double * iP, int nP);
bool is_inside_element(double xyz[3], EntityHandle eh);
diff --git a/tools/mbcslam/case1_test.cpp b/tools/mbcslam/case1_test.cpp
index c15e3c7..4814e0f 100644
--- a/tools/mbcslam/case1_test.cpp
+++ b/tools/mbcslam/case1_test.cpp
@@ -32,7 +32,7 @@ using namespace moab;
double gtol = 0.0001; // this is for geometry tolerance
std::string input_mesh_file("eulerHomme.vtk"); // input file, plain vtk, sphere-cube mesh
double CubeSide = 6.; // the above file starts with cube side 6; radius depends on cube side
-double t = 0.1, delta_t = 0.01; // check the script
+double t = 0.1, delta_t = 0.43; // check the script
ErrorCode manufacture_lagrange_mesh_on_sphere(Interface * mb,
EntityHandle euler_set, EntityHandle & lagr_set)
@@ -184,7 +184,7 @@ int main(int argc, char **argv)
worker.SetRadius(radius);
- worker.SetEntityType(MBQUAD);
+ //worker.SetEntityType(MBQUAD);
worker.SetErrorTolerance(gtol);
std::cout << "error tolerance epsilon_1=" << gtol << "\n";
diff --git a/tools/mbcslam/cslam_par_test.cpp b/tools/mbcslam/cslam_par_test.cpp
index 9887047..788ff21 100644
--- a/tools/mbcslam/cslam_par_test.cpp
+++ b/tools/mbcslam/cslam_par_test.cpp
@@ -2,7 +2,7 @@
* cslam_par_test.cpp
* test to trigger intersection on a sphere in parallel
* it will start from an eulerian mesh (mesh + velocity from Mark Taylor)
- * file: VELO00.h5m; mesh + velo at 850 milibar, in an unstrucured mesh refined from
+ * file: VELO00.h5m; mesh + velo at 850 milibar, in an unstructured mesh refined from
* a cube sphere grid
* the mesh is read in parallel (euler mesh);
*
@@ -190,7 +190,7 @@ void test_intx_in_parallel_elem_based()
worker.SetRadius(Radius);
worker.set_box_error(EPS1);//
- worker.SetEntityType(MBQUAD);
+ //worker.SetEntityType(MBQUAD);
worker.SetErrorTolerance(Radius*1.e-8);
std::cout << "error tolerance epsilon_1="<< Radius*1.e-8 << "\n";
This diff is so big that we needed to truncate the remainder.
https://bitbucket.org/fathomteam/moab/commits/b4039938515d/
Changeset: b4039938515d
Branch: master
User: iulian07
Date: 2013-06-19 05:55:11
Summary: Merge branch 'master' of bitbucket.org:fathomteam/moab
Affected #: 2 files
diff --git a/src/io/NCHelperEuler.cpp b/src/io/NCHelperEuler.cpp
index a01d4b8..85b6887 100644
--- a/src/io/NCHelperEuler.cpp
+++ b/src/io/NCHelperEuler.cpp
@@ -147,7 +147,7 @@ ErrorCode NCHelperEuler::init_mesh_vals(const FileOptions& opts, EntityHandle fi
#ifdef USE_MPI
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
- for (int i = 0; i < 3; i++)
+ for (int i = 0; i < 2; i++)
parData.gPeriodic[i] = globallyPeriodic[i];
parData.partMethod = partMethod;
int pdims[3];
diff --git a/src/io/NCHelperFV.cpp b/src/io/NCHelperFV.cpp
index 566b719..93b1a99 100644
--- a/src/io/NCHelperFV.cpp
+++ b/src/io/NCHelperFV.cpp
@@ -158,7 +158,7 @@ ErrorCode NCHelperFV::init_mesh_vals(const FileOptions& opts, EntityHandle file_
#ifdef USE_MPI
for (int i = 0; i < 6; i++)
parData.gDims[i] = gDims[i];
- for (int i = 0; i < 3; i++)
+ for (int i = 0; i < 2; i++)
parData.gPeriodic[i] = globallyPeriodic[i];
parData.partMethod = partMethod;
int pdims[3];
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