MeshKit  1.0
ProjectShell.cpp
Go to the documentation of this file.
00001 #include "meshkit/ProjectShell.hpp"
00002 #include "meshkit/MKCore.hpp"
00003 #include "meshkit/ModelEnt.hpp"
00004 #include "meshkit/SizingFunction.hpp"
00005 #include "meshkit/RegisterMeshOp.hpp"
00006 #include "moab/ReadUtilIface.hpp"
00007 
00008 
00009 #include <math.h>
00010 #include <sstream>
00011 #include <queue>
00012 #include <algorithm>
00013 #include <string.h>
00014 
00015 #define ERROR(msg)                                        \
00016   do {                                                            \
00017     Error err(0, "%s, line %d: %s", __FILE__, __LINE__, (msg));\
00018     throw err;                                                    \
00019   } while(false)
00020 
00021 #define ERROR1(fmt, arg1)                                          \
00022   do {                                                             \
00023     Error err(0, "%s, line %d: " fmt, __FILE__, __LINE__, arg1);   \
00024     throw err;                                                     \
00025   } while(false)
00026 
00027 #define ERROR2(fmt, arg1,arg2)                                        \
00028   do {                                                                \
00029     Error err(0, "%s, line %d: " fmt, __FILE__, __LINE__, arg1,arg2); \
00030     throw err;                                                        \
00031   } while(false)
00032 
00033 #define ERROR3(fmt, arg1,arg2,arg3)                                        \
00034   do {                                                                     \
00035     Error err(0, "%s, line %d: " fmt, __FILE__, __LINE__, arg1,arg2,arg3); \
00036     throw err;                                                             \
00037   } while(false)
00038 
00039 
00040 namespace MeshKit {
00041   
00042 double epsilon = 1.e-5; // cm, for coincident points in P, the intersection area
00043 
00044 //---------------------------------------------------------------------------//
00045 // Construction Function for Project Shell
00046 ProjectShell::ProjectShell(MKCore *mk_core, const MEntVector &me_vec) 
00047   : MeshScheme(mk_core, me_vec), 
00048     m_numNodes(0),
00049     m_numTriangles(0),
00050     m_xyz(NULL),  // original coordinates
00051     m_triangles(NULL), // original triangles
00052     ps_edges(NULL),
00053     ps_nodes(NULL),
00054     ps_triangles(NULL),
00055     m_xy(NULL), // 2d coordinates after projection
00056     m_redMesh(NULL),
00057     m_blueMesh(NULL),
00058     m_num2dPoints(0),
00059     m_dbg(false)
00060 {
00061   double direction[3] = {1.0, 0.0, 0.0};
00062   m_direction = Vector<3>(direction);
00063 }
00064   
00065 //---------------------------------------------------------------------------//
00066 // return the type of entities this mesh creates
00067 void ProjectShell::mesh_types(std::vector<moab::EntityType> &tps)
00068 {
00069   tps.push_back(moab::MBVERTEX);
00070   tps.push_back(moab::MBEDGE);
00071   tps.push_back(moab::MBTRI);
00072 }
00073   
00074 
00075 //set up the projection parameters: the direction of projection
00076 void ProjectShell::setup_this() 
00077 {
00078   // Need to set the direction here
00079 }
00080 
00081 //---------------------------------------------------------------------------//
00082 // execute function: Generate the projected 2D
00083 void ProjectShell::execute_this(){
00084   double dist1= length(m_direction);
00085   if (dist1==0.0) ERROR("Null direction of projection");
00086   // normalize direction
00087   m_direction/=dist1;
00088   int ret = getMesh();
00089   if(ret) ERROR("Bad input mesh");
00090   
00091   // we have now the 3D mesh
00092   // verify orientation of the triangles; is it manifold?
00093   ret =  checkMeshValidity();
00094   if(ret) ERROR("Bad shell");
00095   
00096   ret =  projectIn2D();
00097   if(ret) ERROR("Cannot project in 2D");
00098   
00099   ret =  computeIntersections();
00100   if(ret) ERROR("Error in computing intersections");
00101 
00102   ret = commitMesh();
00103   if(ret) ERROR("Error in committing the computed shell");
00104 }
00105 
00106 
00107 int ProjectShell::getEdge(ProjectShell::Node & v1, ProjectShell::Node & v2, int & edgeId, ProjectShell::Edge * edgesArr, int sizeArr)
00108 {
00109   // not found -- return 0 
00110   std::list<int> & L = v1.edges;
00111   std::list<int>::iterator i;
00112   for (i=L.begin(); i!=L.end(); i++)
00113     {
00114       int edge = *i;
00115       int ae = abs(edge);
00116       if (ae>sizeArr || ae ==0 ) ERROR("Bad index in getEdge");
00117       
00118       ProjectShell::Edge & pse = edgesArr[ae-1];
00119       //      int sig = (edge<0) ? -1 : 1;
00120       if (pse.v[0] == v1.id && pse.v[1] == v2.id) 
00121         {
00122           edgeId = edge;
00123           return 1;
00124         }
00125       if (pse.v[0] == v2.id && pse.v[1] == v1.id)
00126         {
00127           edgeId = edge;
00128           return 1;
00129         }
00130     }
00131   return 0;// did not find any edge between v1 and v2 
00132 }
00133   
00134 // Destructor
00135 ProjectShell::~ProjectShell()
00136 {
00137   delete [] ps_edges;
00138   delete [] ps_nodes;
00139   delete [] ps_triangles;
00140   delete [] m_xyz;
00141   delete [] m_triangles;
00142   delete [] m_xy;
00143   delete [] m_redMesh;
00144   delete [] m_blueMesh;
00145 }
00146 #if 0
00147 int ProjectShell::getShellMesh()
00148 {
00149     
00150   // get original coordinates of mesh
00151   iBase_EntityHandle *verts = NULL;
00152   int verts_alloc = 0;
00153   
00154   int vert_coords_alloc = 0;
00155   int vertex_coord_size;
00156   
00157   /* check storage order */
00158   int result;
00159   int this_order;
00160   iMesh_Instance m_mesh = NULL;
00161   iBase_EntitySetHandle m_hRootSet = NULL;
00162   
00163   iMesh_getDfltStorage(m_mesh, &this_order, &result);
00164   if (iBase_SUCCESS != result) 
00165     ERROR("failed to get preferred storage order in getMesh");
00166   
00167   
00168   /* now get the vertex coordinates from a vertex array */
00169   /* need to get the vertices in the model */
00170   verts = NULL;
00171   verts_alloc = 0;
00172   iMesh_getEntities(m_mesh, m_hRootSet, iBase_VERTEX, 
00173                     iMesh_POINT, &verts, &verts_alloc, &m_numNodes, &result);
00174   IBERRCHK(result, "Failed to get vertices");
00175   /*
00176     iMesh_getNumOfType          (       iMesh_Instance           instance,
00177     const iBase_EntitySetHandle         entity_set_handle,
00178     const int   entity_type,6
00179     int *       num_type,
00180     int *       err      
00181     ) */
00182   int numEdges=0;
00183   iMesh_getNumOfType    (m_mesh, m_hRootSet, iBase_EDGE, 
00184                          &numEdges, &result);
00185   
00186   // get the triangles and the vertices in one shot
00187   iBase_EntityHandle *triangles = NULL;
00188   int triangles_alloc = 0;
00189   iBase_EntityHandle *vert_adj = NULL;
00190   int  vert_adj_alloc = 0, vert_adj_size;
00191   int * offsets = NULL, offsets_alloc = 0, indices_size;
00192   int * indices = NULL, indices_alloc = 0, offsets_size;
00193   iMesh_getAdjEntIndices( m_mesh, m_hRootSet, 
00194                           iBase_FACE, iMesh_TRIANGLE, iBase_VERTEX,
00195                           &triangles, &triangles_alloc, &m_numTriangles,
00196                           &vert_adj, &vert_adj_alloc, &vert_adj_size,
00197                           &indices, &indices_alloc, &indices_size,
00198                           &offsets, &offsets_alloc, &offsets_size,
00199                           &result );
00200   if (iBase_SUCCESS != result) {
00201     ERROR("failed to get triangles and vertices.");
00202   }
00203   
00204   
00205   /* get the coordinates in one array */
00206   
00207   vert_coords_alloc = 0;
00208   
00209   iMesh_getVtxArrCoords(m_mesh, vert_adj, vert_adj_size, iBase_INTERLEAVED, 
00210                         &m_xyz, &vert_coords_alloc, &vertex_coord_size, &result);
00211   if (iBase_SUCCESS != result || 3*m_numNodes!=vertex_coord_size) {
00212     ERROR("failed to get vertex coordinates of entities in getMeshData.");
00213   }
00214   
00215   
00216   // build list of edges; we cannot rely on iMesh to give them to us
00217   // can we force imesh to give the edges? It does not seem like that
00218   
00219   // the vertices are identified as the index in vert_adj 
00220   // indices are the triangles
00221   // i is index in triangles
00222   //  offsets[i], offsets[i+1] are offsets in indices
00223   // vertices of triangle [i] have indices  indices[offsets[i]+j], j=0:(offsets[i+1]-offsets[i])-1
00224   
00225   // first, determine the edges of triangle i
00226   if (offsets_size - 1 != m_numTriangles)
00227   {
00228       ERROR("bad indexing");
00229   }
00230   int i=0;// index used everywhere
00231   for (i=0; i<m_numTriangles; i++)
00232   {
00233     if (offsets[i+1]-offsets[i] !=3)
00234     {
00235           ERROR("Not a triangle");
00236     }
00237   }
00238   // build edges from triangle information
00239   ps_edges = new ProjectShell::Edge [m_numTriangles*3];// overestimate; probably only half will be created
00240   ps_nodes = new ProjectShell::Node [m_numNodes];
00241   ps_triangles = new ProjectShell::Triangle3D [m_numTriangles];
00242   //
00243   for (i=0; i<m_numNodes; i++)
00244   {
00245     ProjectShell::Node & psn = ps_nodes[i];
00246     psn.id = i;// this can be found by address in the array, but why bother
00247     for (int j=0; j<3; j++)
00248     {
00249       psn.xyz[j] = m_xyz[3*i+j];
00250     }
00251   }
00252   // index in edge: 
00253   int edgeIndex = 0;
00254   int j=0;
00255   for (j=0; j<m_numTriangles; j++)
00256   {
00257     ProjectShell::Triangle3D & curTria = ps_triangles[j];
00258       int ii=0;
00259       for (ii=0; ii<3; ii++)
00260         curTria.v[ii]=indices[offsets[j]+ii];
00261       for (ii=0; ii<3; ii++)
00262         {
00263           ProjectShell::Node & v1= ps_nodes[curTria.v[ii]];
00264           int nextii = ii+1;
00265           if (ii==2)
00266             nextii=0;
00267           ProjectShell::Node & v2= ps_nodes[curTria.v[nextii]];
00268           // is there an edge from v1 to v2, or from v2 to v1
00269           int edgeId;
00270           int exi = getEdge(v1, v2, edgeId, ps_edges, edgeIndex);// will be created if not existing already
00271           if (exi)
00272           {
00273               curTria.e[ii] = edgeId;
00274               ProjectShell::Edge & foundEdge = ps_edges[abs(edgeId)-1];
00275               foundEdge.used++;
00276               if(foundEdge.used>2)
00277               {
00278                 ERROR("Bad indexing in ps_edge");
00279               }
00280               foundEdge.t[1]=j; // mark the triangle it belongs to
00281           }
00282           else
00283           {
00284             int cId = curTria.e[ii]=edgeIndex+1;
00285             ProjectShell::Edge & newEdge = ps_edges[edgeIndex];
00286             newEdge.v[0] = curTria.v[ii];
00287             newEdge.v[1]=curTria.v[nextii];
00288             v1.edges.push_back(cId); // positive means starts at vertex
00289             v2.edges.push_back(-cId);// negative means incident to the vertex
00290             edgeIndex++;
00291             newEdge.used=1;
00292             newEdge.t[0] = j; // index for triangles
00293           }
00294         }
00295   }
00296   m_numEdges = edgeIndex;
00297   // fill up now the triangle neighbor information
00298   for (j=0; j<m_numTriangles; j++)
00299     {
00300       ProjectShell::Triangle3D & tri = ps_triangles[j];
00301       for (int k=0; k<3; k++)
00302         {
00303           int edgeId = tri.e[k];
00304           int ae = abs(edgeId);
00305           ProjectShell::Edge & edge = ps_edges[ae-1];
00306           int t0 = edge.t[0], t1 = edge.t[1];
00307           if (t0==j)
00308             tri.t[k] = t1;
00309           else if (t1==j)
00310             tri.t[k] = t0;
00311         }
00312     }
00313   
00314   free(verts);
00315   free(triangles);
00316   free(vert_adj);
00317   free (indices);
00318   free (offsets);
00319   
00320   //free(vert_coords);
00321   return 0;
00322 }
00323 #endif
00324 int ProjectShell::checkMeshValidity()
00325 {
00326   // basically, see if the edges are all used 2 times, and each is is used in each direction
00327   int i=0;
00328   for (i=0; i<m_numEdges; i++)
00329     {
00330       ProjectShell::Edge & edge = ps_edges[i];
00331       if (edge.used!=2)
00332         {
00333           ERROR("Edge not used exactly 2 times");
00334         }
00335       int edgeId = i+1;
00336       int sum=0;
00337       for (int k=0; k<2; k++)
00338         {
00339           ProjectShell::Triangle3D & tr1=ps_triangles[edge.t[k]];
00340           for (int kk=0; kk<3; kk++)
00341             {
00342               int ec= tr1.e[kk];
00343               if (abs(ec)==edgeId)
00344                 {
00345                   sum+=ec;
00346                   break;
00347                 }
00348               
00349             }
00350         }
00351       if (sum!=0)
00352           ERROR1("Edge %d not used exactly 2 times", edgeId);      
00353     }
00354   // check that all triangles have 3 neighbors
00355   for (i=0; i<m_numTriangles; i++)
00356     {
00357       ProjectShell::Triangle3D & tr = ps_triangles[i];
00358       for (int k=0; k<3; k++)
00359         {
00360           int t = tr.t[k];
00361           if (t==i || t < 0 || t >= m_numTriangles)
00362             ERROR1("triangle %d does not have 3 neighbors", i);
00363         }
00364     }
00365   
00366   return 0; // everything is OK
00367 }
00368 
00369 
00370 int ProjectShell::projectIn2D()
00371 {
00372   // considering the direction, classify each projected triangle as red or blue
00373   // the red ones are positive (inbound), blue ones are negative
00374   // first decide the other 2 vectors
00375   double v[3] = {0.,-1.,0.};
00376   Vector<3> vv(v);
00377   m_dirX = vector_product(m_direction, vv);
00378   
00379   double d1 = length(m_dirX);
00380   if (d1<1.e-5)// consider another direction
00381     {
00382       vv[1]=0.; vv[2]=-1.;
00383       m_dirX = vector_product(m_direction, vv);
00384       d1 = length(m_dirX);
00385       if (d1<1.e-5)
00386         ERROR("cannot find a suitable direction; abort");
00387     }
00388   m_dirX /= d1;
00389   m_dirY = vector_product(m_direction, m_dirX);
00390   // dirY must be already normalized, but why worry
00391   d1 = length(m_dirY);
00392   if (d1==0.)
00393     ERROR("Get out of here, it is hopeless");
00394   m_dirY /= d1;
00395   
00396   // now do the projection in 2d
00397   //   we have 3 vectors, m_dirX, m_dirY, m_direction
00398   // for every point, A, the projection in xy plane will be just
00399   // xA = A . u; yA = A . v  (dirX and dirY)
00400   
00401   // also, it is important to compute the normal 
00402   // some triangles will be reverted and some may be projecting to a line
00403   
00404   // coordinates of the projected nodes on 2D
00405   // what could be a good estimate for total number of points in 2D? 
00406   //  we start with 2d capacity 3* m_numNodes 
00407   m_2dcapacity = m_numNodes*3;
00408   m_num2dPoints = m_numNodes; // directly project 3d nodes in 2d, keep the ids
00409   m_xy = new double [3*m_2dcapacity];
00410   // this array will be parallel with the original mesh
00411   // new nodes will appear after intersection computation
00412   // triangles are characterized by their orientation: negative, positive or 0
00413   
00414   // the neighbors will be preserved
00415   // some edges will be collapsed, and also some triangles
00416   // in those cases, what will be the neighbors?
00417   // flag them with a high number ()
00418   
00419   //m_finalNodes.resize(3*m_numNodes); // the actual size is n_numNodes first
00420   for (int k=0; k<m_numNodes; k++)
00421     {
00422       // double xx= dot( ps_nodes[k].xyz, m_dirX);
00423       // double yy= dot( ps_nodes[k].xyz, m_dirY);
00424       // _finalNodes[k].x =  dot( ps_nodes[k].xyz, m_dirX);
00425       // _finalNodes[k].y =  dot( ps_nodes[k].xyz, m_dirY);
00426       m_xy[3*k] = inner_product( ps_nodes[k].xyz, m_dirX);
00427       m_xy[3*k+1] = inner_product( ps_nodes[k].xyz, m_dirY);
00428       //  m_finalNodes[k].x =  m_xy[2*k];
00429       //  m_finalNodes[k].y =  m_xy[2*k+1];
00430     }
00431   for (int k=0; k<m_2dcapacity; k++)
00432     m_xy[3*k+2] = 0;// the z ccordinate is always 0 !!
00433   // if any edges are collapsed, the corresponding triangles should be collapsed too
00434   // we will collapse the triangles first, then the edges
00435   // when a triangle is collapsed on an edge, it should break in 2 other triangles
00436   // we should really form another triangle array, in 2D 
00437   
00438   // first, loop over edges and see which are eliminated;
00439   
00440   double *edgeLen2D=new double [m_numEdges];
00441   for (int k=0; k<m_numEdges; k++)
00442     {
00443       
00444       int v0 = ps_edges[k].v[0];
00445       int v1 = ps_edges[k].v[1];
00446       edgeLen2D[k] = length(Vector<2>(&m_xy[3*v0])-Vector<2>(&m_xy[3*v1])); 
00447     }
00448   double *triArea = new double [m_numTriangles];
00449   for (int k=0; k<m_numTriangles; k++)
00450     {
00451       const int * pV =&( ps_triangles[k].v[0]);
00452       triArea[k] = area2D(Vector<2>(&m_xy[ 3*pV[0]]),  
00453                           Vector<2>(&m_xy[ 3*pV[1]]),  
00454                           Vector<2>(&m_xy[ 3*pV[2]]));
00455     }
00456   // if an edge is length 0, we must collapse the nodes, and 2 triangles
00457   // first construct a new 2d mesh
00458   // we will have to classify all triangles, edges
00459   // for each triangle we need to know its original triangle
00460   
00461   // first mark all triangles; then count positive, negative and zero area (some tolerance is needed)
00462   int numNeg = 0, numPos = 0, numZero = 0;
00463   // those that are negative will be reverted in the new array
00464   int * newTriId = new int [m_numTriangles];
00465   for (int k=0; k<m_numTriangles ; k++)
00466     {
00467       if (triArea[k] > 0)
00468         {
00469           //numPos++;
00470           newTriId[k] = numPos++;
00471         }
00472       else if (triArea[k] <0)
00473         {
00474           //numNeg ++;
00475           newTriId[k] = numNeg++;
00476         }
00477       else
00478         {
00479           numZero ++;
00480           newTriId[k] = -1;// receive no Id, as it will not be carried along
00481         }
00482     }
00483   // there are 2 groups: red (positive), blue (negative)
00484   // revert the negative ones, and decide the neighbors
00485   
00486   // red triangles are positive, blue are negative
00487   // the negative ones need to be reverted; revert the nodes, first, then edges
00488   // do we really need the edges? or just the neighboring triangles?
00489   // if a neighbor is the other sign or zero, will become boundary marker (large number, m_numTriangles+1) 
00490   // 
00491   // we need to find first 2 triangles (red and blue) that are intersecting
00492   // they will be the seeds
00493   m_redMesh = new ProjectShell::Triangle2D [numPos] ;
00494   m_blueMesh = new  ProjectShell::Triangle2D [numNeg];
00495   m_numPos = numPos;
00496   m_numNeg = numNeg;
00497   // do another loop, to really build the 2D mesh we start with
00498   // we may have potential starting triangles for the marching along, if the sign of one of the neighbors is 
00499   // different
00500   for (int k=0; k<m_numTriangles ; k++)
00501     {
00502       ProjectShell::Triangle3D & orgTria = ps_triangles[k];
00503       if (triArea[k] > 0)
00504         {
00505           //numPos++;
00506           ProjectShell::Triangle2D & redTria= m_redMesh[newTriId[k]];
00507           redTria.oldId = k ; // index 
00508           redTria.area = triArea[k];
00509           for (int j=0; j<3; j++)
00510             {
00511               
00512               // copy the edge information too
00513               redTria.e[j] = orgTria.e[j]; 
00514               // qualify the neighbors
00515               //
00516               int t = orgTria.t[j];
00517               redTria.v[j] = orgTria.v[j];
00518               if (triArea[t]>0)
00519                 {
00520                   redTria.t[j] = newTriId[t];// will be the index in red mesh 
00521                 }
00522               else
00523                 {
00524                   redTria.t[j] = numPos; // marker for boundary
00525                 }
00526             }
00527           
00528         }
00529       else if (triArea[k] <0)
00530         {
00531           //numNeg ++;
00532           ProjectShell::Triangle2D & blueTria= m_blueMesh[newTriId[k]];
00533           blueTria.oldId = k;
00534           blueTria.area =triArea[k]; // this is for debugging I think
00535           for (int j=0; j<3; j++)
00536             {
00537               // copy the edge information too
00538               blueTria.e[j] = orgTria.e[j];
00539               // qualify the neighbors
00540               // 
00541               int t = orgTria.t[j];
00542               blueTria.v[j] = orgTria.v[j];
00543               if (triArea[t]<0)
00544                 {
00545                   blueTria.t[j] = newTriId[t];// will be the index in red mesh 
00546                 }
00547               else
00548                 {
00549                   blueTria.t[j] = numNeg; // marker for boundary
00550                 }
00551             }
00552           
00553         }
00554       else
00555         {
00556           // numZero ++;
00557           // newTriId[k] = -1;// receive no Id, as it will not be carried along
00558           // nothing to do for null triangles
00559         }
00560     }
00561   // revert the blue triangles, so they become positive oriented too
00562   for (int k=0; k<numNeg; k++)
00563     {
00564       // 
00565       ProjectShell::Triangle2D & blueTri = m_blueMesh[k];
00566       int tmp = blueTri.v[1];
00567       blueTri.v[1] = blueTri.v[2];
00568       blueTri.v[2] = tmp;
00569       // this is really stupid: 
00570       // we should have switched triangle 1 and 3, not 2 and 3
00571       // hard to catch
00572       tmp = blueTri.t[0];
00573       blueTri.t[0] = blueTri.t[2];
00574       blueTri.t[2] = tmp;
00575       // reverse the edges too
00576       tmp = blueTri.e[0];
00577       blueTri.e[0] = -blueTri.e[2];
00578       blueTri.e[1] = -blueTri.e[1];
00579       blueTri.e[2] = -tmp; 
00580     }
00581   // at this point, we have red triangles and blue triangles, and 2d nodes array
00582   // we will start creating new triangles, step by step, pointing to the nodes in _finalNodes
00583   m_numCurrentNodes = m_numNodes;
00584   if (m_dbg)
00585     {
00586       m_dbg_out.open("dbg.m"); 
00587       m_dbg_out << "P=[ \n";
00588       for (int i=0; i<m_numNodes; i++)
00589         {
00590           m_dbg_out << m_xy[3*i] << " " << m_xy[3*i+1] << "\n";
00591         }
00592       m_dbg_out << "];\n ";
00593       m_dbg_out << "Ta=[ \n";
00594       for (int k=0; k<m_numPos; k++)
00595         {
00596           ProjectShell::Triangle2D & redTri = m_redMesh[k];
00597           for (int j=0; j<3; j++)
00598             m_dbg_out << redTri.v[j]+1 << " " ;
00599           for (int jj=0; jj<3; jj++)
00600             {
00601               m_dbg_out << redTri.t[jj]+1 << " " ;
00602             }
00603           m_dbg_out << "\n";
00604         }
00605       m_dbg_out << "]; \n";
00606       m_dbg_out << "Tb=[ \n";
00607       for (int kk=0; kk<m_numNeg; kk++)
00608         {
00609           ProjectShell::Triangle2D & blueTri = m_blueMesh[kk];
00610           for (int j=0; j<3; j++)
00611             m_dbg_out << blueTri.v[j]+1 << " " ;
00612           for (int jj=0; jj<3; jj++)
00613             {
00614               m_dbg_out << blueTri.t[jj]+1 << " " ;
00615             }
00616           m_dbg_out << "\n";
00617         }
00618       m_dbg_out << "]; \n";
00619       m_dbg_out.close();
00620       
00621     }
00622   delete [] edgeLen2D;
00623   delete [] triArea;
00624   delete [] newTriId;
00625   return 0; 
00626 }
00627 
00628 int ProjectShell::computeIntersections()
00629 {
00630   // will start at 2 triangles, and advance an intersection front (2 queues, one for red triangles, one for blue ..)
00631   // will mark a red triangle that is processed, and add to the queue
00632   // for each red triangle will find all the blue triangles that are intersected
00633   // find first 2 triangles that intersect: these will be the seeds for intersection
00634   
00635   int startRed=0;
00636   int startBlue = 0;
00637   for (startBlue = 0; startBlue<m_numNeg; startBlue++)
00638     {
00639       double area = 0;
00640       // if area is > 0 , we have intersections
00641       double P[24]; // max 6 points, but it may grow bigger; why worry
00642       int nP = 0;
00643       int n[3];// sides
00644       computeIntersectionBetweenRedAndBlue(/* red */0, startBlue, P, nP, area,n);
00645       if (area>0)
00646         break; // found 2 triangles that intersect; these will be the seeds
00647     }
00648   if (startBlue==m_numNeg)
00649     ERROR("Can't find any triangle");
00650   
00651   // on the red edges, we will keep a list of new points (in 2D)
00652   // they will be used to create or not new points for points P from intersection
00653   // (sometimes the points P are either on sides, or on vertices of blue or even red triangles)
00654   
00655   /*
00656     matlab code: 
00657     function M=InterfaceMatrix(Na,Ta,Nb,Tb);
00658     % INTERFACEMATRIX projection matrix for nonmatching triangular grids 
00659     %   M=InterfaceMatrix(Na,Ta,Nb,Tb); takes two triangular meshes Ta
00660     %   and Tb with associated nodal coordinates in Na and Nb and
00661     %   computes the interface projection matrix M
00662     
00663     bl=[1];                        % bl: list of triangles of Tb to treat
00664     bil=[1];                       % bil: list of triangles Ta to start with
00665     bd=zeros(size(Tb,1)+1,1);      % bd: flag for triangles in Tb treated 
00666     bd(end)=1;                     % guard, to treat boundaries
00667     bd(1)=1;                       % mark first triangle in b list.
00668     M=sparse(size(Nb,2),size(Na,2));
00669     while length(bl)>0
00670     bc=bl(1); bl=bl(2:end);      % bc: current triangle of Tb 
00671     al=bil(1); bil=bil(2:end);   % triangle of Ta to start with
00672     ad=zeros(size(Ta,1)+1,1);    % same as for bd
00673     ad(end)=1;
00674     ad(al)=1; 
00675     n=[0 0 0];                   % triangles intersecting with neighbors
00676     while length(al)>0
00677     ac=al(1); al=al(2:end);    % take next candidate
00678     [P,nc,Mc]=Intersect(Nb(:,Tb(bc,1:3)),Na(:,Ta(ac,1:3)));
00679     if ~isempty(P)             % intersection found
00680     M(Tb(bc,1:3),Ta(ac,1:3))=M(Tb(bc,1:3),Ta(ac,1:3))+Mc;
00681     t=Ta(ac,3+find(ad(Ta(ac,4:6))==0)); 
00682     al=[al t];               % add neighbors 
00683     ad(t)=1;
00684     n(find(nc>0))=ac;        % ac is starting candidate for neighbor  
00685     end
00686     end
00687     tmp=find(bd(Tb(bc,4:6))==0); % find non-treated neighbors
00688     idx=find(n(tmp)>0);          % take those which intersect
00689     t=Tb(bc,3+tmp(idx));
00690     bl=[bl t];                   % and add them
00691     bil=[bil n(tmp(idx))];       % with starting candidates Ta
00692     bd(t)=1;
00693     end
00694   */
00695   std::queue<int> blueQueue; // these are corresponding to Ta,
00696   blueQueue.push( startBlue);
00697   std::queue<int> redQueue;
00698   redQueue.push (startRed);
00699   // the flags are used for marking the triangles already considered
00700   int * blueFlag = new int [m_numNeg+1]; // number of blue triangles + 1, to account for the boundary
00701   int k=0;
00702   for (k=0; k<m_numNeg; k++)
00703     blueFlag[k] = 0;
00704   blueFlag[m_numNeg] = 1; // mark the "boundary"; stop at the boundary
00705   blueFlag[startBlue] = 1; // mark also the first one
00706   // also, red flag is declared outside the loop
00707   int * redFlag = new int [m_numPos+1];
00708   
00709   if (m_dbg)
00710     m_dbg_out.open("patches.m");
00711   while( !blueQueue.empty() )
00712     {
00713       int n[3]; // flags for the side : indices in red mesh start from 0!!! (-1 means not found )
00714       for (k=0; k<3; k++)
00715         n[k] = -1; // a paired red not found yet for the neighbors of blue 
00716       int currentBlue = blueQueue.front();
00717       blueQueue.pop();
00718       for (k=0; k<m_numPos; k++)
00719         redFlag[k] = 0;
00720       redFlag[m_numPos] = 1; // to guard for the boundary
00721       int currentRed = redQueue.front(); // where do we check for redQueue???? 
00722       // red and blue queues are parallel
00723       redQueue.pop();// 
00724       redFlag[currentRed] = 1; // 
00725       std::queue<int> localRed;
00726       localRed.push(currentRed);
00727       while( !localRed.empty())
00728         {
00729           // 
00730           int redT = localRed.front();
00731           localRed.pop();
00732           double P[24], area;
00733           int nP = 0; // intersection points
00734           int nc[3]= {0, 0, 0}; // means no intersection on the side (markers)
00735           computeIntersectionBetweenRedAndBlue(/* red */redT, currentBlue, P, nP, area,nc);
00736           if (nP>0) 
00737             {
00738               // intersection found: output P and original triangles if nP > 2
00739               if (m_dbg && area>0)
00740                 {
00741                   //std::cout << "area: " << area<< " nP:"<<nP << " sources:" << redT+1 << ":" << m_redMesh[redT].oldId <<
00742                   //  " " << currentBlue+1<< ":" << m_blueMesh[currentBlue].oldId << std::endl;
00743                 }
00744               if (m_dbg)
00745                 {
00746                   m_dbg_out << "pa=[\n";
00747                     
00748                   for (k=0; k<nP; k++)
00749                     {
00750                         
00751                       m_dbg_out <<  P[2*k] << "\t " ;
00752                     }
00753                     
00754                   m_dbg_out << "\n";
00755                   for (k=0; k<nP; k++)
00756                     {
00757                         
00758                       m_dbg_out << P[2*k+1] << "\t "; 
00759                     }
00760                     
00761                   m_dbg_out << " ]; \n";
00762                   m_dbg_out << " patch(pa(1,:),pa(2,:),'m');       \n";
00763                 }
00764               // add neighbors to the localRed queue, if they are not marked
00765               for (int nn= 0; nn<3; nn++)
00766                 {
00767                   int neighbor = m_redMesh[redT].t[nn];
00768                   if (redFlag[neighbor] == 0)
00769                     {
00770                       localRed.push(neighbor);
00771                       redFlag[neighbor] =1; // flag it to not be added anymore
00772                     }
00773                   // n(find(nc>0))=ac;        % ac is starting candidate for neighbor
00774                   if (nc[nn]>0) // intersected  
00775                     n[nn] = redT;// start from 0!! 
00776                 }
00777               if (nP>1) // this will also construct triangles, if needed
00778                 findNodes(redT, currentBlue, P, nP);
00779             }
00780             
00781         }
00782       for (int j=0; j<3; j++)
00783         {
00784           int blueNeigh = m_blueMesh[currentBlue].t[j];
00785           if (blueFlag[blueNeigh]==0 && n[j] >=0 ) // not treated yet and marked as a neighbor
00786             {
00787               // we identified triangle n[j] as intersecting with neighbor j of the blue triangle
00788               blueQueue.push(blueNeigh);
00789               redQueue.push(n[j]);
00790               if (m_dbg)
00791                 std::cout << "new triangles pushed: blue, red:" << blueNeigh+1 << " " << n[j]+1 << std::endl;
00792               blueFlag[blueNeigh] = 1;
00793             }
00794         }
00795     }
00796   delete [] redFlag;
00797   redFlag = NULL;
00798     
00799   delete [] blueFlag; // get rid of it
00800   blueFlag = NULL;
00801   if (m_dbg)
00802     m_dbg_out.close();
00803   return 0;
00804 }
00805 
00806  
00807 int borderPointsOfXinY(double * X, double * Y, double * P)
00808 {
00809   // 2 triangles, 3 corners, is the corner of X in Y?
00810   // Y must have a positive area
00811   /*
00812     function P=PointsOfXInY(X,Y);
00813     % POINTSOFXINY finds corners of one triangle within another one
00814     %   P=PointsOfXInY(X,Y); computes for the two given triangles X
00815     %   and Y (point coordinates are stored column-wise, in counter clock
00816     %   order) the corners P of X which lie in the interior of Y.
00817       
00818     k=0;P=[];
00819     v0=Y(:,2)-Y(:,1); v1=Y(:,3)-Y(:,1);  % find interior points of X in Y
00820     d00=v0'*v0; d01=v0'*v1; d11=v1'*v1;  % using baricentric coordinates
00821     id=1/(d00*d11-d01*d01);
00822     for i=1:3 
00823     v2=X(:,i)-Y(:,1); d02=v0'*v2; d12=v1'*v2; 
00824     u=(d11*d02-d01*d12)*id; v=(d00*d12-d01*d02)*id;
00825     if u>=0 & v>=0 & u+v<=1            % also include nodes on the boundary
00826     k=k+1; P(:,k)=X(:,i);
00827     end;
00828     end;
00829   */
00830   int extraPoint = 0;
00831   for (int i=0; i<3; i++)
00832     {
00833       // compute twice the area of all 3 triangles formed by a side of Y and a corner of X; if one is negative, stop
00834       double A[2];
00835       for (int k=0; k<2; k++)
00836         A[k] = X[2*i+k];
00837       int inside=1;
00838       for (int j=0; j<3; j++)
00839         {
00840           double B[2], C[2];
00841           for (int k=0; k<2; k++)
00842             {
00843               B[k] = Y[2*j+k];
00844               int j1 = (j+1)%3;
00845               C[k] = Y[2*j1+k];
00846             }
00847             
00848           double area2 = (B[0]-A[0])*(C[1]-A[1])-(C[0]-A[0])*(B[1]-A[1]);
00849           if (area2<0.)
00850             {
00851               inside = 0; break;
00852             }
00853         }
00854       if (inside)
00855         {
00856           P[extraPoint*2  ] = A[0];
00857           P[extraPoint*2+1] = A[1];
00858           extraPoint++;
00859         }
00860     }
00861   return extraPoint;
00862 }
00863   
00864 /* 
00865    function P=SortAndRemoveDoubles(P);
00866    % SORTANDREMOVEDOUBLES sort points and remove duplicates
00867    %   P=SortAndRemoveDoubles(P); orders polygon corners in P counter
00868    %   clock wise and removes duplicates
00869      
00870    ep=10*eps;                           % tolerance for identical nodes
00871    m=size(P,2); 
00872    if m>0                              
00873    c=sum(P')'/m;                      % order polygon corners counter 
00874    for i=1:m                          % clockwise
00875    d=P(:,i)-c; ao(i)=angle(d(1)+sqrt(-1)*d(2));
00876    end;
00877    [tmp,id]=sort(ao); 
00878    P=P(:,id);
00879    i=1;j=2;                           % remove duplicates
00880    while j<=m
00881    if norm(P(:,i)-P(:,j))>ep
00882    i=i+1;P(:,i)=P(:,j);j=j+1;
00883    else
00884    j=j+1;
00885    end;
00886    end;
00887    P=P(:,1:i);
00888    end;
00889 */
00890   
00891 int swap (double * p , double * q)
00892 {
00893   double tmp = *p;
00894   *p = *q;
00895   *q = tmp;
00896   return 0;
00897 }
00898 int SortAndRemoveDoubles(double * P, int & nP)
00899 {
00900   if (nP<2)
00901     return 0; // nothing to do
00902   // center of gravity for the points
00903   double c[2] = {0., 0.};
00904   int k=0;
00905   for (k=0; k<nP; k++)
00906     {
00907       c[0]+=P[2*k];
00908       c[1]+=P[2*k+1];
00909     }
00910   c[0]/=nP;
00911   c[1]/=nP;
00912   double angle[12]; // could be at most 12 points; much less usually
00913   for (k=0; k<nP; k++)
00914     {
00915       double x = P[2*k] - c[0], y = P[2*k+1] - c[1] ;
00916       if ( x!= 0. || y!=0.)
00917         angle[k] = atan2(y, x);
00918       else
00919         {
00920           angle[k] = 0;
00921           // this means that the points are on a line, or all coincident // degenerate case
00922         }
00923     } 
00924   // sort according to angle; also eliminate close points
00925   int sorted = 1;
00926   do
00927     {
00928       sorted = 1;
00929       for(k=0; k<nP-1; k++)
00930         {
00931           if (angle[k]>angle[k+1])
00932             {
00933               sorted = 0;
00934               swap ( angle+k, angle+k+1);
00935               swap ( P+(2*k), P+(2*k+2));
00936               swap ( P+(2*k+1), P+(2*k+3));
00937             }
00938         }
00939     }
00940   while (!sorted);
00941   // eliminate doubles
00942     
00943   int i=0, j=1; // the next one; j may advance faster than i
00944   // check the unit 
00945   // these are cm; 2 points are the same if the distance is less than epsilon
00946   while (j<nP)
00947     {
00948       double d2 = length(Vector<2>(&P[2*i])-Vector<2>(&P[2*j]));
00949       if (d2 > epsilon)
00950         {
00951           i++;  
00952           P[2*i] = P[2*j];
00953           P[2*i+1] = P[2*j+1]; 
00954         }
00955       j++;
00956     }
00957   // test also the last point with the first one (index 0)
00958     
00959   double d2 = length(Vector<2>(P)-Vector<2>(&P[2*i])); // check the first and last points (ordered from -pi to +pi)
00960   if (d2 > epsilon)
00961     {
00962       nP = i+1;
00963     }
00964   else
00965     nP = i; // effectively delete the last point (that would have been the same with first)
00966   if (nP==0)
00967     nP=1; // we should be left with at least one point we already tested if nP is 0 originally
00968   return 0;
00969 }
00970 
00971 // this method computed intersection between 2 triangles: will output n points, area, affected sides
00972 int ProjectShell::computeIntersectionBetweenRedAndBlue(int red, int blue, double * P, int & nP, double & area,  int mark[3]) {
00973   // the points will be at most 9; they will describe a convex patch, after the points will be ordered and 
00974   // collapsed (eliminate doubles)
00975   // the area is not really required
00976   ProjectShell::Triangle2D & redTri = m_redMesh[red];
00977   ProjectShell::Triangle2D & blueTri = m_blueMesh[blue];
00978   double redTriangle[6];// column wise
00979   double blueTriangle[6];
00980   for (int i=0; i<3; i++)
00981     {
00982       // node i, 2 coordinates
00983       for (int k=0; k<2; k++)
00984         {
00985           int node = redTri.v[i];
00986           redTriangle[2*i+k] = m_xy[3*node+k];
00987             
00988           node = blueTri.v[i];
00989           blueTriangle[2*i+k] =  m_xy[3*node+k];
00990         }
00991     }
00992   /* Matlab source code:
00993      function [P,n,M]=Intersect(X,Y);
00994      % INTERSECT intersection of two triangles and mortar contribution
00995      %   [P,n,M]=Intersect(X,Y); computes for the two given triangles X and
00996      %   Y (point coordinates are stored column-wise, in counter clock
00997      %   order) the points P where they intersect, in n the indices of
00998      %   which neighbors of X are also intersecting with Y, and the local
00999      %   mortar matrix M of contributions of the P1 elements X on the P1
01000      %   element Y. The numerical challenges are handled by including
01001      %   points on the boundary and removing duplicates at the end.
01002        
01003      [P,n]=EdgeIntersections(X,Y);
01004      P1=PointsOfXInY(X,Y);
01005      if size(P1,2)>1                      % if two or more interior points
01006      n=[1 1 1];                         % the triangle is candidate for all 
01007      end                                  % neighbors
01008      P=[P P1];
01009      P=[P PointsOfXInY(Y,X)];
01010      P=SortAndRemoveDoubles(P);           % sort counter clock wise
01011      M=zeros(3,3);
01012      if size(P,2)>0
01013      for j=2:size(P,2)-1                % compute interface matrix
01014      M=M+MortarInt(P(:,[1 j j+1]),X,Y);
01015      end;
01016      patch(P(1,:),P(2,:),'m')           % draw intersection for illustration
01017      %  H=line([P(1,:) P(1,1)],[P(2,:),P(2,1)]);
01018      %  set(H,'LineWidth',3,'Color','m');
01019      pause(0)
01020      end;
01021        
01022        
01023   */
01024   //we do not really need the mortar matrix
01025     
01026   //int n[3]={0, 0, 0};// no intersection of side red with blue
01027   //double area= 0.;
01028   // X corresponds to blue, Y to red
01029   nP=0; // number of intersection points
01030   int ret = edgeIntersections(blueTriangle, redTriangle, mark, P, nP);
01031   if (ret!=0) ERROR("Something went wrong");
01032     
01033   int extraPoints = borderPointsOfXinY(blueTriangle, redTriangle, &(P[2*nP]));
01034   if (extraPoints>1)
01035     {
01036       mark[0] = mark[1] = mark[2]=1;
01037     }                                                
01038   nP+=extraPoints;
01039   extraPoints =  borderPointsOfXinY(redTriangle, blueTriangle, &(P[2*nP]));
01040   nP+=extraPoints;
01041     
01042   // now sort and orient the points in P, such that they are forming a convex polygon
01043   // this will be the foundation of our new mesh
01044   //
01045   SortAndRemoveDoubles (P, nP); // nP should be at most 6 in the end 
01046   // if there are more than 3 points, some area will be positive
01047   area = 0.;
01048   if (nP>=3)
01049     {
01050       for (int k=1; k<nP-1; k++)
01051         area += area2D(Vector<2>(P), Vector<2>(&P[2*k]), Vector<2>(&P[2*k+2]));
01052     }
01053   return 0; // no error
01054 }
01055 
01056 int ProjectShell::findNodes(int red, int blue, double * iP, int nP)
01057 {
01058   // first of all, check against red and blue vertices
01059   //
01060   if (m_dbg)
01061     {
01062       std::cout<< "red, blue, nP, P " << red << " " << blue << " " << nP <<"\n";
01063       for (int n=0; n<nP; n++)
01064         std::cout << " \t" << iP[2*n] << "\t" << iP[2*n+1] << "\n";
01065         
01066     }
01067   int * foundIds = new int [nP];
01068   for (int i=0; i<nP; i++)
01069     {
01070       double * pp = &iP[2*i];// iP+2*i
01071       int found = 0; 
01072       // first, are they on vertices from red or blue?
01073       ProjectShell::Triangle2D & redTri = m_redMesh[red];
01074       int j=0;
01075       for (j=0; j<3&& !found; j++)
01076         {
01077           int node = redTri.v[j];
01078           double d2 = length(Vector<2>(pp)-Vector<2>(m_xy+(3*node)));   
01079           if (m_dbg && i==0)
01080             std::cout<< "  red node " << j << " " << node << " " << m_xy[3*node] << " " << m_xy[3*node+1] << " d2:" << d2 <<  " \n";
01081           if (d2<epsilon)
01082             {
01083               foundIds[i] = node; // no new node
01084               found = 1;
01085             }
01086         }
01087       ProjectShell::Triangle2D & blueTri = m_blueMesh[blue];
01088       for (j=0; j<3 && !found; j++)
01089         {
01090           int node = blueTri.v[j];
01091           double d2 = length(Vector<2>(pp)-Vector<2>(m_xy+(3*node)));   
01092           if (m_dbg && i==0) 
01093             std::cout<< "  blue node " << j << " " << node << " " << m_xy[3*node] << " " << m_xy[3*node+1] << " d2:" << d2 <<  " \n";
01094           if (d2<epsilon)
01095             {
01096               foundIds[i] = node; // no new node
01097               found = 1;
01098             }
01099         }
01100       if (!found)
01101         {
01102           // find the edge it belongs, first
01103           //
01104           for (j=0; j<3; j++)
01105             {
01106               int edge = redTri.e[j];
01107               int ae = abs(edge);
01108               int v1 = ps_edges[ae-1].v[0];
01109               int v2 = ps_edges[ae-1].v[1];
01110               double area = area2D (Vector<2>(&m_xy[3*v1]), 
01111                                     Vector<2>(&m_xy[3*v2]), 
01112                                     Vector<2>(pp));
01113               if (m_dbg)
01114                 std::cout << "   edge " << j << ": " << edge << " " << v1 << " " << v2 << "  area : " << area << "\n"; 
01115               if ( fabs(area) < epsilon*epsilon  )
01116                 {
01117                   // found the edge; now find if there is a point in the list here
01118                   std::list<int> & expts = ps_edges[ae-1].extraNodes;
01119                   // if the points pp is between extra points, then just give that id
01120                   // if not, create a new point, (check the id)
01121                   //
01122                   std::list<int>::iterator it;
01123                   for ( it = expts.begin(); it!=expts.end() && !found; it++)
01124                     {
01125                       int pnt = *it;
01126                       double d2 = length(Vector<2>(pp)-Vector<2>(&m_xy[3*pnt]));
01127                       if (d2<epsilon)
01128                         {
01129                           found = 1;
01130                           foundIds[i] = pnt;
01131                         }
01132                     }
01133                   if (!found)
01134                     {
01135                       // create a new point in 2d (at the intersection)
01136                       foundIds [i] = m_num2dPoints;
01137                       expts.push_back(m_num2dPoints);
01138                       if (m_2dcapacity < m_num2dPoints+1)
01139                         {
01140                           // Underestimated the number of intersection points
01141                           //std::cout << " underestimated capacity for 2d array\n";
01142                           // double the capacity of m_xy array
01143                             
01144                           double * new_xy = new double [6* m_2dcapacity];
01145                           int jj=0;
01146                           for (jj=0; jj<3*m_2dcapacity; jj++)
01147                             {
01148                               new_xy[jj] = m_xy[jj];
01149                                 
01150                             }
01151                           for (jj=3*m_2dcapacity-1; jj<6*m_2dcapacity; jj+=3)
01152                             new_xy[jj] = 0.;
01153                           // make 0 the z coordinate
01154                           m_2dcapacity *= 2;
01155                           delete [] m_xy;
01156                           m_xy = new_xy;
01157                         }
01158                       m_xy[3*m_num2dPoints] = pp[0];
01159                       m_xy[3*m_num2dPoints+1] = pp[1];
01160                       m_num2dPoints++;
01161                       if (m_dbg)
01162                         {
01163                           std::cout<< " new 2d " << m_num2dPoints - 1 << " : " << pp[0] << " " << pp[1] <<  "on edge " << ae << "\n";
01164                         }
01165                       found = 1;
01166                     }
01167                 } 
01168             }
01169         }
01170       if (!found)
01171         {
01172           ERROR3("Point (%g,%g) is not on a red triangle %d", pp[0], pp[1], red);
01173         }
01174     }
01175   // now we can build the triangles, from P array, with foundIds
01176   if (nP>=3)
01177     {
01178       // what are the triangles 
01179       FinalTriangle ftr;
01180       ftr.redTriangle = red;
01181       ftr.blueTriangle = blue;
01182       ftr.v[0] = foundIds[0];
01183       for (int i=1; i<=nP-2; i++)
01184         {
01185           // triangle 0, i, i+1
01186           ftr.v[1]=foundIds[i];
01187           ftr.v[2] = foundIds[i+1];
01188           m_finalMesh.push_back(ftr);
01189           if (m_dbg)
01190             {
01191               std::cout << " triangle " << ftr.v[0] << " " << ftr.v[1] << " " << ftr.v[2] << "\n"; 
01192             }
01193         }
01194     }
01195   delete [] foundIds;
01196   foundIds = NULL;
01197   return 0;
01198 }
01199   
01200 #if 0  
01201 int ProjectShell::commitProjectedMesh()
01202 {
01203   // here we will create new vertices, and use the coordinates in 2D
01204   //  m_num2dPoints is the number of nodes (some are not used, because they were probably   // collapsed
01205   //  we do not specifically merge red and blue nodes, but we are looking first for 
01206   //  nodes in red , then blue, then on red edges; some will not appear, according
01207   //  to the tolerance
01208   //
01209   //iMesh_Instance mesh = this->mk_core()->imesh_instance();
01210   iMesh_Instance mesh = NULL;
01211   iBase_EntityHandle * newVerts = NULL; // no vertices yet
01212   // iBase_INTERLEAVED
01213   int err= 0;
01214   int size1, size2;
01215   iMesh_createVtxArr(mesh,
01216                      /*in*/ m_num2dPoints,
01217                      /*in*/ iBase_INTERLEAVED,
01218                      /*in*/ m_xy,
01219                      /*in*/ m_num2dPoints*3,
01220                      /*inout*/ &newVerts,
01221                      /*inout*/ &size1,
01222                      /*inout*/ &size2,
01223                      /*out*/ &err); 
01224   if (err!=0)
01225     {
01226       std::cout<<"can't create vertices\n";
01227       exit (1);
01228     }
01229   // size of 
01230   //  then entity set, then elements (triangles)
01231   //
01232   int numTriangles = m_finalMesh.size();
01233   long int * adjacency = new long int [3*numTriangles];
01234   iBase_EntityHandle * conn = (iBase_EntityHandle *)adjacency;
01235   for (int L =0; L<numTriangles; L++)
01236     {
01237       FinalTriangle & tria = m_finalMesh[L];
01238         
01239         
01240       for (int k=0; k<3; k++)
01241         {
01242           int indexInV =  tria.v[k];
01243           conn[3*L+k] = newVerts[indexInV];
01244         }
01245     }
01246   int numElements = numTriangles;
01247   iBase_EntitySetHandle orig_set;
01248   iMesh_createEntSet(mesh, 0, &orig_set, &err);
01249   int n = numTriangles;
01250   int junk1 = n, junk2 = n, junk3 = n, junk4 = n;
01251   int * stat = new int [numElements];
01252   int* ptr2 = stat;
01253   int ierr;
01254   iBase_EntityHandle *           new_entity_handles = NULL;
01255   iMesh_createEntArr( mesh,
01256                       iMesh_TRIANGLE,
01257                       conn, 3*numElements,
01258                       &new_entity_handles, &junk1, &junk2,
01259                       &ptr2, &junk3, &junk4,
01260                       &ierr );
01261   if (ierr!=0)
01262     {
01263       std::cout<<" can't create triangles\n";
01264       exit (1);
01265     }
01266   iMesh_addEntArrToSet          (  mesh,
01267                                    new_entity_handles,
01268                                    numElements,
01269                                    orig_set,
01270                                    &ierr);
01271   if (ierr!=0)
01272     {
01273       std::cout<< " can't add to entity set \n";
01274       exit(1);
01275     }    
01276   // now, look at the tags : red is positive, blue is negative oriented triangle
01277   iBase_TagHandle       pos_tag_handle; 
01278     
01279   const char * tagName1 = "Positive";
01280   iMesh_createTag       ( mesh,
01281                           tagName1,
01282                           /*  size ? */ 1,
01283                           iBase_INTEGER ,
01284                           &pos_tag_handle,
01285                           &ierr,
01286                           strlen(tagName1) ) ;   
01287     
01288   if (ierr!=0)
01289     {
01290       std::cout<< " can't create positive tag \n";
01291       exit(1);
01292     }
01293   int * tagValues = new int [numElements];
01294   for (int e=0; e<numElements; e++)
01295     {
01296       int redId = m_finalMesh[e].redTriangle;
01297       tagValues[e] = m_redMesh[ redId ].oldId;
01298         
01299     }
01300   // positive or negative triangles
01301   iMesh_setIntArrData   (  mesh,
01302                            new_entity_handles,
01303                            numElements,
01304                            pos_tag_handle,
01305                            tagValues,
01306                            numElements,
01307                            &ierr); 
01308   if (ierr!=0)
01309     {
01310       std::cout<< " can't create positive tag field \n";
01311       exit(1);
01312     }
01313   // now blue tag
01314   // now, look at the tags : red is positive, blue is negative oriented triangle
01315   iBase_TagHandle       neg_tag_handle; 
01316     
01317   const char * tagName2 = "Negative";
01318   iMesh_createTag       ( mesh,
01319                           tagName2,
01320                           /*  size ? */ 1,
01321                           iBase_INTEGER ,
01322                           &neg_tag_handle,
01323                           &ierr,
01324                           strlen(tagName1) ) ;   
01325     
01326   if (ierr!=0)
01327     {
01328       std::cout<< " can't create negative tag \n";
01329       exit(1);
01330     }
01331   for (int e=0; e<numElements; e++)
01332     {
01333       int blueId = m_finalMesh[e].blueTriangle;
01334       tagValues[e] = m_blueMesh[ blueId ].oldId;
01335         
01336     }
01337   // positive or negative triangles
01338   iMesh_setIntArrData   (  mesh,
01339                            new_entity_handles,
01340                            numElements,
01341                            neg_tag_handle,
01342                            tagValues,
01343                            numElements,
01344                            &ierr); 
01345   if (ierr!=0)
01346     {
01347       std::cout<< " can't create negative tag field \n";
01348       exit(1);
01349     }
01350   delete [] tagValues;
01351   delete [] adjacency;
01352   return 0;
01353 }
01354 #endif
01355 int edgeIntersections(double * red, double * blue, int mark[3], double * points, int & nPoints)
01356 {
01357   /* EDGEINTERSECTIONS computes edge intersections of two triangles
01358      [P,n]=EdgeIntersections(X,Y) computes for the two given triangles  * red 
01359      and blue ( stored column wise )
01360      (point coordinates are stored column-wise, in counter clock
01361      order) the points P where their edges intersect. In addition,
01362      in n the indices of which neighbors of red  are also intersecting
01363      with blue are given.
01364   */
01365     
01366   // points is an array with 12 slots   (12 * 2 doubles) 
01367   nPoints = 0; 
01368   mark[0] = mark[1] = mark[2]=0 ; // no neighbors of red involved yet
01369   /*for i=1:3                            % find all intersections of edges
01370     for j=1:3
01371     b=Y(:,j)-X(:,i);               
01372     A=[X(:,mod(i,3)+1)-X(:,i) -Y(:,mod(j,3)+1)+Y(:,j)];
01373     if rank(A)==2                   % edges not parallel
01374     r=A\b;
01375     if r(1)>=0 & r(1)<=1 & r(2)>=0 & r(2)<=1,  % intersection found
01376     k=k+1; P(:,k)=X(:,i)+r(1)*(X(:,mod(i,3)+1)-X(:,i)); n(i)=1;
01377     end;
01378     end; 
01379     end;
01380     end;*/
01381   for (int i=0; i<3; i++)
01382     {
01383       for (int j=0; j<3; j++)
01384         {
01385           double b[2];
01386           double a[2][2]; // 2*2
01387           int iPlus1 = (i+1)%3;
01388           int jPlus1 = (j+1)%3;
01389           for (int k=0; k<2; k++)
01390             {
01391               b[k] = blue[2*j+k] - red[2*i+k];
01392               // row k of a: a(k, 0), a(k, 1)
01393               a[ k ][0]    = red  [ 2*iPlus1 + k] - red [2*i + k];
01394               a[ k ][1]    = blue [ 2 * j + k] - blue [2*jPlus1 +k];
01395                 
01396             }
01397           double delta = a[0][0]*a[1][1] - a[0][1]*a[1][0];
01398           if (delta != 0.)
01399             {
01400               // not parallel
01401               double alfa = (b[0]*a[1][1]-a[0][1]*b[1])/delta;
01402               double beta = (-b[0] * a[1][0] + b[1] * a[0][0])/delta;
01403               if (0<= alfa && alfa <=1. && 0<= beta && beta <= 1.)
01404                 {
01405                   // the intersection is good
01406                   for (int k=0; k<2; k++)
01407                     {
01408                       points[2*nPoints+k] = red[2*i+k] + alfa*(red[2*iPlus1+k]-red[2*i+k]);
01409                     }
01410                   mark[i] = 1; // so neighbor number i will be considered too.
01411                   nPoints++;
01412                 }
01413             }
01414             
01415         }
01416     }
01417   return 0;
01418 }
01419 
01420 }//namespace MeshKit
01421 
01422 
01423 
01424 
01425 
01426 
01427 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines