00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef synaps_topology_TopSbdCurve3d_h
00022 #define synaps_topology_TopSbdCurve3d_h
00023
00024
00025 #include <vector>
00026 #include <list>
00027 #include <math.h>
00028 #include <synaps/upol/bezier/BEZIER.m>
00029 #include <synaps/mpol/MPol.h>
00030 #include <synaps/mpol/MPolDse.h>
00031 #include <synaps/msolve/sbdslv.h>
00032 #include <synaps/topology/Cell3dN.h>
00033 #include <synaps/topology/Octree.h>
00034
00035 __BEGIN_NAMESPACE_SYNAPS
00036
00037 template<class C,
00038 class CELL= Cell3dN<5,C>,
00039 class MSLV= SbdSolver<typename CELL::coeff_t, SBDSLV_RDL> >
00040 struct TopSbdCurve3d
00041 {
00042
00043
00044 typedef TopSbdCurve3d<C, CELL,MSLV> self_t;
00045
00046 typedef typename CELL::coeff_t coeff_t;
00047
00048 typedef typename CELL::MPOL_T MPOL_T;
00049 typedef MPolDse<C, mpoldse::bernstein<C> > mpolb_t;
00050
00051 typedef topology::box3d<coeff_t> box_t;
00052 typedef topology::point<C> point_t;
00053 typedef topology::point_graph<C> pointgraph_t;
00054
00055
00056 typedef CELL cell_t;
00057 typedef OCTREE::Node<cell_t*> node_t;
00058 typedef OCTREE::Octree<cell_t*> octree_t;
00059
00060 typedef intersection3d<C> intersection_t;
00061 typedef shared_object<intersection_t*> intersection_ptr;
00062
00063 typedef Seq<intersection3d<coeff_t>*> intersectionlist_t;
00064
00065 typedef Seq<intersection_ptr> intersection_ptr_list;
00066 typedef topology::point_graph<coeff_t, intersection_ptr_list> pointgraph_int_t;
00067
00068 protected:
00069
00070
00071 C m_dEpsilon;
00072 int m_nMinSbd;
00073 octree_t m_cOcTree;
00074
00075 int m_nCount;
00076
00077 public:
00078
00079
00080 TopSbdCurve3d():m_dEpsilon(1e-3), m_nMinSbd(0),m_nCount(0) {}
00081 TopSbdCurve3d(coeff_t e):m_dEpsilon(e), m_nMinSbd(0),m_nCount(0) {}
00082 TopSbdCurve3d(coeff_t e, int d):m_dEpsilon(e), m_nMinSbd(d),m_nCount(0) {}
00083
00084 ~TopSbdCurve3d();
00085
00086
00087
00088
00089 void computeinter_face(intersectionlist_t& inters, cell_t& c, int nAxis, int nFace);
00090
00091 void computeintersections(intersectionlist_t& inters, cell_t& c);
00092
00093 template<class POINTGRAPH>
00094 void run(POINTGRAPH& g, CELL* rc);
00095
00096 protected:
00097
00098
00099
00100 void subdivide(node_t* n, std::vector<node_t*>& aNodeList);
00101
00102 void subdivide(CELL* cell, std::vector<CELL*>& subCells);
00103
00104
00105 void deleteNode(node_t* n, std::vector<intersection_t*>& inters);
00106 };
00107
00109
00111
00112
00113
00115
00117 template< class X, class CELL, class MSLV>
00118 TopSbdCurve3d<X, CELL, MSLV>::~TopSbdCurve3d()
00119 {
00120
00121 m_nCount++;
00122 std::vector<intersection_t*> inters;
00123 deleteNode(m_cOcTree.GetRoot(), inters);
00124 for (int i = 0; i < inters.size(); i++)
00125 delete inters[i];
00126 }
00127
00128
00130
00131
00132
00133
00134
00135
00136
00137
00138
00139
00140
00141
00142
00144
00145 template< class C, class CELL, class MSLV>
00146 void TopSbdCurve3d<C, CELL, MSLV>::deleteNode(node_t* n, std::vector<intersection_t*>& inters)
00147 {
00148 if (n->GetF_NEchild())
00149 deleteNode(n->F_NEchild, inters);
00150 if (n->GetF_NWchild())
00151 deleteNode(n->F_NWchild, inters);
00152 if (n->GetF_SEchild())
00153 deleteNode(n->F_SEchild, inters);
00154 if (n->GetF_SWchild())
00155 deleteNode(n->F_SWchild, inters);
00156 if (n->GetB_NEchild())
00157 deleteNode(n->B_NEchild, inters);
00158 if (n->GetB_NWchild())
00159 deleteNode(n->B_NWchild, inters);
00160 if (n->GetB_SEchild())
00161 deleteNode(n->B_SEchild, inters);
00162 if (n->GetB_SWchild())
00163 deleteNode(n->B_SWchild, inters);
00164
00165 for (int i = 0; i < n->GetCell()->m_aInterList.size(); i++)
00166 {
00167 if (n->GetCell()->m_aInterList[i]->nCount != m_nCount)
00168 {
00169 inters.push_back(n->GetCell()->m_aInterList[i]);
00170 n->GetCell()->m_aInterList[i]->nCount=m_nCount;
00171 }
00172 }
00173
00174 delete n->GetCell();
00175 delete n;
00176 }
00177
00178 template<class C, class CELL, class MSLV>
00179 void TopSbdCurve3d<C, CELL, MSLV>::subdivide(node_t* n, std::vector<node_t*>& aNodeList)
00180 {
00181 using namespace OCTREE;
00182 cell_t* c = n->GetCell();
00183 std::vector<cell_t*> aSubCells;
00184
00185
00186 subdivide(c,aSubCells);
00187
00188
00189 aNodeList.push_back(new node_t(B_SW, n, aSubCells[0]));
00190 aNodeList.push_back(new node_t(F_SW, n, aSubCells[1]));
00191 aNodeList.push_back(new node_t(B_SE, n, aSubCells[2]));
00192 aNodeList.push_back(new node_t(F_SE, n, aSubCells[3]));
00193 aNodeList.push_back(new node_t(B_NW, n, aSubCells[4]));
00194 aNodeList.push_back(new node_t(F_NW, n, aSubCells[5]));
00195 aNodeList.push_back(new node_t(B_NE, n, aSubCells[6]));
00196 aNodeList.push_back(new node_t(F_NE, n, aSubCells[7]));
00197 }
00198
00199
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00218
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00236
00237 template < class C, class CELL, class MSLV >
00238 void TopSbdCurve3d<C, CELL, MSLV>::subdivide(CELL* cell,
00239 std::vector<CELL*>& subCells)
00240 {
00241 typedef typename CELL::intersection_t intersection_t;
00242 typedef typename CELL::intersectionlist_t intersectionlist_t;
00243 typedef typename CELL::MPOL_T MPOL_T;
00244
00245 #ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
00246
00247 C t[] = {(C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
00248 (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5),
00249 (C(rand())/C(RAND_MAX)-C(0.5))*C(0.095)+C(0.5)};
00250 t[0]+= (t[0]>0.5?0.005:-0.005);
00251 t[1]+= (t[1]>0.5?0.005:-0.005);
00252 t[2]+= (t[2]>0.5?0.005:-0.005);
00253 C mid[] = {cell->m_cBox.x0*(C(1)-t[0])+cell->m_cBox.x1*t[0], cell->m_cBox.y0*(C(1)-t[1])+cell->m_cBox.y1*t[1], cell->m_cBox.z0*(C(1)-t[2])+cell->m_cBox.z1*t[2]};
00254 #else
00255 C mid[] = {(cell->m_cBox.x0+cell->m_cBox.x1)/C(2), (cell->m_cBox.y0+cell->m_cBox.y1)/C(2), (cell->m_cBox.z0+cell->m_cBox.z1)/C(2)};
00256 #endif
00257 C box[] = {cell->m_cBox.x0, cell->m_cBox.x1, cell->m_cBox.y0, cell->m_cBox.y1, cell->m_cBox.z0, cell->m_cBox.z1};
00258
00259
00260 subCells.push_back(new CELL(box[0], mid[0], box[2], mid[1], box[4], mid[2]));
00261 subCells.push_back(new CELL(mid[0], box[1], box[2], mid[1], box[4], mid[2]));
00262 subCells.push_back(new CELL(box[0], mid[0], mid[1], box[3], box[4], mid[2]));
00263 subCells.push_back(new CELL(mid[0], box[1], mid[1], box[3], box[4], mid[2]));
00264 subCells.push_back(new CELL(box[0], mid[0], box[2], mid[1], mid[2], box[5]));
00265 subCells.push_back(new CELL(mid[0], box[1], box[2], mid[1], mid[2], box[5]));
00266 subCells.push_back(new CELL(box[0], mid[0], mid[1], box[3], mid[2], box[5]));
00267 subCells.push_back(new CELL(mid[0], box[1], mid[1], box[3], mid[2], box[5]));
00268
00269
00270 for (int i = 0; i < 8; i++)
00271 {
00272 subCells[i]->m_nDepth = cell->m_nDepth+1;
00273 }
00274
00275
00276 for (int i = 0; i < cell->m_aInterList.size(); i++)
00277 {
00278 intersection_t* inter = cell->m_aInterList[i];
00279
00280
00281
00282 int nidx = inter->nPlane;
00283 int *cidx=subcellidx+(nidx<<3);
00284 int axisidx0 = plane_axisidx[nidx<<1], axisidx1 = plane_axisidx[(nidx<<1)+1];
00285
00286 assert(_is_verysmall(inter->pnt[nidx]-box[nidx<<1]) || _is_verysmall(box[(nidx<<1)+1]-inter->pnt[nidx]));
00287 assert(inter->pnt[axisidx0] >= box[axisidx0<<1] && inter->pnt[axisidx0] <= box[(axisidx0<<1)+1] &&
00288 inter->pnt[axisidx1] >= box[axisidx1<<1] && inter->pnt[axisidx1] <= box[(axisidx1<<1)+1]);
00289
00290 if (inter->pnt[axisidx0] < mid[axisidx0])
00291 {
00292 if (inter->pnt[axisidx1] < mid[axisidx1])
00293 {
00294 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx])
00295 subCells[cidx[0]]->m_aInterList.push_back(inter);
00296 else
00297 subCells[cidx[1]]->m_aInterList.push_back(inter);
00298 }
00299 else
00300 {
00301 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx])
00302 subCells[cidx[2]]->m_aInterList.push_back(inter);
00303 else
00304 subCells[cidx[3]]->m_aInterList.push_back(inter);
00305 }
00306 }
00307 else
00308 {
00309 if (inter->pnt[axisidx1] < mid[axisidx1])
00310 {
00311 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx])
00312 subCells[cidx[4]]->m_aInterList.push_back(inter);
00313 else
00314 subCells[cidx[5]]->m_aInterList.push_back(inter);
00315 }
00316 else
00317 {
00318 if (inter->pnt[nidx]-box[nidx<<1] < box[(nidx<<1)+1]-inter->pnt[nidx])
00319 subCells[cidx[6]]->m_aInterList.push_back(inter);
00320 else
00321 subCells[cidx[7]]->m_aInterList.push_back(inter);
00322 }
00323 }
00324 }
00325
00326
00327 for (int i = 0; i < 8*MAX_FACEIDX; i++)
00328 {
00329 int nCellIdx=_cellidx[i], nInheritIdx = _inheritidx[i];
00330
00331 assert((subCells[nCellIdx]->m_aFaceComputed[nInheritIdx]) == false);
00332
00333
00334 if (cell->m_aFaceComputed[nInheritIdx])
00335 subCells[nCellIdx]->m_aFaceComputed[nInheritIdx] = true;
00336 }
00337
00338
00339
00340 for (int i = 0; i < 5; i++)
00341 {
00342 MPOL_T splited[8];
00343 int vr_idx[3] = {VARIDX(cell->m_F[i], 0), VARIDX(cell->m_F[i], 1), VARIDX(cell->m_F[i], 2)};
00344
00345 subCells[7]->m_F[i] = cell->m_F[i];
00346 #ifdef SYNAPS_SBD3D_RNDFACTORWHENSUBDIVISION
00347
00348 mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2], t[2]);
00349
00350
00351 mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1], t[1]);
00352 mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1], t[1]);
00353
00354
00355 mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0], t[0]);
00356 mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0], t[0]);
00357 mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0], t[0]);
00358 mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0], t[0]);
00359 #else
00360
00361 mpoldse::split(subCells[3]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[2]);
00362
00363
00364 mpoldse::split(subCells[1]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[1]);
00365 mpoldse::split(subCells[5]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[1]);
00366
00367
00368 mpoldse::split(subCells[0]->m_F[i].rep(), subCells[1]->m_F[i].rep(), vr_idx[0]);
00369 mpoldse::split(subCells[2]->m_F[i].rep(), subCells[3]->m_F[i].rep(), vr_idx[0]);
00370 mpoldse::split(subCells[4]->m_F[i].rep(), subCells[5]->m_F[i].rep(), vr_idx[0]);
00371 mpoldse::split(subCells[6]->m_F[i].rep(), subCells[7]->m_F[i].rep(), vr_idx[0]);
00372 #endif
00373 }
00374
00375
00376
00377 for (int i = 0; i < 8*MAX_FACEIDX; i++)
00378 {
00379 intersectionlist_t aInters;
00380 int nCellIdx=_cellidx[i], nFaceIdx = _faceidx[i], nShareIdx = _shareidx[i];
00381
00382 if (nFaceIdx >= 0)
00383 {
00384 computeinter_face(aInters, *subCells[nCellIdx], nFaceIdx>>1, nFaceIdx%2);
00385
00386
00387 cell->mergelist(subCells[nCellIdx]->m_aInterList, aInters);
00388 cell->mergelist(subCells[nShareIdx]->m_aInterList, aInters);
00389
00390 subCells[nCellIdx]->m_aFaceComputed[nFaceIdx] = true;
00391
00392
00393
00394
00395
00396
00397 assert(subCells[nShareIdx]->m_aFaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] == false);
00398
00399
00400 subCells[nShareIdx]->m_aFaceComputed[nFaceIdx+(nFaceIdx%2?-1:1)] = true;
00401 }
00402 }
00403 }
00404
00405 template <class C, class CELL, class MSLV >
00406 void TopSbdCurve3d<C, CELL, MSLV>::computeintersections(intersectionlist_t& inters, CELL& c)
00407 {
00408
00409 for (int v = 0; v < 3; v++)
00410 {
00411 for (int f = 0; f < 2; f++)
00412 {
00413
00414 if (!c.m_aFaceComputed[v*2+f])
00415 {
00416 computeinter_face(inters, c, v, f);
00417 c.m_aFaceComputed[v*2+f] = true;
00418 }
00419 }
00420 }
00421 }
00422
00423 template < class C, class CELL, class MSLV >
00424 void TopSbdCurve3d<C, CELL, MSLV>::computeinter_face(intersectionlist_t& inters,
00425 CELL& c,
00426 int nAxis, int nFace)
00427 {
00428 point_t l(3u);
00429
00430
00431 MPOL_T facepol[2];
00432 face(facepol[0].rep(), c.m_F[0].rep(), nAxis, nFace);
00433 face(facepol[1].rep(), c.m_F[1].rep(), nAxis, nFace);
00434
00435
00436 std::list<MPOL_T> eqnList;
00437 eqnList.push_back(facepol[0]);
00438 eqnList.push_back(facepol[1]);
00439
00440
00441 MSLV solver(facepol, facepol+2, 1e-4);
00442 std::vector<C> results;
00443 int nSol, nVar;
00444
00445 solver.solve(results, nSol, nVar, eqnList.begin(), eqnList.end());
00446 assert(nVar==2);
00447
00448
00449
00450 for (int s = 0; s < nSol; s++)
00451 {
00452 C x[3] = {(C)nFace,(C)nFace,(C)nFace};
00453
00454 switch(nAxis)
00455 {
00456 case PLANEYZ:
00457 x[1] = (C)(results[nVar*s]);
00458 x[2] = (C)(results[nVar*s+1]);
00459 break;
00460 case PLANEXZ:
00461 x[0] = (C)(results[nVar*s]);
00462 x[2] = (C)(results[nVar*s+1]);
00463 break;
00464 case PLANEXY:
00465 x[0] = (C)(results[nVar*s]);
00466 x[1] = (C)(results[nVar*s+1]);
00467 break;
00468 }
00469
00470
00471
00472
00473 int nMulti = 1;
00474
00475
00476
00477
00478 C val[3] = {-1,-1,-1};
00479 eval(val[0], c.m_T[0], x);
00480 eval(val[1], c.m_T[1], x);
00481 eval(val[2], c.m_T[2], x);
00482
00483
00484
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 val[nAxis] /= sqrt(val[0]*val[0]+val[1]*val[1]+val[2]*val[2]);
00495
00496 if (_is_verysmall(std::abs(val[nAxis])))
00497 {
00498 nMulti = 2;
00499
00500 }
00501
00502
00503
00504
00505 l[0] = c.m_cBox.x0*((C)1-x[0]) + c.m_cBox.x1*x[0];
00506 l[1] = c.m_cBox.y0*((C)1-x[1]) + c.m_cBox.y1*x[1];
00507 l[2] = c.m_cBox.z0*((C)1-x[2]) + c.m_cBox.z1*x[2];
00508
00509 intersection_t *ni = new intersection_t(l, nAxis, nMulti);
00510 inters.push_back(ni);
00511 ni->nDepth = c.m_nDepth;
00512 }
00513 }
00514
00515
00516 template< class C, class CELL, class MSLV>
00517 template< class POINTGRAPH>
00518 void TopSbdCurve3d<C, CELL, MSLV>::run(POINTGRAPH& g, CELL* rc)
00519 {
00520
00521 node_t* root = m_cOcTree.GetRoot();
00522
00523
00524
00525 computeintersections(rc->m_aInterList,*rc);
00526 std::cout << "no. of roots at root: " << rc->m_aInterList.size() << std::endl;
00527 root->SetCell(rc);
00528
00529
00530
00531 std::list<node_t*> aNodeQueue;
00532 std::vector<node_t*> aSingular;
00533 pointgraph_int_t g_tmp;
00534 int nCellProc = 0;
00535
00536 aNodeQueue.push_front(root);
00537
00538 while(aNodeQueue.size())
00539 {
00540
00541 node_t* n = *aNodeQueue.begin();
00542 aNodeQueue.pop_front();
00543
00544 bool bSubdivide = false;
00545 bool bBuildGraph = false;
00546
00547
00548 if (regularity(*n->GetCell()))
00549 {
00550 if (is_variety_in(*n->GetCell()))
00551 {
00552 if (n->GetCell()->m_nDepth < m_nMinSbd)
00553 {
00554
00555 bSubdivide = true;
00556 }
00557 else
00558 {
00559 bBuildGraph = true;
00560 }
00561 }
00562 else
00563 {
00564
00565 memset(n->GetCell()->m_aFaceComputed, true, sizeof(bool)*6);
00566 }
00567 }
00568 else
00569 {
00570 if (n->GetCell()->testsize(m_dEpsilon))
00571 {
00572
00573 if (is_variety_in(*n->GetCell()))
00574 {
00575 bBuildGraph = true;
00576 aSingular.push_back(n);
00577 }
00578 }
00579 else
00580 {
00581
00582 if (is_variety_in(*n->GetCell()))
00583 bSubdivide = true;
00584 }
00585 }
00586
00587
00588
00589 if (bSubdivide)
00590 {
00591
00592
00593 std::vector<node_t*> aNodeList;
00594 subdivide(n, aNodeList);
00595 for (int i = 0; i < (int)aNodeList.size(); i++)
00596 aNodeQueue.push_front(aNodeList[i]);
00597 }
00598 else if (bBuildGraph)
00599 {
00600 build_graph(g_tmp,*n->GetCell());
00601 }
00602
00603
00604 nCellProc++;
00605 }
00606
00607
00608
00609 int nVert = (int)rep(rep(g_tmp).Vertices).size();
00610 int nEdges = 0;
00611
00612 for(int i=0; i<(int)rep(rep(g_tmp).Edges).size(); ++i)
00613 nEdges += (int)rep(rep(g_tmp).Edges)[i].size();
00614
00615 if (nVert)
00616 {
00617 int nVertIdx = 0;
00618 int nBase = (int)rep(rep(g).Vertices).size();;
00619 m_nCount++;
00620 for(int i=0; i<nVert; ++i)
00621 {
00622 if (rep(rep(rep(g_tmp).Vertices)[i])->nCount != m_nCount)
00623 {
00624 rep(rep(rep(g_tmp).Vertices)[i])->nCount = m_nCount;
00625 rep(rep(rep(g_tmp).Vertices)[i])->nFlag = nVertIdx++;
00626 g.push_back(rep(rep(rep(g_tmp).Vertices)[i])->pnt);
00627 }
00628 }
00629
00630 for(int i=0; i<(int)rep(rep(g_tmp).Edges).size(); ++i)
00631 {
00632 for(int j=0; j<(int)rep(rep(g_tmp).Edges)[i].size(); ++j)
00633 {
00634 g.insert(rep(rep(rep(g_tmp).Vertices)[i])->nFlag+nBase,
00635 rep(rep(rep(g_tmp).Vertices)[rep(rep(g_tmp).Edges)[i][j]])->nFlag+nBase);
00636 }
00637 }
00638
00639 std::cout << "No. of Vertices: " << (int)rep(rep(g).Vertices).size() - nBase << std::endl;
00640 std::cout << "No. of Edges: " << nEdges << std::endl;
00641 }
00642
00643 std::cout << "Total Cell Processed: " << nCellProc << std::endl;
00644 std::cout << "Sing Cells: " << aSingular.size() << std::endl;
00645 }
00646
00647
00648
00649 namespace topology {
00650
00651 template<class C, class U, class CELL, class SLV>
00652 void assign(point_graph<C>& g,
00653 const MPol<U>& p1,
00654 const MPol<U>& p2,
00655 TopSbdCurve3d<C,CELL,SLV> slv,
00656 double x0, double x1,
00657 double y0, double y1,
00658 double z0, double z1)
00659 {
00660 typedef TopSbdCurve3d<C,CELL,SLV> MTH;
00661 typedef typename MTH::cell_t cell_t;
00662 typedef typename MTH::mpolb_t mpolb_t;
00663
00664
00665
00666 mpolb_t p1dse, p2dse;
00667
00668 toBernstein(p1, p1dse, x0,x1,y0,y1,z0,z1);
00669 toBernstein(p2, p2dse, x0,x1,y0,y1,z0,z1);
00670
00671
00672
00673
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696 mpolb_t p1_dx, p1_dy, p1_dz;
00697 mpolb_t p2_dx, p2_dy, p2_dz;
00698
00699 diff(p1_dx, p1dse, VARIDX(pd1dse, 0));
00700 diff(p1_dy, p1dse, VARIDX(pd1dse, 1));
00701 diff(p1_dz, p1dse, VARIDX(pd1dse, 2));
00702
00703 diff(p2_dx, p2dse, VARIDX(pd1dse, 0));
00704 diff(p2_dy, p2dse, VARIDX(pd2dse, 1));
00705 diff(p2_dz, p2dse, VARIDX(pd2dse, 2));
00706
00707
00708
00709
00710
00711
00712
00713
00714
00715 mpolb_t c1, c2, c3;
00716 using arithm::operator*;
00717 c1 = p1_dy * p2_dz - p1_dz * p2_dy;
00718 c2 = p2_dx * p1_dz - p2_dz * p1_dx;
00719 c3 = p1_dx * p2_dy - p1_dy * p2_dx;
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729 cell_t* rc = new cell_t(x0,x1,y0,y1,z0,z1);
00730 rc->m_F[0] = p1dse;
00731 rc->m_F[1] = p2dse;
00732 rc->m_T[0] = c1;
00733 rc->m_T[1] = c2;
00734 rc->m_T[2] = c3;
00735
00736 slv.run(g,rc);
00737 }
00738 }
00739
00740
00741 template < class C >
00742 bool is_variety_in(Cell3dN<5,C>& c)
00743 {
00744 return c.containsVariety(0) && c.containsVariety(1);
00745 }
00746
00747
00748 template < class POINTGRAPH, class CELL>
00749 void build_graph(POINTGRAPH& g, CELL& c)
00750 {
00751
00752 typedef typename CELL::coeff_t coeff_t;
00753 typedef typename CELL::MPOL_T MPOL_T;
00754 typedef intersection3d<coeff_t> intersection_t;
00755 typedef topology::point<coeff_t> point_t;
00756 typedef shared_object<intersection_t*> intersection_ptr;
00757
00758 if (c.m_nRegularIdx == -2)
00759 regularity(c);
00760
00761 if (c.m_nRegularIdx == -1)
00762 {
00763 if (c.m_aInterList.size())
00764 {
00765 buildsingular:
00766 ;
00767 }
00768 }
00769 else
00770 {
00771 Seq<intersection_t*> inters = c.m_aInterList;
00772
00773
00774
00775
00776
00777 switch(c.m_nRegularIdx)
00778 {
00779 case 0:
00780 sort(inters.begin(), inters.end(), XMAJOR());
00781 break;
00782 case 1:
00783 sort(inters.begin(), inters.end(), YMAJOR());
00784 break;
00785 case 2:
00786 sort(inters.begin(), inters.end(), ZMAJOR());
00787 break;
00788 }
00789
00790
00791 bool bConnectWithLast = false;
00792 int nMulti = 1;
00793 int nLastIdx = -1;
00794 for (int v = 0; v < (int)inters.size(); v++)
00795 {
00796 int nBIdx = g.push_back(intersection_ptr(inters[v]));
00797
00798 if (bConnectWithLast)
00799 g.insert(nLastIdx, nBIdx);
00800
00801 if (!bConnectWithLast && nMulti < inters[v]->nMulti)
00802 bConnectWithLast = false;
00803 else
00804 bConnectWithLast = !bConnectWithLast;
00805
00806 if (nMulti < inters[v]->nMulti)
00807 {
00808 nMulti++;
00809 v--;
00810 }
00811 else
00812 {
00813 nMulti = 1;
00814 }
00815
00816 nLastIdx = nBIdx;
00817 }
00818
00819 if (inters.size() > 2)
00820 {
00821
00822
00823
00824
00825
00826
00827 goto buildsingular;
00828 }
00829 }
00830 }
00831
00832 template < class C >
00833 bool regularity(Cell3dN<5,C>& c)
00834 {
00835
00836
00837 for (int i = 0; i < 3; i++)
00838 {
00839 C* coeff = c.m_T[i].begin();
00840 C lastC = coeff[0];
00841 bool bSignChange = false;
00842
00843 for (int j = 1; j < c.m_T[i].rep().env.sz(); j++)
00844 {
00845 C cirC = coeff[j];
00846 if (_is_singular(cirC))
00847 {
00848 continue;
00849 }
00850 else
00851 {
00852 if (lastC*cirC < 0)
00853 {
00854 bSignChange = true;
00855 lastC = cirC;
00856 break;
00857 }
00858
00859 lastC = cirC;
00860 }
00861 }
00862
00863 if (_is_singular(lastC))
00864 continue;
00865 else if (!bSignChange)
00866 {
00867 c.m_nRegularIdx = i;
00868 return TRUE;
00869 }
00870 }
00871 c.m_nRegularIdx = -1;
00872 return FALSE;
00873 }
00874
00875 template<typename C, typename MPOL, typename MPOLB>
00876 void toBernstein(const MPOL& p, MPOLB& bp,
00877 const C& a0, const C& b0,
00878 const C& a1, const C& b1,
00879 const C& a2, const C& b2)
00880 {
00881 MPOLB tmp(p);
00882
00883 mpoldse::restrict(tmp.rep(),VARIDX(tmp,0),a0,b0);
00884 mpoldse::restrict(tmp.rep(),VARIDX(tmp,1),a1,b1);
00885 mpoldse::restrict(tmp.rep(),VARIDX(tmp,2),a2,b2);
00886
00887 bp = tmp;
00888 }
00889
00890 template < class POINTGRAPH, class CELL>
00891 void build_box(POINTGRAPH& g, CELL* c)
00892 {
00893 typedef typename CELL::coeff_t coeff_t;
00894 typedef topology::point<coeff_t> point_t;
00895 int nBIdx;
00896 point_t l(3u);
00897
00898
00899 for (int i = 0; i < 8; i++)
00900 {
00901 l[0] = ((i&0x1) ? c->m_cBox.x1:c->m_cBox.x0);
00902 l[1] = ((i&0x2) ? c->m_cBox.y1:c->m_cBox.y0);
00903 l[2] = ((i&0x4) ? c->m_cBox.z1:c->m_cBox.z0);
00904
00905
00906 nBIdx = g.push_back(l);
00907 }
00908
00909 nBIdx -= 7;
00910
00911 g.insert(nBIdx,nBIdx+1);
00912 g.insert(nBIdx,nBIdx+2);
00913 g.insert(nBIdx+1,nBIdx+3);
00914 g.insert(nBIdx+2,nBIdx+3);
00915
00916 g.insert(nBIdx+4,nBIdx+5);
00917 g.insert(nBIdx+4,nBIdx+6);
00918 g.insert(nBIdx+5,nBIdx+7);
00919 g.insert(nBIdx+6,nBIdx+7);
00920
00921 g.insert(nBIdx,nBIdx+4);
00922 g.insert(nBIdx+1,nBIdx+5);
00923 g.insert(nBIdx+2,nBIdx+6);
00924 g.insert(nBIdx+3,nBIdx+7);
00925 }
00926
00927 __END_NAMESPACE_SYNAPS
00928
00929 #endif // synaps_topology_TopSbdCurve3d_h
00930